数値計算② ルンゲ-クッタ法【Pythonと物理#6】

※Amazonのアソシエイトとして、当メディア(Nラボ備忘録)は適格販売により収入を得ています。

今回のターゲット：ルンゲ-クッタ法（Runge–Kutta Method）
ルンゲ-クッタ法（Runge–Kutta Method）
ルンゲ-クッタ法で数値計算をしてみる
Pythonで計算する
実行結果
最後に

今回のターゲット：ルンゲ-クッタ法（Runge–Kutta Method）

　前回はオイラー法を勉強してみました。

数値計算① オイラー法【Pythonと物理#5】

この記事は、オイラー法について紹介しています。初学者向けの内容となっています。

今回はオイラー法よりも精度の高い ルンゲ-クッタ法 を勉強してみました。

勉強の際に参考にさせていただいた書籍もご紹介させていただきます。具体的なPythonのコードも載っており、大変勉強になる１冊です。

リンク

ルンゲ-クッタ法（Runge–Kutta Method）

ルンゲ-クッタ法は以下の式で定義されます。

$$
y_{n+1} = y_{n} + \cfrac{h}{6}\left(k_{1} + 2k_{2} + 2k_{3} + k_{4}\right)
$$

ただし、$t_{n+1} = t_{n} + h$ です。
また、$k_{1}$、$k_{2}$、$k_{3}$、$k_{4}$は以下のように表されます。

$$
\begin{eqnarray}
&k_{1}& = f(t_{n}, y_{n})\\
&k_{2}& = f\left(t_{n} + \cfrac{h}{2}, y_{n} + \cfrac{h}{2}k_{1}\right)\\
&k_{3}& = f\left(t_{n} + \cfrac{h}{2}, y_{n} + \cfrac{h}{2}k_{2}\right)\\
&k_{4}& = f\left(t_{n} + h, y_{n} + hk_{3}\right)
\end{eqnarray}
$$

中々難解ですね。あくまでもイメージですが、$k_{1}$、$k_{2}$、$k_{3}$、$k_{4}$を以下のように理解しました。

　$k_{1}$は、$t_{n}$ における傾き（オイラー法で求まるものと同じ）
　$k_{2}$は、$t_{n}$ から傾き $k_{1}$ で $h/2$ 進んだ時の、$t_{n} + h/2$ における傾き
　$k_{3}$は、$t_{n}$ から傾き $k_{2}$ で $h/2$ 進んだ時の、$t_{n} + h/2$ における傾き
　$k_{4}$は、$t_{n} + k_{2}$ から傾き $k_{3}$ で $h/2$ 進んだ時の、$t_{n} + h$ における傾き

また、定義式をあらためて見てみます。

$$
y_{n+1} = y_{n} + \cfrac{h}{6}\left(k_{1} + 2k_{2} + 2k_{3} + k_{4}\right)
$$

よく見ると、$1/6$ やら、$k_{2}$ と $k_{3}$ の係数に $2$ がついています。
$1/6$は$k_{1}$、$k_{2}\times2$、$k_{3}\times2$、$k_{4}$、合計6つの数値の平均を取っています。
また、$k_{2}$ と $k_{3}$ の係数に $2$ をつけることで、重み付けしています。

今回勉強したルンゲ-クッタ法は、4次の公式です。そのため、精度は $O\left(h^{4}\right)$ となります。後述する計算結果を見ていただければ、オイラー法に比べてだいぶ精度はよくなっていることがわかると思います。

ルンゲ-クッタ法で数値計算をしてみる

それではルンゲクッタ法を使って数値計算をしてみます。今回は、下の式をルンゲ-クッタ法で数値計算してみます。

$$
\cfrac{dy}{dt} = y
$$

これは普通に厳密解 $y$ を求めることもできるので、求めてみると以下のようになります。

$$
y = y_{0}\exp\left(t\right)
$$

ただし、 $y_{0}$ は定数です。

Pythonで計算する

コードは下のようにしました。詳細はいつも通りコメントアウトで書いています。
今回は厳密解を求められる常微分方程式だったので、厳密解も同じグラフに描画するコードにしています。また、前回同じ式をオイラー法で計算しているので、そちらも同じグラフに描画するコードにしています。

import numpy as np
import matplotlib.pyplot as plt

#数値計算したい微分方程式の定義
def func(y):
    return y

#ルンゲ-クッタ法
def RungeKutta_method(f, y, h):
    k1 = f(y)
    k2 = f(y + k1*h/2)
    k3 = f(y + k2*h/2)
    k4 = f(y + k3*h/2)
    return y + (k1 + 2*k2 + 2*k3 + k4)*h/6

#オイラー法
def euler_method(f, y, h):
    return y + f(y) * h

#常微分方程式を計算
def ode_solve_RKM(f, y0, T, N):
    t = np.linspace(0, T, N)
    h = t[1] - t[0]
    y = [y0, ]
    for i in range(N - 1):
        y.append(RungeKutta_method(f, y[-1], h))
    return t, np.array(y)
def ode_solve_EM(f, y0, T, N):
    t = np.linspace(0, T, N)
    h = t[1] - t[0]
    y = [y0, ]
    for i in range(N - 1):
        y.append(euler_method(f, y[-1], h))
    return t, np.array(y)
    
#初期条件などを設定して、上で定義した関数たちを呼び出す
y0 = 0.01
T = 10
N =101
t, y_RKM = ode_solve_RKM(func, y0, T, N) #ルンゲ-クッタ法
t, y_EM = ode_solve_EM(func, y0, T, N) #オイラー法
y_exact = y0 * np.exp(t) #math.exp()だと配列ではないので計算不可


#以下は可視化のためにグラフを書くコード
fig, ax = plt.subplots()
ax.plot(t, y_RKM, color='orange', label='Runge-Kutta Method')
ax.plot(t, y_EM, color='crimson', label='Euler Method')
ax.plot(t, y_exact, color='cornflowerblue', label='Exact')
    
ax.set_xlabel('t')
ax.set_ylabel('y')
ax.legend()
    
ax.tick_params(direction='in', axis='x')
ax.tick_params(direction='in', axis='y')

#以下の範囲設定は適当です
ax.set_xlim(0, 5)
ax.set_ylim(0, 1)

plt.show()