리차드슨 오차 추정
빌드업
미분방정식을 푸는 메소드의 퍼포먼스를 확인하는 방법으로 참값과 비교할 수 있다면 가장 좋겠지만, 당장 참값을 구하기 귀찮은 경우부터 시작해서 아예 트루 솔루션을 구하기 곤란한 경우도 많다. 이땐 $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$ 를 작게 잡을수록 줄어든다.
정의 1
쉽게 말해 리차드슨 오차 추정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)
Atkinson. (1989). An Introduction to Numerical Analysis(2nd Edition): p372. ↩︎