LoginSignup
17
14

More than 5 years have passed since last update.

@sci_Haru(Haru_M D)

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

Last updated at Posted at 2017-10-29

はじめに

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法の中間に位置するようです。

17
14
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
17
14