リチャードソン誤差推定
ビルドアップ
微分方程式を解くメソッドのパフォーマンスを確認する方法として、真値と比較するのが最も良いが、すぐに真値を求めるのが面倒な場合から、そもそもトゥルーソリューションを求めるのが難しい場合も多い。こんな時は、$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}}$を使わない理由は全くない。
実装
以下は、リチャードソン誤差推定を通じてステップサイズを変えながら自動的に問題を解くPythonコードだ。
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. ↩︎