Python
数値計算
物理
科学技術計算
計算物理学

[Pythonによる科学・技術計算] 常微分方程式,予測子-修正子法,アダムス-バッシュフォース-モルトン法,陽・陰解法

More than 1 year has passed since last update.

はじめに

1階常微分方程式の初期値問題,
$$ \frac{dy}{dx} = f(x,y), y(0)=a {\tag {1}}$$

を予測子-修正子法(補遺)で解くPythonコードを作成します。

予測子を得るために
$$y_{n+1} = y_n +\frac{h}{24}(55 f_n-59 f_{n-1}+37 f_{n-2}-9f_{n-3})+O(h^5){\tag {2}}$$

のアダムス・バッシュフォース(Adams-Bashforth)の4段公式を,

修正子は
$$y_{n+1} = y_n +\frac{h}{24}(9 f_{n+1}+ 19 f_{n}-5 f_{n-1}+9f_{n-2})+O(h^5){\tag {3}}$$

のアダムス・モルトン(Adams-Moulton)の4段公式を使います。

また,初期値以外に必要な$y$の値を4次のルンゲ・クッタ法で求めます。

この予測子・修正子法は,単純な陽解法であるルンゲ・クッタ法よりも計算速度が早く,
また大変安定な方法とされています。

どんな人向け?

● 予測子・修正子法を勉強したい人
● ルンゲ・クッタ法,オイラー法などの陽解法以外の勉強をしたい人


内容

$$ \frac{dy}{dx} = -y, \
y(0)=1 {\tag {4}}$$

を予測子-修正子法でx=0からx=20まで解きます。
刻み幅はdelta_x = 0.2とします。

次いで数値解を厳密解と比較します。
厳密解である特殊解は$y(x)=e^x$です。

また,他の数値解法であるオイラー法,アダムス・バッシュフォース法,4次のルンゲクッタ法,そしてさらに高段の予測子・修正子法と比較します。


計算コード

"""
予測子-修正子法
<アダムス・バッシュフォース4段法 (+4次のルンゲクッタ法)+ アダムス・モルトン4段法 O(h^5)
Oct. 2017
"""

import numpy as np
import matplotlib.pyplot as plt


# dy/dxを定義
def dydx (x,y):
    return y


#
x_min = 0
x_max = 20

delta_x = 0.2


y0=1



## y1,y2,y3を4次のルンゲクッタ法で求める
x=x_min
y=y0

x_lis=[x]
y_lis=[y0]
error_lis_pc=[0.0]
ii = 0
for ii in range(3):

    k1 = delta_x*dydx(x,y)
    k2 = delta_x*dydx(x+0.5*delta_x, y+0.5*k1)
    k3 = delta_x*dydx (x+0.5*delta_x, y+0.5*k2)
    k4 = delta_x*dydx(x+delta_x, y+k3)
    y += (k1+2*k2+2*k3+k4)/6     # O(h^5)
    x += delta_x
    x_lis.append(x)
    y_lis.append(y)
    error_lis_pc.append((y-np.exp(x))/np.exp(x))
#


n=2
eps = 1e-8
Iter_max = 50
while x <= x_max : 
    n = n+1



    ## 予測子-修正子法 ##    
    # アダムス-バッシュフォース4段法
    y_pred = y+delta_x*(55*dydx(x_lis[n],y_lis[n])-59*dydx(x_lis[n-1],y_lis[n-1])+37*dydx(x_lis[n-2],y_lis[n-2])-9*dydx(x_lis[n-3],y_lis[n-3]))/24


    #アダムス-モルトン4段法
    y_corr = y+(delta_x)*(9*dydx(x_lis[n]+delta_x,y_pred)+19*dydx(x_lis[n],y_lis[n])\
                          -5*dydx(x_lis[n-1],y_lis[n-1])+dydx(x_lis[n-2],y_lis[n-2]))/24
    delta_y = np.abs(y_pred-y_corr)

    iter=0
    yy = y_corr
    while  (delta_y > eps) and (iter <= Iter_max):
        iter +=1
        y_corr2 = y+(delta_x)*(9*dydx(x_lis[n]+delta_x,yy)+19*dydx(x_lis[n],y_lis[n])-5*dydx(x_lis[n-1],y_lis[n-1])+dydx(x_lis[n-2],y_lis[n-2]))/24
        #delta_y = np.abs((y_corr2-yy)/yy)
        delta_y = np.abs(y_corr2-yy)
        yy = y_corr2


    y = y_corr2

    x += delta_x

    x_lis.append(x)
    y_lis.append(y)
    error_lis_pc.append((y-np.exp(x))/np.exp(x))



#plot
#比較のための厳密解exp(x)をリストexact_solへセット
exact_sol=[]
xx=np.linspace(x_min,x_max,101)
for i in xx:
    exact_sol.append(np.exp(float(i)))
#


plt.plot(x_lis, y_lis, linewidth=1,label='PC') # 数値解
plt.plot(xx, exact_sol, color='red',linewidth=1, label='exact') #真の解

plt.xlabel('x', fontsize=24)
plt.ylabel('y', fontsize=24,rotation='horizontal')
plt.legend(loc='upper right')
plt.show()


plt.close()

# show error
plt.plot(x_lis, error_lis_pc, linewidth=1,label='error') # 数値解

plt.xlabel('x', fontsize=24)
plt.ylabel('y', fontsize=24,rotation='horizontal')
plt.legend(loc='upper right')
plt.show()



結果

スクリーンショット 2017-10-29 9.36.44.png

スクリーンショット 2017-10-29 9.36.49.png


他の手法との比較

(1)単純オイラー法

import numpy as np
import matplotlib.pyplot as plt

"""常備微分方程式: オイラー法"""
def dydx (x,y):
    return  y


#
x_min = 0
x_max = 20

delta_x = 0.2

x=x_min

y0=1
y=y0

x_lis=[]
y_lis=[]
error_lis=[]

x_lis.append(x_min)
y_lis.append(y0)
error_lis.append(0.0)

while x <= x_max : 
    y += dydx(x,y)*delta_x
    x += delta_x

    x_lis.append(x)
    y_lis.append(y)
    error_lis.append((y-np.exp(x))/np.exp(x))
#    print (x,y)

#plot
#比較のための厳密解(exp(-x**2/2))をリストexact_solへセット
exact_sol=[]
xx=np.linspace(x_min,x_max,101)
for i in xx:
    exact_sol.append(np.exp(i))
#


plt.plot(x_lis, y_lis, linewidth=1,label='Euler') # 数値解
plt.plot(xx, exact_sol, color='red',linewidth=1, label='exact') #真の解

plt.xlabel('x', fontsize=24)
plt.ylabel('y', fontsize=24,rotation='horizontal')
plt.legend(loc='upper right')
plt.show()


plt.close()

# show error
plt.plot(x_lis, error_lis, linewidth=1,label='error') # 数値解


plt.xlabel('x', fontsize=18)
plt.ylabel('Relative error', fontsize=18)
plt.legend(loc='upper right')
plt.show()


スクリーンショット 2017-10-29 11.31.19.png


(2) アダムス-バッシュフォース4段 陽解法



import numpy as np
import matplotlib.pyplot as plt


def dydx (x,y):
    return y


#
x_min = 0
x_max = 20

delta_x = 0.2


y0=1

## y1,y2,y3を4次のルンゲクッタ法で求める
x=x_min
y=y0

x_lis=[x]
y_lis=[y0]
error_lis_AB4=[0.0]

ii = 0
for ii in range(3):

    k1 = delta_x*dydx(x,y)
    k2 = delta_x*dydx(x+0.5*delta_x, y+0.5*k1)
    k3 = delta_x*dydx (x+0.5*delta_x, y+0.5*k2)
    k4 = delta_x*dydx(x+delta_x, y+k3)
    y += (k1+2*k2+2*k3+k4)/6     # O(h^5)
    x += delta_x
    x_lis.append(x)
    y_lis.append(y)
    error_lis_AB4.append((y-np.exp(x))/np.exp(x))
#


x=x_min+3*delta_x
y=y3



n=2
while x <= x_max : 
    n = n+1
    y += delta_x*(55*dydx(x_lis[n],y_lis[n])-59*dydx(x_lis[n-1],y_lis[n-1])+37*dydx(x_lis[n-2],y_lis[n-2])-9*dydx(x_lis[n-3],y_lis[n-3]))/24
    x += delta_x

    x_lis.append(x)
    y_lis.append(y)
    error_lis_AB4.append((y-np.exp(x))/np.exp(x))


#plot
#比較のための厳密解をリストexact_solへセット
exact_sol=[]
xx=np.linspace(x_min,x_max,101)
for i in xx:
    exact_sol.append(np.exp(float(i)))
#


plt.plot(x_lis, y_lis, linewidth=1,label='Adams-Bashforth_4') # 数値解
plt.plot(xx, exact_sol, color='red',linewidth=1, label='exact') #真の解

plt.xlabel('x', fontsize=24)
plt.ylabel('y', fontsize=24,rotation='horizontal')
plt.legend(loc='upper right')
plt.show()


plt.close()

# show error
plt.plot(x_lis, error_lis_AB4, linewidth=1,label='error') # 数値解


plt.xlabel('x', fontsize=18)
plt.ylabel('Relative error', fontsize=18)
plt.legend(loc='upper right')
plt.show()

スクリーンショット 2017-10-29 11.32.24.png


(3) 古典的ルンゲクッタ法 (4次のルンゲクッタ法)

import numpy as np
import matplotlib.pyplot as plt
"""
4次のルンゲ-クッタ法による1階常微分方程式
The fourth-order Runge-Kutta method for the 1st order ordinary differencial equation

dy/dx = -y を解く
"""
def dydx (x,y):
    return y


#
x_min = 0
x_max = 20

delta_x = 0.2

x=x_min

y0=1
y=y0

x_lis=[]
y_lis=[]
error_lis3=[]

x_lis.append(x_min)
y_lis.append(y0)
error_lis3.append(0.0)


while x <= x_max : 
    x += delta_x
    k1 = delta_x*dydx(x,y)
    k2 = delta_x*dydx(x+0.5*delta_x, y+0.5*k1)
    k3 = delta_x*dydx (x+0.5*delta_x, y+0.5*k2)
    k4 = delta_x*dydx(x+delta_x, y+k3)
    y += (k1+2*k2+2*k3+k4)/6     # O(h^5)
    x_lis.append(x)
    y_lis.append(y)
    error_lis3.append((y-np.exp(x))/np.exp(x))

#plot
#比較のための厳密解をリストexact_solへセット
exact_sol=[]
xx=np.linspace(x_min,x_max,101)
for i in xx:
    exact_sol.append(np.exp(float(i)))
#


plt.plot(x_lis, y_lis, linewidth=1,label='RK4') # 数値解
plt.plot(xx, exact_sol, color='red',linewidth=1, label='exact') #真の解

plt.xlabel('x', fontsize=24)
plt.ylabel('y', fontsize=24,rotation='horizontal')
plt.legend(loc='upper right')
plt.show()


plt.close()

# show error
plt.plot(x_lis, error_lis3, linewidth=1,label='error') # 数値解

plt.xlabel('x', fontsize=24)
plt.ylabel('y', fontsize=24,rotation='horizontal')
plt.legend(loc='upper right')
plt.show()

スクリーンショット 2017-10-29 11.34.31.png


(4) 予測子-修正子法: アダムス・バッシュフォース4段法 (+4次のルンゲクッタ法)+ アダムス・モルトン5段法

修正子のアダムス・モルトンの段数を1つ上げたもの。 局所打ち切り誤差はO(h^6)。

import numpy as np
import matplotlib.pyplot as plt


def dydx (x,y):
    return y


#
x_min = 0
x_max = 20

delta_x = 0.2


y0=1



## y1,y2,y3を4次のルンゲクッタ法で求める
x=x_min
y=y0

x_lis=[x]
y_lis=[y0]
error_lis_pc=[0.0]
ii = 0
for ii in range(3):

    k1 = delta_x*dydx(x,y)
    k2 = delta_x*dydx(x+0.5*delta_x, y+0.5*k1)
    k3 = delta_x*dydx (x+0.5*delta_x, y+0.5*k2)
    k4 = delta_x*dydx(x+delta_x, y+k3)
    y += (k1+2*k2+2*k3+k4)/6     # O(h^5)
    x += delta_x
    x_lis.append(x)
    y_lis.append(y)
    error_lis_pc.append((y-np.exp(x))/np.exp(x))
#




n=2
eps = 1e-12
Iter_max = 20
while x <= x_max : 
    n = n+1



    ## 予測子-修正子法 ##    
    # アダムス-バッシュフォース4段法
    y_pred = y+delta_x*(55*dydx(x_lis[n],y_lis[n])-59*dydx(x_lis[n-1],y_lis[n-1])+37*dydx(x_lis[n-2],y_lis[n-2])-9*dydx(x_lis[n-3],y_lis[n-3]))/24

    #アダムス-モルトン5段法
    y_corr = y+(delta_x)*(251*dydx(x_lis[n]+delta_x,y_pred)+646*dydx(x_lis[n],y_lis[n])\
                          -264*dydx(x_lis[n-1],y_lis[n-1])+106*dydx(x_lis[n-2],y_lis[n-2])-19*dydx(x_lis[n-3],y_lis[n-3]))/720
    delta_y = np.abs((y_pred-y_corr)/y_pred)

    iter=0
    yy = y_corr
    while  (delta_y > eps) and (iter <= Iter_max):
        iter +=1
        y_corr2 =  y+(delta_x)*(251*dydx(x_lis[n]+delta_x,yy)+646*dydx(x_lis[n],y_lis[n])\
                          -264*dydx(x_lis[n-1],y_lis[n-1])+106*dydx(x_lis[n-2],y_lis[n-2])-19*dydx(x_lis[n-3],y_lis[n-3]))/720
        delta_y = np.abs((y_corr2-yy)/yy)
        yy = y_corr2


    y = y_corr2

    x += delta_x

    x_lis.append(x)
    y_lis.append(y)
    error_lis_pc.append((y-np.exp(x))/np.exp(x))



#plot
#比較のための厳密解をリストexact_solへセット
exact_sol=[]
xx=np.linspace(x_min,x_max,101)
for i in xx:
    exact_sol.append(np.exp(float(i)))
#


plt.plot(x_lis, y_lis, linewidth=1,label='PC') # 数値解
plt.plot(xx, exact_sol, color='red',linewidth=1, label='exact') #真の解

plt.xlabel('x', fontsize=24)
plt.ylabel('y', fontsize=24,rotation='horizontal')
plt.legend(loc='upper right')
plt.show()


plt.close()

# show error
plt.plot(x_lis, error_lis_pc, linewidth=1,label='error') # 数値解


plt.xlabel('x', fontsize=18)
plt.ylabel('Relative error', fontsize=18)
plt.legend(loc='upper right')
plt.show()

スクリーンショット 2017-10-29 11.35.37.png


誤差の比較

オイラー法(Euler),アダムス・バッシュフォース4段(AB4),4次のルンゲクッタ法(RK4),予測子・修正子4段5位(PC),予測子・修正子5段6位(PC2)の相対誤差をプロットする。

PCはルンゲ・クッタ型と同程度の精度を持っている。PC2はさらに高精度である。

スクリーンショット 2017-10-29 10.02.32.png
図1

スクリーンショット 2017-10-29 10.02.40.png
図2. 図1からオイラー法の結果を削除したもの。

スクリーンショット 2017-10-29 10.02.56.png
図3. 図2からアダムス・バッシュフォース法の結果を削除したもの。


参考文献

[1] ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ, 丹慶 勝市ほか(訳), 技術評論社, 1993.

[2] Qiita記事, [Pythonによる科学・技術計算]4次のルンゲ-クッタ法による1次元ニュートン方程式の解法


補遺

以下,式(1)を
$$y_{n+1} = y_n+\int_{x_n}^{x_{n+1}}{f(x,y(x))dx} {\tag {5}}$$

と書き直したものを考えます。

(1)アダムス・バッシュフォース法

式(5)の被積分関数を$x=x_{n-1}$までの導関数値$f_{n-1}$までの情報から多項式を作り,$[x_n,x_{n+1}]$における関数を外挿によって求めます。たとえば一次式の場合,$[x_n,x_{n+1}]$では

$$f = \frac{(f_{n+1}-f_{n})}{h}(x-x_{n})+f_n {\tag {6}}$$

となります。これを式(5)に代入するとアダムス・バッシュフォースの2段公式

$$y_{n+1} = y_n +h(\frac{3}{2}f_n -\frac{1}{2}f_{n-1})+O(h^3){\tag {7}}$$
が得られます。

多項式を3次関数にしたものが式(2)の
$$y_{n+1} = y_n +\frac{h}{24}(55 f_n-59 f_{n-1}+37 f_{n-2}-9f_{n-3})+O(h^5)$$
でした。

後に述べるアダムス・モルトン法とは異なり,$y_{n+1}$を求めるのに必要な情報は$y_{n}$や$y_{n-1}$などでこれらはすでに計算できているものとされます。このような解法は陽解法と呼ばれます。

(2) アダムス・モルトン法

アダムス・バッシュフォース法は(5)の積分において被積分関数を領域外から外挿して求めていました。外挿は不安定性が伴います。そこで内挿を考えたものがアダムス・モルトン法です。線形内挿の場合は,

$$y_{n+1} = y_n +h(\frac{1}{2}f_{n+1} -\frac{1}{2}f_{n})+O(h^3){\tag {8}}$$

となることが容易に示されます。アダムス・バッシュフォース法と異なり,右辺の$f_{n+1} = f(x_
{n+1},y_{n+1}$に$y_{n+1}$があります。左辺の$y_{n+1}$を求めるためには$y_{n+1}$の情報が必要です。この場合,式(8)は漸化式としてみるのではなくむしろ$y_{n+1}$についての非線形方程式とみなされるべきです。

$$y_{n+1} - y_n -h(\frac{1}{2}f_{n+1} -\frac{1}{2}f_{n})+O(h^3)=0 {\tag {8'}}$$

このような解法を陰解法と呼びます。

内挿関数を3次関数にした場合が式(3)の

$$y_{n+1} = y_n +\frac{h}{24}(9 f_{n+1}+ 19 f_{n}-5 f_{n-1}+9f_{n-2})+O(h^5)$$

でした。

(3) 予測子・修正子法

上記の(1)と(2)の良いところが合わさっています。
式(8)の$f_{n+1}$の計算に必要な$y_{n+1}$をアダムス・バッシュフォースの式(7)の解を用い,それによって$y_{n+1}$を修正します。得られた$y_{n+1}$(修正解)を再び(8)の右辺の$f_{n+1}$の計算の予測子として用い,予測子解と修正子解との差が数値誤差(公式の位数程度)まで小さくなるまで反復計算します。

最良か?

文献[1]では,
常微分方程式の初期値問題を解くための実用的な方法として

  1. ルンゲ・クッタ法
  2. リチャードソン補外法の派生型であるBulirsch-Stoer法
  3. 予測子・修正子法

の3つを挙げそれらの有効性を比較検討し,「予測子・修正子法はかなり洗練させたコードにしないとBulirsch-Stoer法に太刀打ちできないようだ」と述べています。ルンゲ・クッタ法とBulirsch-Stoer法の中間に位置するようです。