はじめに
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()
結果
他の手法との比較
(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()
(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()
(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()
(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()
誤差の比較
オイラー法(Euler),アダムス・バッシュフォース4段(AB4),4次のルンゲクッタ法(RK4),予測子・修正子4段5位(PC),予測子・修正子5段6位(PC2)の相対誤差をプロットする。
PCはルンゲ・クッタ型と同程度の精度を持っている。PC2はさらに高精度である。
図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]では,
常微分方程式の初期値問題を解くための実用的な方法として
- ルンゲ・クッタ法
- リチャードソン補外法の派生型であるBulirsch-Stoer法
- 予測子・修正子法
の3つを挙げそれらの有効性を比較検討し,「予測子・修正子法はかなり洗練させたコードにしないとBulirsch-Stoer法に太刀打ちできないようだ」と述べています。ルンゲ・クッタ法とBulirsch-Stoer法の中間に位置するようです。