常微分方程式の数値解法 — オイラー法からルンゲ=クッタ法まで

人工衛星の軌道、化学反応の進行、ニューラルネットワークの学習 — これらはすべて微分方程式で記述されますが、そのほとんどは解析的に解けません。解析的に解ける微分方程式は、実はごく一部の特殊なケースに限られています。たとえば、地球と月と太陽の3体問題は解析解を持たず、NASAの月探査ミッションの軌道計算はすべて数値解法に頼っています。

しかし、微分の「定義」に立ち返ると、数値的に解く方法が見えてきます。

$$ y'(t) \approx \frac{y(t + h) – y(t)}{h} $$

この素朴な近似がオイラー法の出発点です。しかし、実用的な精度を得るにはより洗練された手法が必要です。精度だけでなく安定性も重要な課題です。一部の問題では、計算が発散してしまう手法があり、適切な手法選択が不可欠になります。本記事では、常微分方程式（ODE）の数値解法を精度と安定性の観点から体系的に解説します。

ODE数値解法を理解すると、以下のような場面で活用できます。

• 軌道力学: 人工衛星・惑星の軌道計算（N体問題）。はやぶさのイオンエンジンによる軌道は数値積分で計算されました
• 化学工学: 反応速度論に基づく化学反応シミュレーション。燃焼のような速い反応と遅い反応が混在するスティッフ系が典型的です
• 制御工学: 状態空間モデル $\dot{\bm{x}} = A\bm{x} + B\bm{u}$ のシミュレーション。リアルタイム制御では計算速度が厳しく制約されます
• 天気予報: 大気の運動方程式（ナビエ・ストークス方程式の離散化）の時間積分。数十億個のODEを効率的に解く必要があります

本記事の内容

• 前進オイラー法と局所打切り誤差
• 後退オイラー法と陰解法
• ルンゲ=クッタ法（2次・4次）の導出
• 安定性領域とスティッフ問題
• 適応的刻み幅制御
• Pythonでの実装と精度比較

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• テイラー展開の基礎 — 打切り誤差の解析
• 連立常微分方程式の理論 — 連立系のベクトル表記
• 補間と近似 — 多段法の基礎

初期値問題と数値解法の基本

問題設定

初期値問題（IVP: initial value problem）は

$$ y'(t) = f(t, y), \quad y(t_0) = y_0 $$

で定義されます。ここで $f$ は右辺関数と呼ばれ、微分方程式の「構造」を決めるものです。初期条件 $y(t_0) = y_0$ は出発点を指定します。リプシッツ条件（$|f(t, y_1) – f(t, y_2)| \leq L|y_1 – y_2|$）が満たされれば、解の存在と一意性が保証されます（ピカール=リンデレフの定理）。

数値解法の基本的な考え方は、連続的な時間軸を離散的な点に置き換えることです。時刻を $t_n = t_0 + nh$（$h$ は刻み幅）で離散化し、$y_n \approx y(t_n)$ を逐次的に計算します。出発点 $y_0$ から始めて、$y_1, y_2, \ldots$ を順に求めていくことで、解の軌跡を離散的に追跡します。

手法の分類

1段法は自己起動可能で扱いやすく、多段法は計算効率に優れますが起動に別の方法（通常RK法）が必要です。陽解法と陰解法の違いは、$y_{n+1}$ を求める際に $y_{n+1}$ 自身の情報を必要とするか否かにあります。陰解法は各ステップで方程式を解く必要があるため計算コストが高いですが、安定性に優れるという重要な利点があります。

まずは最も基本的で直感的なオイラー法から始めましょう。

オイラー法

前進オイラー法（陽的オイラー法）

オイラー法の発想は非常に素朴です。現在地点 $(t_n, y_n)$ における接線の傾き $f(t_n, y_n)$ がわかっているので、その方向に刻み幅 $h$ だけ「まっすぐ進む」のです。直線近似なので曲がった解軌道からは少しずつずれますが、$h$ を小さくすればそのずれを抑えられます。

数学的には、テイラー展開の1次近似

$$ y(t + h) = y(t) + hy'(t) + O(h^2) $$

から、$y’ = f(t, y)$ を代入して

$$ \begin{equation} y_{n+1} = y_n + hf(t_n, y_n) \end{equation} $$

幾何学的には、$(t_n, y_n)$ における接線の方向に $h$ だけ進むことに対応します。

局所打切り誤差

1ステップの打切り誤差（local truncation error, LTE）は

$$ \tau_{n+1} = y(t_{n+1}) – y(t_n) – hf(t_n, y(t_n)) = \frac{h^2}{2}y”(\xi_n) $$

LTE が $O(h^2)$ なので、前進オイラー法は1次の方法（order 1）と言います。ここで「LTEが $O(h^2)$」なのに「1次の方法」と呼ぶのは混乱しやすい点です。理由は、LTEは「1ステップの誤差」であり、「大域誤差」はLTEの蓄積で決まるからです。

大域誤差

区間 $[a, b]$ を $N = (b – a)/h$ ステップで進むと、各ステップの LTE $O(h^2)$ が $N = O(1/h)$ 回蓄積します。最悪の場合、誤差は $N \times O(h^2) = O(h)$ となります。

$$ |y(t_N) – y_N| \leq C h $$

$O(h)$ です。刻み幅を半分にすると精度は約2倍になります。一般に、LTEが $O(h^{p+1})$ の方法は $p$ 次の方法と呼ばれ、大域誤差は $O(h^p)$ になります。この「1次下がる」関係は、誤差の蓄積ステップ数が $O(1/h)$ であることに起因しています。

後退オイラー法（陰的オイラー法）

テイラー展開を $t_{n+1}$ から後ろ向きに行うと

$$ \begin{equation} y_{n+1} = y_n + hf(t_{n+1}, y_{n+1}) \end{equation} $$

右辺に $y_{n+1}$ が含まれているため、一般には非線形方程式を解く必要があります（ニュートン法など）。計算コストは高いですが、安定性に優れるためスティッフ問題に不可欠です。

前進オイラー法と後退オイラー法の違いを直感的に理解するには、「未来を予測する」のか「未来に適応する」のかという視点が有用です。前進オイラー法は「現在の傾きで未来を予測する」ので、傾きが急激に変わる場面では的外れになりがちです。後退オイラー法は「未来の傾きを使って現在から進む」ので、急変する場面でも適切に対応できます。ただし、「未来の傾き」を知るためには方程式を解く必要があるため、計算コストが増すのです。

オイラー法の精度は1次で、実用には不十分なことが多いです。$h$ を半分にしても精度が2倍にしかならないため、高精度を得るには膨大な計算が必要になります。テイラー展開をもう少し活用すれば、同じ関数評価回数でより高精度な方法が得られます。

ルンゲ=クッタ法

基本思想

オイラー法が「1点の傾き情報だけで次の値を決める」のに対し、ルンゲ=クッタ法は「区間内の複数の点で傾きを調べ、その情報を総合して次の値を決める」という戦略を取ります。偵察隊を何人か先に派遣して道の状況を調べ、その報告を総合して最適な進路を決めるようなイメージです。

ルンゲ=クッタ（RK）法の基本思想は、$[t_n, t_{n+1}]$ の区間内で $f$ を複数回評価し、その加重平均を傾きとして使うことです。

$$ y_{n+1} = y_n + h\sum_{i=1}^{s} b_i k_i $$

ここで $k_i = f(t_n + c_i h, y_n + h\sum_j a_{ij} k_j)$ は段（stage）と呼ばれます。パラメータ $a_{ij}, b_i, c_i$ を適切に選ぶことで、高い精度を達成します。

2次のルンゲ=クッタ法（修正オイラー法）

最も単純な改良として、区間の中点で傾きを評価する方法を考えましょう。オイラー法は $t_n$ での傾きだけを使いますが、この傾きは $t_n$ から $t_{n+1}$ に進む間に変化する可能性があります。ならば、区間の中点 $t_n + h/2$ での傾きを使えば、区間全体をより良く代表する傾きが得られるはずです。

$$ k_1 = f(t_n, y_n) $$

$$ k_2 = f\left(t_n + \frac{h}{2}, y_n + \frac{h}{2}k_1\right) $$

$$ \begin{equation} y_{n+1} = y_n + hk_2 \end{equation} $$

テイラー展開で確認すると、LTE は $O(h^3)$ で、2次の方法です。

直感的には、「まず半分進んでみて、その点での傾きを使って全ステップを進む」ということです。これにより、曲率の情報が反映されて精度が向上します。

ホイン法（改良オイラー法）

もう一つの2次RK法で、台形則に対応します。

$$ k_1 = f(t_n, y_n), \quad k_2 = f(t_n + h, y_n + hk_1) $$

$$ y_{n+1} = y_n + \frac{h}{2}(k_1 + k_2) $$

区間の両端での傾きの平均を使います。ホイン法は「予測子-修正子法」としても理解できます。まずオイラー法で $y_{n+1}$ を「予測」し（$\tilde{y}_{n+1} = y_n + hk_1$）、次にその予測値での傾き $k_2$ と元の傾き $k_1$ の平均で「修正」するのです。

修正オイラー法とホイン法はどちらも2次のRK法ですが、内部の評価点が異なります。修正オイラー法は中点法則、ホイン法は台形法則に対応しており、数値積分との深い関連が見えます。

古典的4次ルンゲ=クッタ法（RK4）

最も広く使われるODE数値解法であり、多くの教科書で「ルンゲ=クッタ法」と言えばこれを指します。4段4次の精度を持ちます。

$$ 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 + hk_3) $$

$$ \begin{equation} y_{n+1} = y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{equation} $$

重み $1:2:2:1$ はシンプソン則と同じパターンです。LTE は $O(h^5)$ で、大域誤差は $O(h^4)$ です。

RK4の直感

$k_1$ は区間の左端、$k_2, k_3$ は中点、$k_4$ は右端での傾きの近似値です。中点の傾きを2回評価して2倍の重みを付けているのは、シンプソン則が中点に大きな重みを付けるのと同じ理由です。

ブッチャー配列

RK法のパラメータを整理するブッチャー配列（Butcher tableau）は

$$ \begin{array}{c|c} \bm{c} & A \\ \hline & \bm{b}^T \end{array} $$

RK4の場合

$$ \begin{array}{c|cccc} 0 & & & & \\ 1/2 & 1/2 & & & \\ 1/2 & 0 & 1/2 & & \\ 1 & 0 & 0 & 1 & \\ \hline & 1/6 & 1/3 & 1/3 & 1/6 \end{array} $$

$A$ が狭義下三角行列（対角以下のみ非零）であれば陽的RK法で、$k_i$ を $k_1, k_2, \ldots$ の順に逐次計算できます。$A$ の対角成分や上三角部分にも非零要素がある場合は陰的RK法で、$k_1, \ldots, k_s$ を同時に（連立方程式を解いて）決定する必要があります。

RK法の次数と段数の関係

$s$ 段の陽的RK法で達成可能な最大次数は、$s$ が小さいうちは $p = s$ ですが、$s \geq 5$ では $p < s$ となります。具体的には、5次の方法には最低6段が必要です。この「次数の壁」は、テイラー展開の照合条件（order conditions）が段数とともに急増するためです。4段4次のRK4は、この意味で非常に効率的な方法と言えます。

ルンゲ=クッタ法の精度を理解したところで、次は安定性の問題を見てみましょう。高精度な方法でも、安定性を満たさなければ計算は破綻します。

安定性とスティッフ問題

なぜ安定性が重要か

精度と安定性は、ODE数値解法の「車の両輪」です。精度の高い方法でも、安定性条件を満たさなければ数値解は発散して意味を持ちません。安定性の問題は特にスティッフな方程式で深刻になります。

テスト方程式

ODE数値解法の安定性は、テスト方程式（ダールキストのテスト方程式とも呼ばれます） $y’ = \lambda y$（$\lambda \in \mathbb{C}$, $\text{Re}(\lambda) < 0$）で解析します。この方程式の厳密解は $y(t) = y_0 e^{\lambda t}$ で、$\text{Re}(\lambda) < 0$ なので $t \to \infty$ で $y \to 0$ に収束します。数値解法が「安定」であるとは、数値解もこの減衰を正しく再現することを意味します。

前進オイラー法を適用すると

$$ y_{n+1} = (1 + h\lambda)y_n = (1 + h\lambda)^n y_0 $$

$|1 + h\lambda| < 1$ でないと解が発散します。つまり、増幅因子 $R(z) = 1 + z$（$z = h\lambda$）の絶対値が1未満でなければならないのです。これが安定性条件です。安定性条件が破れると、真の解は減衰しているにもかかわらず、数値解が指数的に増大するという致命的な状況が生じます。

安定性領域

安定性領域（stability region）$\mathcal{S}$ は、$z = h\lambda$ の複素平面上で $|R(z)| < 1$（$R(z)$ は増幅因子）を満たす領域です。

スティッフ問題

スティッフ（stiff）な問題とは、系の中に非常に異なる時定数が混在する問題です。

たとえば $y’ = -1000y + 1000t + 1$（$\lambda = -1000$）では、過渡応答は $e^{-1000t}$ で急速に減衰しますが、定常解 $y = t$ はゆっくり変化します。

前進オイラー法で安定に計算するには $h < 2/1000 = 0.002$ が必要ですが、定常解のスケールからは $h \approx 0.1$ で十分な精度が得られます。つまり、安定性のために精度上は不要に小さい刻み幅を使わざるを得ません。この「安定性による刻み幅の制約が精度の要求よりもはるかに厳しい」状況こそが、スティッフ問題の本質的な困難さです。

A安定な方法（安定性領域が左半平面全体を含む）はこの問題を回避できます。A安定であれば、$\text{Re}(\lambda) < 0$ である限り、どんなに大きな $|\lambda|$ に対しても任意の刻み幅 $h$ で安定に計算できます。後退オイラー法と台形法はA安定で、スティッフ問題に適しています。RK4は安定性領域が有界であるため、A安定ではなくスティッフ問題には不向きです。

ただし、A安定な方法の中にも優劣があります。台形法はA安定ですが、$\text{Re}(\lambda)$ が非常に負の場合に数値解が振動することがあります（この問題を避けるL安定性という概念もあります）。後退オイラー法やRadau IIA法はL安定で、強いスティッフ性にも対応できます。

陰的RK法のA安定性

スティッフ問題に対しては、陰的RK法（例: Radau IIA法）やBDF法（後退差分法）が使われます。SciPyの solve_ivp の method='Radau' や method='BDF' はこれらに対応しています。

実用上のスティッフ問題の例としては、電気回路のシミュレーション（時定数が数桁異なるRCL回路）、化学反応ネットワーク（高速な中間反応と低速な全体反応の混在）、半導体デバイスのシミュレーション（キャリアの緩和時間と外部駆動の時間スケールの乖離）などが挙げられます。

安定性に加えて、実用的な数値計算では刻み幅を自動調整する機能が重要です。解が急変する領域では細かい刻み幅が必要で、緩やかに変化する領域では粗い刻み幅で十分だからです。

適応的刻み幅制御

動機

固定刻み幅では、解が急激に変化する区間で精度が不足し、緩やかに変化する区間で不必要に小さい刻み幅を使って計算が無駄になります。たとえば惑星探査機の軌道計算では、惑星への最接近時には極めて細かい刻み幅が必要ですが、巡航中は粗い刻み幅で十分です。適応的刻み幅制御は、局所誤差の推定に基づいて刻み幅を自動的に調整し、必要な精度を保ちながら計算コストを最小化します。

埋め込み型RK法

適応的刻み幅制御の鍵は、「現在の刻み幅での誤差がどの程度か」を効率的に見積もることです。ドルマン=プリンス法（Dormand-Prince, RK45）は、4次と5次のRK解を同時に計算し、その差を局所誤差の推定に使います。両方の方法が同じ $k_i$ を共有するため、追加の関数評価はほとんど不要です。

$$ e_{n+1} = \hat{y}_{n+1} – y_{n+1} = h\sum_{i=1}^{s} (b_i – \hat{b}_i) k_i $$

5次の方法 $\hat{y}_{n+1}$ を「真の値」と見なし、4次の方法 $y_{n+1}$ との差を誤差とします。

刻み幅の更新

局所誤差が許容値 $\varepsilon$ を超えたらステップをやり直し、刻み幅を縮小します。

$$ h_{\text{new}} = h \cdot \min\left(S \cdot \left(\frac{\varepsilon}{\|e_{n+1}\|}\right)^{1/5}, h_{\text{max}}\right) $$

安全係数 $S \approx 0.9$ は過度な変更を防ぎます。SciPyの solve_ivp(method='RK45') はこの手法を実装しています。

適応的刻み幅制御は「ユーザーが刻み幅を選ぶ必要がない」という実用上の大きな利点を持ちます。固定刻み幅の場合、小さすぎれば計算が無駄に遅くなり、大きすぎれば精度が不足するか安定性が破綻します。適切な刻み幅は問題によって異なり、しかも同じ問題でも時刻によって変わります。適応制御はこの選択を自動化し、ユーザーは許容誤差（rtol, atol）を指定するだけでよくなります。

理論を一通り押さえたところで、Pythonで各手法の挙動を実際に確認してみましょう。

Pythonによる実装

以下のコードでは、指数減衰 $y’ = -y$ とスティッフ問題 $y’ = -50(y – \cos t)$ を題材に、各手法の精度・収束次数・安定性を可視化します。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 数値解法の実装
def euler_forward(f, t_span, y0, h):
    t = np.arange(t_span[0], t_span[1] + h, h)
    y = np.zeros((len(t), len(np.atleast_1d(y0))))
    y[0] = np.atleast_1d(y0)
    for i in range(len(t) - 1):
        y[i+1] = y[i] + h * np.atleast_1d(f(t[i], y[i]))
    return t, y

def rk2_midpoint(f, t_span, y0, h):
    t = np.arange(t_span[0], t_span[1] + h, h)
    y = np.zeros((len(t), len(np.atleast_1d(y0))))
    y[0] = np.atleast_1d(y0)
    for i in range(len(t) - 1):
        k1 = np.atleast_1d(f(t[i], y[i]))
        k2 = np.atleast_1d(f(t[i] + h/2, y[i] + h/2 * k1))
        y[i+1] = y[i] + h * k2
    return t, y

def rk4_classic(f, t_span, y0, h):
    t = np.arange(t_span[0], t_span[1] + h, h)
    y = np.zeros((len(t), len(np.atleast_1d(y0))))
    y[0] = np.atleast_1d(y0)
    for i in range(len(t) - 1):
        k1 = np.atleast_1d(f(t[i], y[i]))
        k2 = np.atleast_1d(f(t[i] + h/2, y[i] + h/2 * k1))
        k3 = np.atleast_1d(f(t[i] + h/2, y[i] + h/2 * k2))
        k4 = np.atleast_1d(f(t[i] + h, y[i] + h * k3))
        y[i+1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4)
    return t, y

fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# --- (a) 精度比較: y' = -y, y(0) = 1 ---
ax = axes[0, 0]
f_test = lambda t, y: -y
y_exact = lambda t: np.exp(-t)
t_span = (0, 5)

for h, ls in [(0.5, "-"), (0.2, "--"), (0.1, ":")]:
    t_e, y_e = euler_forward(f_test, t_span, [1.0], h)
    ax.plot(t_e, y_e[:, 0], f"b{ls}", linewidth=1.5,
            label=f"Euler $h={h}$", alpha=0.7)

t_rk, y_rk = rk4_classic(f_test, t_span, [1.0], 0.5)
ax.plot(t_rk, y_rk[:, 0], "r-", linewidth=2, label="RK4 $h=0.5$")

t_exact = np.linspace(0, 5, 200)
ax.plot(t_exact, y_exact(t_exact), "k-", linewidth=2.5, label="Exact")

ax.set_xlabel("$t$", fontsize=12)
ax.set_ylabel("$y(t)$", fontsize=12)
ax.set_title("(a) $y' = -y$: Euler vs RK4", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# --- (b) 収束次数の確認 ---
ax = axes[0, 1]
hs = np.array([0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005])
errs_euler = []
errs_rk2 = []
errs_rk4 = []

for h in hs:
    t_e, y_e = euler_forward(f_test, (0, 1), [1.0], h)
    errs_euler.append(abs(y_e[-1, 0] - np.exp(-1)))

    t_m, y_m = rk2_midpoint(f_test, (0, 1), [1.0], h)
    errs_rk2.append(abs(y_m[-1, 0] - np.exp(-1)))

    t_r, y_r = rk4_classic(f_test, (0, 1), [1.0], h)
    errs_rk4.append(abs(y_r[-1, 0] - np.exp(-1)))

ax.loglog(hs, errs_euler, "bo-", linewidth=2, markersize=6, label="Euler (order 1)")
ax.loglog(hs, errs_rk2, "gs-", linewidth=2, markersize=6, label="RK2 (order 2)")
ax.loglog(hs, errs_rk4, "r^-", linewidth=2, markersize=6, label="RK4 (order 4)")

# 参考線
ax.loglog(hs, 0.5 * hs**1, "b:", linewidth=1, alpha=0.5, label="$O(h)$")
ax.loglog(hs, 0.2 * hs**2, "g:", linewidth=1, alpha=0.5, label="$O(h^2)$")
ax.loglog(hs, 0.05 * hs**4, "r:", linewidth=1, alpha=0.5, label="$O(h^4)$")

ax.set_xlabel("Step size $h$", fontsize=12)
ax.set_ylabel("Error at $t = 1$", fontsize=12)
ax.set_title("(b) Convergence Order Verification", fontsize=13)
ax.legend(fontsize=9, ncol=2)
ax.grid(True, alpha=0.3, which="both")

# --- (c) 安定性領域 ---
ax = axes[1, 0]
x_range = np.linspace(-4, 2, 400)
y_range = np.linspace(-3, 3, 400)
X, Y = np.meshgrid(x_range, y_range)
Z = X + 1j * Y

# 前進オイラー: |1 + z| < 1
R_euler = np.abs(1 + Z)

# RK4: |1 + z + z^2/2 + z^3/6 + z^4/24| < 1
R_rk4 = np.abs(1 + Z + Z**2/2 + Z**3/6 + Z**4/24)

# 後退オイラー: |1/(1-z)| < 1 → |1-z| > 1
R_be = np.abs(1 / (1 - Z))

ax.contour(X, Y, R_euler, levels=[1], colors="blue", linewidths=2)
ax.contourf(X, Y, R_euler, levels=[0, 1], colors=["lightblue"], alpha=0.3)

ax.contour(X, Y, R_rk4, levels=[1], colors="red", linewidths=2)
ax.contourf(X, Y, R_rk4, levels=[0, 1], colors=["lightyellow"], alpha=0.3)

ax.contour(X, Y, R_be, levels=[1], colors="green", linewidths=2, linestyles="--")

ax.plot(0, 0, "ko", markersize=4)
ax.axhline(0, color="gray", linewidth=0.5)
ax.axvline(0, color="gray", linewidth=0.5)

ax.text(-1, 0.3, "Euler", fontsize=11, color="blue", fontweight="bold")
ax.text(-2.8, 2.2, "RK4", fontsize=11, color="red", fontweight="bold")
ax.text(1.5, 2.2, "Backward\nEuler", fontsize=10, color="green", fontweight="bold")

ax.set_xlabel("Re($h\\lambda$)", fontsize=12)
ax.set_ylabel("Im($h\\lambda$)", fontsize=12)
ax.set_title("(c) Stability Regions", fontsize=13)
ax.set_aspect("equal")
ax.set_xlim(-4, 2)
ax.set_ylim(-3, 3)
ax.grid(True, alpha=0.3)

# --- (d) スティッフ問題での比較 ---
ax = axes[1, 1]

# スティッフ問題: y' = -50(y - cos(t))
def stiff_rhs(t, y):
    return -50 * (y - np.cos(t))

def stiff_exact(t):
    # 定常応答
    return (50 * np.cos(t) + np.sin(t)) / (50**2 + 1) * 50

# RK45（適応的）
sol_rk45 = solve_ivp(stiff_rhs, [0, 2], [0], method="RK45",
                      t_eval=np.linspace(0, 2, 500), rtol=1e-8)

# Radau（陰的、スティッフ向け）
sol_radau = solve_ivp(stiff_rhs, [0, 2], [0], method="Radau",
                       t_eval=np.linspace(0, 2, 500), rtol=1e-8)

# 前進オイラー（大きいh）
try:
    t_eu, y_eu = euler_forward(lambda t, y: np.array(stiff_rhs(t, y)),
                               (0, 2), [0.0], 0.05)
    ax.plot(t_eu, y_eu[:, 0], "b-", linewidth=1, alpha=0.6,
            label="Euler $h=0.05$ (unstable)")
except:
    pass

# 前進オイラー（安定なh）
t_eu2, y_eu2 = euler_forward(lambda t, y: np.array(stiff_rhs(t, y)),
                              (0, 2), [0.0], 0.01)
ax.plot(t_eu2, y_eu2[:, 0], "c-", linewidth=1, alpha=0.6,
        label="Euler $h=0.01$ (stable)")

ax.plot(sol_rk45.t, sol_rk45.y[0], "r-", linewidth=2, label="RK45 (adaptive)")
ax.plot(sol_radau.t, sol_radau.y[0], "g--", linewidth=2, label="Radau (implicit)")

ax.set_xlabel("$t$", fontsize=12)
ax.set_ylabel("$y(t)$", fontsize=12)
ax.set_title("(d) Stiff Problem: $y' = -50(y - \\cos t)$", fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("numerical_ode.png", dpi=150, bbox_inches="tight")
plt.show()

この4つの図から、ODE数値解法の精度と安定性に関する重要な知見が読み取れます。

1. 左上（精度比較）: 同じ刻み幅 $h = 0.5$ でも、前進オイラー法は真の解から大きく逸れていますが、RK4はほぼ厳密解と一致しています。オイラー法で同等の精度を得るには、はるかに小さい刻み幅（$h = 0.1$ 以下）が必要です

2. 右上（収束次数の検証）: 対数-対数プロットで、各手法の誤差が理論的な収束次数（オイラー: $O(h)$, RK2: $O(h^2)$, RK4: $O(h^4)$）の参考線と平行になっています。RK4の傾き4は、刻み幅を半分にすると誤差が $1/16$ になることを意味し、非常に効率的です

3. 左下（安定性領域）: 前進オイラー法（青）の安定性領域は中心 $(-1, 0)$, 半径1の円内部で最も小さく、RK4（赤）は実軸上で約 $[-2.8, 0]$ まで広がっています。後退オイラー法（緑の破線）の安定性領域は原点 $(1, 0)$ の円の外部で、左半平面全体を含む（A安定）ことが確認できます

4. 右下（スティッフ問題）: $\lambda = -50$ のスティッフ問題に対して、前進オイラー法（$h = 0.05$）は安定性条件 $h < 2/50 = 0.04$ を満たさず振動しています。$h = 0.01$ では安定ですが、物理的に意味のある時間スケールに比べて不必要に小さい刻み幅を使わざるを得ず非効率です。適応的RK45とRadau法は自動的に適切な刻み幅を選び、効率よく解を計算しています。この比較から、スティッフ問題ではA安定な陰解法が不可欠であることが実感できます

手法選択のガイドライン

実際の問題でどの手法を選ぶべきかの実践的な指針をまとめます。

まとめ

本記事では、常微分方程式の数値解法を精度と安定性の観点から体系的に解説しました。

• 前進オイラー法は1次の精度で、概念的に重要だが実用には不十分。テイラー展開から自然に導出され、すべての高次手法の出発点となる
• ルンゲ=クッタ法は区間内で $f$ を複数回評価し、高精度を達成する。特にRK4は4段4次の精度で広く使われており、「次数の壁」を考えると効率の良い方法である
• 安定性領域は $h\lambda$ の複素平面上で定義され、陽解法の安定性には刻み幅の上限がある。安定性条件を満たさない計算は、精度以前に結果自体が無意味になる
• スティッフ問題では陰解法（後退オイラー法、Radau法）やA安定な方法が不可欠。異なる時間スケールが混在する系は工学的に非常に多い
• 適応的刻み幅制御は埋め込み型RK法による局所誤差推定に基づき、効率と精度を両立する。実用上は solve_ivp のような適応型ソルバーを使うのが最も安全

ODE数値解法は、数値解析の中でも最も実用的な分野の一つです。衛星軌道から化学反応まで、自然界の動的な現象を理解し予測するための不可欠な道具となっています。

次のステップとして、以下の記事も参考にしてください。

• 連立常微分方程式の理論 — 連立系の解析的理論
• 常微分方程式の安定性解析 — 解析的な安定性理論
• 数値線形代数の基礎 — 陰解法での連立方程式
• 数値積分の理論 — 数値積分との関連