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

#はじめに
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次のルンゲ・クッタ法](https://qiita.com/sci_Haru/items/b8d5c9cfe64c4366630a)で求めます。

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

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

***

# 内容

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

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

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

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


***

#計算コード


```py3

"""
予測子-修正子法
<アダムス・バッシュフォース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](https://qiita-image-store.s3.amazonaws.com/0/192457/97728d50-cc5b-32ba-9039-ee94b6c04f4f.png "スクリーンショット 2017-10-29 9.36.44.png")




![スクリーンショット 2017-10-29 9.36.49.png](https://qiita-image-store.s3.amazonaws.com/0/192457/1f2b7b62-da33-4cc2-2ec9-fc1388278a43.png "スクリーンショット 2017-10-29 9.36.49.png")


***

#他の手法との比較

##(1)単純オイラー法

```py3

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](https://qiita-image-store.s3.amazonaws.com/0/192457/d6deffe6-40de-5253-4717-593968c480c9.png)

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

```py3



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](https://qiita-image-store.s3.amazonaws.com/0/192457/a012612f-e77b-9bc3-7f65-387ee671f59b.png)

***

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

```py3

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](https://qiita-image-store.s3.amazonaws.com/0/192457/e83b4fc9-a497-6d8a-ecf3-cbf8166191b2.png)

***

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

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


```py3
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](https://qiita-image-store.s3.amazonaws.com/0/192457/b2063fc6-64c9-9150-0ebe-b1f4d9340508.png)





***
#誤差の比較

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

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


![スクリーンショット 2017-10-29 10.02.32.png](https://qiita-image-store.s3.amazonaws.com/0/192457/d000b6bf-8a17-7278-476e-19b86efc5770.png "スクリーンショット 2017-10-29 10.02.32.png")
図1



![スクリーンショット 2017-10-29 10.02.40.png](https://qiita-image-store.s3.amazonaws.com/0/192457/761cdf16-e761-d4d6-4fec-7c221de91265.png "スクリーンショット 2017-10-29 10.02.40.png")
図2. 図1からオイラー法の結果を削除したもの。





![スクリーンショット 2017-10-29 10.02.56.png](https://qiita-image-store.s3.amazonaws.com/0/192457/4f97fd20-e05a-042d-a5df-7c20e37e621c.png "スクリーンショット 2017-10-29 10.02.56.png")
図3. 図2からアダムス・バッシュフォース法の結果を削除したもの。

***
#参考文献
[1] [ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ](https://www.amazon.co.jp/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%A1%E3%83%AA%E3%82%AB%E3%83%AB%E3%83%AC%E3%82%B7%E3%83%94%E3%83%BB%E3%82%A4%E3%83%B3%E3%83%BB%E3%82%B7%E3%83%BC-%E6%97%A5%E6%9C%AC%E8%AA%9E%E7%89%88%E2%80%95C%E8%A8%80%E8%AA%9E%E3%81%AB%E3%82%88%E3%82%8B%E6%95%B0%E5%80%A4%E8%A8%88%E7%AE%97%E3%81%AE%E3%83%AC%E3%82%B7%E3%83%94-William-H-Press/dp/4874085601),  丹慶 勝市ほか(訳), 技術評論社, 1993.


[2] Qiita記事， [[Pythonによる科学・技術計算]4次のルンゲ-クッタ法による1次元ニュートン方程式の解法](https://qiita.com/sci_Haru/items/b8d5c9cfe64c4366630a)
***


#補遺
以下，式(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法の中間に位置するようです。

***