📂数値解析

リチャードソン誤差推定

ビルドアップ

微分方程式を解くメソッドのパフォーマンスを確認する方法として、真値と比較するのが最も良いが、すぐに真値を求めるのが面倒な場合から、そもそもトゥルーソリューションを求めるのが難しい場合も多い。こんな時は、$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}}$を使わない理由は全くない。

実装

以下は、リチャードソン誤差推定を通じてステップサイズを変えながら自動的に問題を解く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)