# ルンゲクッタ法による常微分方程式の数値解析

ここでは，オイラー法よりも計算精度が良く，広く用いられている4次のルンゲ・クッタ（Runge-Kutta）法を学ぶ．

オイラー法のときと同様に，2階常微分方程式を1階化した

$$
\begin{cases}
 \frac{dx(t)}{dt} = v(t), \\
 \frac{dv(t)}{dt} = f(x(t), v(t))
\end{cases}
\tag{2, 再掲}
$$

を考える．式（2）の数値解 $x_i$, $v_i$，$i=0, 1, 2, \ldots,$は，４次のルンゲ・クッタ法により次のように計算される．

$$
\begin{eqnarray}
x_{i+1} &= x_i + \frac{h}{6} \left(k_{x1} + 2k_{x2} + 2k_{x3} + k_{x4} \right),
\tag{9}
\\
%
v_{i+1} &= v_i + \frac{h}{6} \left(k_{v1} + 2k_{v2} + 2k_{v3} + k_{v4} \right)
\tag{10}
\end{eqnarray}
$$

ここで，

$$
\begin{eqnarray}
&k_{x1} = v_i, \quad
&k_{v1} = f(x_i, v_i),\\
%
&k_{x2} = v_i + \frac{h k_{v1}}{2}, \quad
&k_{v2} = f\left(x_i + \frac{h k_{x1}}{2}, v_i + \frac{h k_{v1}}{2} \right),\\
%
&k_{x3} = v_i + \frac{h k_{v2}}{2}, \quad
&k_{v3} = f\left(x_i + \frac{h k_{x2}}{2}, v_i + \frac{h k_{v2}}{2} \right),\\
%
&k_{x4} = v_i + h k_{v3}, \quad
&k_{v4} = f(x_i + h k_{x3}, v_i + h k_{v3}),
\end{eqnarray}
$$

であり，$h$は積分刻みである．

では，ルンゲクッタ法をPythonで実装してみよう．ここでは

$$
\begin{align}
  f = -x(t) \tag{5}\\
\end{align}
$$

とする．また，初期条件を$x(0)=1.0, y(0)=0$とする．

このとき，ルンゲクッタ法を実装すると下記のようになる．

<font color="red">
（TODO） 下記のコードを理解し，セルを実行せよ
</font>


```
# 必要なライブラリのインポート
import numpy as np
from matplotlib import pyplot as plt

# 関数定義
def f(x, v):
  return -x # dv/dt = f(x, y) = -x の実装

# 計算条件を変数として定義
x = 1.0 # xの初期値
v = 0 # yの初期値
t_ini = 0 # 初期時刻
h = 0.2 # 積分刻み
t_max = 10 + h # 終了時刻

# 結果保存用の配列定義
ts = np.arange(start=t_ini, stop=t_max, step=h) # tの時系列
xs = np.zeros(len(ts)) # xの時系列
vs = np.zeros(len(ts)) # yの時系列

# 初期値の代入
xs[0] = x
vs[0] = v

# ルンゲクッタ法
for i in range(1, len(ts)): # i=1 から i=len(ts) までループ
  kx1 = v
  kv1 = f(x, v)
  kx2 = v + h * kv1 / 2.0
  kv2 = f(x + h * kx1 / 2.0, v + h * kv1 / 2.0)
  kx3 = v + h * kv2 / 2.0
  kv3 = f(x + h * kx2 / 2.0, v + h * kv2 / 2.0)
  kx4 = v + h * kv3
  kv4 = f(x + h * kx3, v + h * kv3)
  # v, xの更新
  v += (kv1 + 2.0 * kv2 + 2.0 * kv3 + kv4) * h / 6.0
  x += (kx1 + 2.0 * kx2 + 2.0 * kx3 + kx4) * h / 6.0
  # 結果の保存
  xs[i] = x
  vs[i] = v
```

ルンゲ・クッタ法で計算された数値解が配列`xs`, `ys`に格納される．また，数値計算の各時刻は`ts`に格納される．

<font color="red">
（TODO） 下記のコードを実行し，`ts`を横軸，`xs`を縦軸にプロットせよ．
</font>


```
plt.plot(ts, xs, '.')
plt.xlabel('t')
plt.ylabel('x')
```

<font color="red">
（TODO） 下記のコードを実行し，`ts`を横軸，`vs`を縦軸にプロットせよ．
</font>


```
plt.plot(ts, vs, '.')
plt.xlabel('t')
plt.ylabel('v')
```

オイラー法のときと同様に，得られた数値解と厳密解を比較してみよう．

<font color="red">
（TODO） 下記のコードを実行し，得られた厳密解とオイラー法の数値解を同時にplotせよ．
</font>


```
xs_exact = np.cos(ts) # 厳密解

# 厳密解とオイラー法の数値解をplot
plt.plot(ts, xs_exact, '-r', label='Exact')
plt.plot(ts, xs, '.g', label='RK') # 既にルンゲ・クッタ法で計算済み
plt.xlabel('t')
plt.ylabel('x')
plt.legend() # labelを表示する
```

<font color="red">
（TODO） 次のコードを実行し，ルンゲ・クッタ法で得た数値解の配列`xs`と，厳密解の配列`xs_exact`のMAEを計算せよ．
</font>


```
np.mean(np.abs(xs - xs_exact))
```

出力された結果から，今回のルンゲ・クッタ法では，$4.0\times10^{-5}$程度のMAEとなるはずだ．よって，ルンゲ・クッタ法はオイラー法よりも高精度な数値解が得られることが確かめられた．

# ここまでのまとめ

数値解析手法であるオイラー法とルンゲ・クッタ法をPythonを用いて実装し，その動作を確認した．両手法ともに，常微分方程式の1階化ができれば，適用できることがわかっただろう．

ここまでの内容を一通り理解したら，[実験のTopページ](https://github.com/yyamnk/numerical-methods-py3/blob/master/uu_experiment.md) に戻り，レポート課題に取り組んでほしい．

fb8a58f2258d2d9d66c2fc8bdabba48aa73f231a
ipynb
Jupyter Notebook
exp_python3.ipynb
yyamnk/numerical-methods-py3
c2ba7d516b81348cb57744c78d9238f11e3104ce
MIT
exp_python3.ipynb
yyamnk/numerical-methods-py3
c2ba7d516b81348cb57744c78d9238f11e3104ce
MIT
exp_python3.ipynb
yyamnk/numerical-methods-py3
c2ba7d516b81348cb57744c78d9238f11e3104ce
MIT
2020-07-01T06:31:13.000Z
2020-07-01T06:31:13.000Z
Qwen/Qwen-72B
1. YES
2. YES