常微分方程式（ODE）の数値解法として、ルンゲ=クッタ法（Runge-Kutta method）は最も広く使われている手法の一つです。特に4次のルンゲ=クッタ法（RK4）は、実装の簡単さと精度のバランスに優れており、科学技術計算の定番です。
本記事では、オイラー法の復習から始めて、ルンゲ=クッタ法がどのように構成されるかを理解し、Pythonで実装して精度を比較します。
本記事の内容
- 常微分方程式の初期値問題
- オイラー法の復習と限界
- ルンゲ=クッタ法の導出
- 4次ルンゲ=クッタ法（RK4）の詳細
- Pythonでの実装と精度比較
初期値問題
次の形の常微分方程式の初期値問題を考えます。
$$ \frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0 $$
この方程式を解析的に解くことが困難な場合、数値的に近似解を求めます。時間を離散化し、$t_n = t_0 + nh$（$h$ は時間刻み幅）として、$y(t_n)$ の近似値 $y_n$ を逐次的に計算します。
オイラー法の復習
最も単純な数値解法がオイラー法です。
$$ y_{n+1} = y_n + h f(t_n, y_n) $$
これは $y(t)$ のテイラー展開の1次項までを用いた近似です。
$$ y(t + h) = y(t) + h y'(t) + \frac{h^2}{2} y”(t) + \cdots \approx y(t) + h f(t, y) $$
オイラー法の局所打ち切り誤差は $O(h^2)$、全体の誤差（大域誤差）は $O(h)$ です。つまり1次精度であり、高い精度を得るには非常に小さい $h$ が必要です。
ルンゲ=クッタ法のアイデア
ルンゲ=クッタ法は、1ステップの中で $f$ を複数回評価し、それらの加重平均を取ることで精度を向上させるアイデアに基づいています。
大雑把に言うと、「区間 $[t_n, t_{n+1}]$ の複数の点で傾きを評価し、それらを適切にブレンドすることで、より正確な傾きを推定する」手法です。
2次のルンゲ=クッタ法（中点法）
まず2次の方法を導出して考え方を理解しましょう。
$$ \begin{align} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + h, y_n + h k_1) \\ y_{n+1} &= y_n + \frac{h}{2}(k_1 + k_2) \end{align} $$
$k_1$ は区間の始点での傾き、$k_2$ はオイラー法で区間の終点に進んだ先での傾きです。両者の平均を取ることで、区間全体の傾きをより正確に近似しています。局所打ち切り誤差は $O(h^3)$、大域誤差は $O(h^2)$ です。
4次ルンゲ=クッタ法（RK4）
最も広く使われる4次のルンゲ=クッタ法は、以下の公式で定義されます。
$$ \begin{align} k_1 &= f(t_n, y_n) \\ k_2 &= f\left(t_n + \frac{h}{2}, \, y_n + \frac{h}{2} k_1\right) \\ k_3 &= f\left(t_n + \frac{h}{2}, \, y_n + \frac{h}{2} k_2\right) \\ k_4 &= f(t_n + h, \, y_n + h k_3) \\ y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{align} $$
各 $k_i$ の意味は次の通りです。
- $k_1$: 区間の始点での傾き
- $k_2$: $k_1$ を使って中点に進んだ先での傾き
- $k_3$: $k_2$ を使って中点に進んだ先での傾き（$k_2$ の修正版）
- $k_4$: $k_3$ を使って終点に進んだ先での傾き
最終的な更新は、これらの傾きの加重平均です。中点の傾き $k_2, k_3$ に2倍の重みがかかっているのは、シンプソンの公式と同じ思想です。
RK4の局所打ち切り誤差は $O(h^5)$、大域誤差は $O(h^4)$ です。つまり刻み幅 $h$ を半分にすると、誤差は $1/16$ に減少します。
Pythonでの実装
オイラー法、2次ルンゲ=クッタ法、4次ルンゲ=クッタ法を実装し、精度を比較します。
import numpy as np
import matplotlib.pyplot as plt
def euler_method(f, t_span, y0, h):
"""オイラー法"""
t = np.arange(t_span[0], t_span[1] + h, h)
y = np.zeros(len(t))
y[0] = y0
for i in range(len(t) - 1):
y[i+1] = y[i] + h * f(t[i], y[i])
return t, y
def rk2_method(f, t_span, y0, h):
"""2次ルンゲ=クッタ法（改良オイラー法）"""
t = np.arange(t_span[0], t_span[1] + h, h)
y = np.zeros(len(t))
y[0] = y0
for i in range(len(t) - 1):
k1 = f(t[i], y[i])
k2 = f(t[i] + h, y[i] + h * k1)
y[i+1] = y[i] + h / 2 * (k1 + k2)
return t, y
def rk4_method(f, t_span, y0, h):
"""4次ルンゲ=クッタ法"""
t = np.arange(t_span[0], t_span[1] + h, h)
y = np.zeros(len(t))
y[0] = y0
for i in range(len(t) - 1):
k1 = f(t[i], y[i])
k2 = f(t[i] + h/2, y[i] + h/2 * k1)
k3 = f(t[i] + h/2, y[i] + h/2 * k2)
k4 = f(t[i] + h, y[i] + h * k3)
y[i+1] = y[i] + h / 6 * (k1 + 2*k2 + 2*k3 + k4)
return t, y
# テスト問題: dy/dt = -y, y(0) = 1（解析解: y = e^(-t)）
f = lambda t, y: -y
y_exact = lambda t: np.exp(-t)
t_span = (0, 5)
y0 = 1.0
h = 0.5
# 各手法で計算
t_euler, y_euler = euler_method(f, t_span, y0, h)
t_rk2, y_rk2 = rk2_method(f, t_span, y0, h)
t_rk4, y_rk4 = rk4_method(f, t_span, y0, h)
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: 解の比較
t_fine = np.linspace(0, 5, 500)
axes[0].plot(t_fine, y_exact(t_fine), 'k-', linewidth=2, label='Exact: $e^{-t}$')
axes[0].plot(t_euler, y_euler, 'ro--', markersize=6, label=f'Euler (h={h})')
axes[0].plot(t_rk2, y_rk2, 'bs--', markersize=6, label=f'RK2 (h={h})')
axes[0].plot(t_rk4, y_rk4, 'g^--', markersize=6, label=f'RK4 (h={h})')
axes[0].set_xlabel('t')
axes[0].set_ylabel('y')
axes[0].set_title("Comparison of Numerical Methods ($dy/dt = -y$)")
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 右: 刻み幅と誤差の関係（収束次数の確認）
h_values = [1.0, 0.5, 0.25, 0.125, 0.0625]
errors_euler = []
errors_rk2 = []
errors_rk4 = []
for h_val in h_values:
t_e, y_e = euler_method(f, t_span, y0, h_val)
t_r2, y_r2 = rk2_method(f, t_span, y0, h_val)
t_r4, y_r4 = rk4_method(f, t_span, y0, h_val)
errors_euler.append(abs(y_e[-1] - y_exact(t_span[1])))
errors_rk2.append(abs(y_r2[-1] - y_exact(t_span[1])))
errors_rk4.append(abs(y_r4[-1] - y_exact(t_span[1])))
axes[1].loglog(h_values, errors_euler, 'ro-', label='Euler (Order 1)')
axes[1].loglog(h_values, errors_rk2, 'bs-', label='RK2 (Order 2)')
axes[1].loglog(h_values, errors_rk4, 'g^-', label='RK4 (Order 4)')
# 参考線
h_ref = np.array(h_values)
axes[1].loglog(h_ref, 0.5 * h_ref, 'r--', alpha=0.3, label='$O(h)$')
axes[1].loglog(h_ref, 0.3 * h_ref**2, 'b--', alpha=0.3, label='$O(h^2)$')
axes[1].loglog(h_ref, 0.1 * h_ref**4, 'g--', alpha=0.3, label='$O(h^4)$')
axes[1].set_xlabel('Step size h')
axes[1].set_ylabel('Global error at t=5')
axes[1].set_title('Convergence Order')
axes[1].legend(fontsize=8)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
右のグラフの対数プロットにおいて、各手法の誤差が理論通りの傾きで減少していることが確認できます。RK4は刻み幅 $h$ を半分にするだけで誤差が $1/16$ に減少するため、少ない計算量で高精度な結果が得られます。
まとめ
本記事では、ルンゲ=クッタ法について解説しました。
- ルンゲ=クッタ法は1ステップ内で複数回 $f$ を評価し、加重平均を取ることで高精度を実現する
- 4次ルンゲ=クッタ法（RK4）は大域誤差 $O(h^4)$ であり、実用上最も広く使われている
- オイラー法と比較して、同じ刻み幅でも格段に高い精度が得られる
- RK4の $k_1, k_2, k_3, k_4$ はそれぞれ区間の始点・中点・中点・終点での傾きの評価に対応する