📂수치해석

리차드슨 오차 추정

빌드업

미분방정식을 푸는 메소드의 퍼포먼스를 확인하는 방법으로 참값과 비교할 수 있다면 가장 좋겠지만, 당장 참값을 구하기 귀찮은 경우부터 시작해서 아예 트루 솔루션을 구하기 곤란한 경우도 많다. 이땐 $y_{h} (x_{n} )$ 과 스탭사이즈 $h$ 를 두배로 늘렸을 때의 $y_{2h} (x_{n} )$ 을 비교함으로써 오차를 추정할 수 있다.

예를 들어 사다리꼴 메소드라면 $h$ 를 고침으로써 아래의 두 식을 얻는다. $$ Y(x_{n}) - y_{h} (x_{n} ) = D ( x_{n} ) h^2 + O ( h^3 ) $$

$$ Y(x_{n}) - y_{2h} (x_{n} ) = 4 D ( x_{n} ) h^2 + O ( h^3 ) $$ 간단한 조작으로써 $D( x_{n} )$ 를 없애보면 $$ Y(x_{n}) = {{1 } \over {3}} [4 y_{h} (x_{n} ) - y_{2h} (x_{n} ) ] + O ( h^3 ) $$ 우리가 원하는 형태로 정리하면 $$ Y(x_{n}) - y_{h} (x_{n} )= {{1 } \over {3}} [ y_{h} (x_{n} ) - y_{2h} (x_{n} ) ] + O ( h^3 ) $$ $O(h^3)$ 는 꽤 작으므로 버리고나면 $$ Y(x_{n}) - y_{h} (x_{n} ) \approx {{1 } \over {3}} [ y_{h} (x_{n} ) - y_{2h} (x_{n} ) ] $$ 즉, 우리는 실제 오차 $\displaystyle Y(x_{n}) - y_{h} (x_{n} )$ 의 추정치로써 $\displaystyle {{1 } \over {3}} [ y_{h} (x_{n} ) - y_{2h} (x_{n} ) ] $ 를 트루 솔루션 없이 계산할 수 있다. 굳이 말할 것도 없겠지만 이 오차는 $h$ 를 작게 잡을수록 줄어든다.

정의 ¹

쉽게 말해 리차드슨 오차 추정^{richardson Error Estimation}이란, 참값이 없는 상태에서 오차를 예상하는 방법을 말한다.

이러한 계산은 비단 사다리꼴 메소드에만 사용할 게 아니라 수치해석 전반에서 요긴하게 쓰일 수 있다. 또한 위의 계산은 어디까지나 예시기 때문에 간단하게 $h$ 와 $2h$ 로 계산했지만, 실제로 리차드슨 오차 추정을 해야할 상황이 온다면 $h$ 와 $\displaystyle {{h} \over {2}}$ 를 쓰지 않을 이유가 전혀 없다.

구현

아래는 리차드슨 오차 추정을 통해 스탭사이즈를 바꿔가며 오토매틱하게 문제를 푸는 파이썬 코드다.

def Detrap(f,x_0,y_0,x_end,tol,h,h_min,h_max,ier) :
    print("x_n\t\th\t\ty_n\t\t\tY_n-y_n\t\ttrunc")
    x_1 = 0
    x_2 = 0
    y_1 = 0
    y_2 = 0
 
    #2
    loop = 1
    ier = 0
 
    while 1 :
    #for i in range(5) :
        #4
        y_guess = y_0 + h * f(x_0, y_0)
        y_guess = y_0 + (h / 2) \ast (f(x_0, y_0) + f(x_0 + h, y_guess))
        y_e = y_0 + (h / 2) \ast (f(x_0, y_0) + f(x_0 + h, y_guess)) # y_e : y_{h} (x_{0} + (h))
 
        y_guess = y_0 + (h / 2) * f(x_0, y_0)
        y_guess = y_0 + (h / 4) \ast (f(x_0, y_0) + f(x_0 + h / 2, y_guess))
        y_ = y_0 + (h / 4) \ast (f(x_0, y_0) + f(x_0 + h / 2, y_guess))  # y_ : y_{h/2} (x_{0} + (h/2))
 
        y_guess = y_ + (h / 2) * f(x_0 + h / 2, y_)
        y_guess = y_ + (h / 4) \ast (f(x_0 + h / 2, y_) + f(x_0 + h, y_guess))
        y__ = y_ + (h / 4) \ast (f(x_0 + h / 2, y_) + f(x_0 + h, y_guess))  # y__ : y_{h/2} (x_{0} + h)
 
        #5
        trunc = (4/3)*(y__ - y_e)
 
        #6
        if ( tol*h/4 <= abs(trunc) <= tol*h ) or loop==2 :
            x_1 = x_0 + h
            y_1 = y_e
            print("%0.4f\t%0.4f\t%0.6f\t%1.2e\t%1.2e" % (x_1,h,y_1,Y(x_1)-y_1,trunc))
            break
 
        #7
        D_3y = (f(x_0 + h,y__)- 2*f( x_0 + (h/2),y_) + f(x_0,y_0) ) / (h**2 / 4)
        if D_3y==0 :
            h = h_max
            loop = 2
        else :
            h = (6*tol/abs(D_3y))**(1/2)
 
        #8
        if h<h_min :
            ier = 2
            return
 
        if h>h_max :
            h=h_max
            ier = 1
            loop = 2
 
    while 1 :
        linepointer = [0.]
 
        while 1 :
            #11
            x_2 = x_1 + h
            y_guess = y_0 + 2*h*f(x_1,y_1)
            y_2 = y_1 + (h / 2) \ast (f(x_1, y_1) + f(x_2, y_guess))
 
            #12
            trunc = -(y_2 - y_guess)/6
 
            #13
            if abs(trunc)<((tol*h)/4) or (tol*h)<abs(trunc) :
                break
 
            #14
            if round(x_2,4) in [0.2725,1.9595,7.6959] :
                print("x_n\t\th\t\ty_n\t\t\tY_n-y_n\t\ttrunc")
                print("-------------------------------------------------")
            print("%0.4f\t%0.4f\t%0.6f\t%1.2e\t%1.2e" % (x_2,h,y_2,Y(x_2)-y_2,trunc))
 
            #15
            x_1 = x_2
            y_0 = y_1
            y_1 = y_2
 
            if not x_1<x_end :
                return
 
        #17
        x_0 = x_1
        y_0 = y_1
        h_0 = h
        h = ( (tol*( h_0**3 ) ) / ( 2 * abs(trunc)) )**(1/2)
 
        #18
        h = min(h,2*h_0)
 
        #19
        if h<h_min :
            ier = 2
            return
 
        if h>h_max :
            ier = 1
            h = h_max
 
        #20
        y_guess = y_0 + h * f(x_0,y_0)
        y_guess = y_0 + (h / 2) \ast (f(x_0, y_0) + f(x_0 + h, y_guess))
        y_1 = y_0 + (h / 2) \ast (f(x_0, y_0) + f(x_0 + h, y_guess))
        x_1 = x_0 + h
 
        #21
        if round(x_1,4) in [0.2725,1.9595,7.6959] :
            print("x_n\t\th\t\ty_n\t\t\tY_n-y_n\t\ttrunc")
            print("-------------------------------------------------")
 
        print("%0.4f\t%0.4f\t%0.6f\t%1.2e\t%1.2e" % (x_1,h,y_1,Y(x_1)-y_1,trunc))
 
        #22
        if not x_1<x_end :
            return
 
#Example
def f(x,y) :
    return (1/(1 + x**2)) - 2 * ( y**2 )
 
def Y(x) :
    return x/(1 + x**2)
 
Detrap(f,x_0 = 0,y_0 = 0,x_end = 10,tol=0.0005,h=0.1,h_min=0.001,h_max=1.0,ier=0)