サイトアイコン アマチュア無線局JS2IIU

Pythonで攻略する流体力学の数値計算 第8回：常微分方程式の解法(1) ― 初期値問題と数値積分

JS2IIU

こんにちはJS2IIUです。Pythonで攻略する流体力学シミュレーション連載の第8回です。

前回は、数値計算の安定性と精度、そして計算が爆発しないための絶対的なルールであるCFL条件について学びました。計算の網目（格子）と情報の伝播速度の関係がいかに重要か、実感いただけたのではないでしょうか。

今回からは、ナビエ・ストークス方程式の時間発展、すなわち「次の瞬間に流体がどう動くか」を計算するための具体的な手法へと踏み込みます。その基礎となるのが、常微分方程式（Ordinary Differential Equation, ODE）の初期値問題を解くための数値積分アルゴリズムです。

流体の方程式は非常に複雑な偏微分方程式ですが、時間を離散化して一歩ずつ進めるという視点に立てば、各時刻における状態の更新は常微分方程式を解くプロセスに他なりません。今回は、最もシンプルな「オイラー法」から、実用的な「ルンゲ・クッタ法」まで、Pythonを用いた実装とともに解説します。今回もよろしくお願いします。

1. 初期値問題とは何か

物理現象の多くは、「現在の状態」が分かれば「わずかな時間後の状態」を予測できる、という形で記述されます。数学的には、従属変数 \(y\) の時間変化率（微分）が、現在の時刻 \(t\) と状態 \(y\) の関数として与えられる形式です。

$$ \frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0 $$

ここで \(y(t_0) = y_0\) が「初期値」です。この式を解いて、任意の時刻 \(t\) における \(y(t)\) を求めることが目標です。流体計算に当てはめるなら、\(y\) は各格子点における「流速」や「圧力」、\(f\) は前回の差分法で学んだ空間微分項や圧力勾配項に相当します。

この微分方程式を積分によって解くため、このプロセスを数値積分と呼びます。
物理現象の多くは、「今この瞬間の状態」が分かれば、「ほんの少し先の未来の状態」を予測できるという性質を持っています。たとえば、振り子の運動、惑星の軌道、流体の流れなど、自然界のさまざまな現象は「現在の状態」から「次の状態」へと時間発展していきます。

このような現象を数学的に記述するのが初期値問題（Initial Value Problem, IVP）です。初期値問題とは、ある従属変数 \(y\) の時間変化率（微分）が、現在の時刻 \(t\) と状態 \(y\) の関数として与えられ、さらに「最初の状態（初期値）」が指定されている問題を指します。

一般的な形は次の通りです。

$$
\frac{dy}{dt} = f(t, y), \quad y(t_0) = y_0
$$

ここで、

この初期値問題を解くとは、「与えられた初期値から出発して、任意の時刻 \(t\) における \(y(t)\) を求める」ことを意味します。

流体計算の文脈で考えると、\(y\) は各格子点における「流速」や「圧力」などの物理量、\(f\) は前回の差分法で学んだ空間微分項や圧力勾配項など、流体の運動を決める要素に相当します。

実際の物理現象は解析的に解けることは稀なので、コンピュータ上ではこの微分方程式を「時間ごとに一歩ずつ」積分していく必要があります。この「微分方程式を数値的に積分して未来を予測する」プロセスが数値積分（numerical integration）です。

数値積分は、初期値問題を解くための最も基本的かつ重要な道具であり、流体シミュレーションをはじめとする多くの科学技術計算の根幹をなしています。

2. オイラー法：直線で未来を予測する

最も基本的でシンプルな数値積分の手法が、第6回や第7回の実験でも登場したオイラー法（Euler method）です。

オイラー法の考え方はとても直感的です。時刻 \(t_n\) における変化率（勾配）\(f(t_n, y_n)\) が、次のステップ \(t_{n+1} = t_n + \Delta t\) までずっと変わらず一定であると仮定し、そのまま直線的に未来を予測します。

図形的に言えば、「現在地での傾き（接線）」をそのまま延長して、次の時刻の値を決めるイメージです。

数式で表すと、
$$
y_{n+1} = y_n + f(t_n, y_n) \, \Delta t
$$
となります。

この式は、「今の値 \(y_n\) に、変化率 \(f(t_n, y_n)\) を \(\Delta t\) だけ掛けた分だけ進める」という非常に単純なものです。

オイラー法の特徴と直感

メリットとデメリット

このように、オイラー法は「簡単だが精度や安定性には限界がある」手法です。実際のシミュレーションでは、計算コストと精度・安定性のバランスを考えて、より高次の手法（例：ルンゲ・クッタ法）と使い分けることが重要です。

3. ルンゲ・クッタ法：複数の地点を「下見」して精度を高める

オイラー法の「直線的な予測」では、関数が曲線的に変化する場合に大きな誤差が生じやすいという欠点がありました。これを克服するために考案されたのが、ルンゲ・クッタ法（Runge-Kutta method）です。

特に「4次のルンゲ・クッタ法（RK4）」は、精度と計算コストのバランスが非常に優れており、科学技術計算の世界標準として広く使われています。

ルンゲ・クッタ法（RK4）の考え方

RK4の本質は、「未来を予測する際に、途中の複数の地点で“下見”をして、その平均的な勾配を使って進む」ことです。

オイラー法が「今の勾配だけで一直線に進む」のに対し、RK4は「途中で何度も勾配を測り直し、より曲線の形に沿った進み方をする」イメージです。

4つの“下見”の意味

1ステップの間に、次の4つの勾配（\(k_1, k_2, k_3, k_4\)）を計算します。

  1. \(k_1 = f(t_n, y_n)\)：現在地での勾配（オイラー法と同じ）
  2. \(k_2 = f(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} k_1)\)：\(k_1\) を使って半歩進んだ仮の地点での勾配
  3. \(k_3 = f(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} k_2)\)：さらにその仮の地点での勾配
  4. \(k_4 = f(t_n + \Delta t, y_n + \Delta t k_3)\)：1ステップ先まで進んだ仮の地点での勾配

これら4つの勾配を「中間の勾配を重視」するように重み付き平均し、次の値を決定します。

$$
y_{n+1} = y_n + \frac{\Delta t}{6} (k_1 + 2k_2 + 2k_3 + k_4)
$$

この重み付けは、曲線の形状をより正確に捉えるための工夫です。特に(k_2)や(k_3)の「中間点での勾配」を重視することで、関数の曲がり具合をしっかり反映できます。

精度・安定性・実用上の意義

このように、RK4は「複数の下見で曲線の形をしっかり捉え、直線的なオイラー法よりもはるかに信頼性の高い未来予測ができる」手法です。流体シミュレーションをはじめ、幅広い分野で標準的に使われています。

RK4のアイデアは、「一気に \(\Delta t\) 先を予測するのではなく、途中の地点で何度か勾配を下見して、その平均的な勾配を使って進む」というものです。

具体的には、1ステップの間に4つの勾配（\(k_1, k_2, k_3, k_4\)）を計算します。

  1. \(k_1 = f(t_n, y_n)\) ：現在の地点の勾配
  2. \(k_2 = f(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} k_1)\) ：\(k_1\) を使って半歩進んだ地点の勾配
  3. \(k_3 = f(t_n + \frac{\Delta t}{2}, y_n + \frac{\Delta t}{2} k_2)\) ：\(k_2\) を使って半歩進んだ地点の勾配
  4. \(k_4 = f(t_n + \Delta t, y_n + \Delta t k_3)\) ：\(k_3\) を使って1歩進んだ地点の勾配

最後に、これらを「中間の勾配を重視」するように重み付き平均して、次の値を決定します。

$$ y_{n+1} = y_n + \frac{\Delta t}{6} (k_1 + 2k_2 + 2k_3 + k_4) $$

RK4は4次精度、つまり時間刻み \(\Delta t\) を半分にすると、誤差は \(1/2^4 = 1/16\) にまで減少します。オイラー法とは比較にならないほど強力な武器です。

4. Pythonによる実装：振り子の振動で精度を比較する

それでは、PythonとPyTorchを用いて、オイラー法とルンゲ・クッタ法の性能を比較してみましょう。

題材として、「単振動（振り子）」を扱います。この系はエネルギーが保存されるため、数値誤差によってエネルギーが増減する様子（解の軌跡がズレる様子）が視覚的に非常に分かりやすいからです。

状態変数 \(\mathbf{y} = [\theta, \omega]^T\) （角度と角速度）とし、次の方程式を解きます。
$$ \frac{d\theta}{dt} = \omega, \quad \frac{d\omega}{dt} = -\sin(\theta) $$

Python
import torch
import matplotlib.pyplot as plt

def f(t, state):
    """
    振り子の方程式 (d_theta/dt = omega, d_omega/dt = -sin(theta))
    state: [theta, omega]
    """
    theta, omega = state[0], state[1]
    d_theta = omega
    d_omega = -torch.sin(theta)
    return torch.stack([d_theta, d_omega])

def euler_step(f, t, y, dt):
    """オイラー法による1ステップの更新"""
    return y + f(t, y) * dt

def runge_kutta_4_step(f, t, y, dt):
    """4次ルンゲ・クッタ法による1ステップの更新"""
    k1 = f(t, y)
    k2 = f(t + dt/2, y + dt/2 * k1)
    k3 = f(t + dt/2, y + dt/2 * k2)
    k4 = f(t + dt, y + dt * k3)
    return y + (dt/6) * (k1 + 2*k2 + 2*k3 + k4)

# --- シミュレーション設定 ---
dt = 0.2  # 比較のため、あえて粗い時間刻みを設定
t_max = 20.0
steps = int(t_max / dt)

# 初期条件: 角度1.0 rad, 角速度0
y0 = torch.tensor([1.0, 0.0])

# 結果の格納用
res_euler = [y0]
res_rk4 = [y0]

y_e, y_rk = y0.clone(), y0.clone()
t = 0.0

# --- 時間進展計算 ---
for _ in range(steps):
    y_e = euler_step(f, t, y_e, dt)
    y_rk = runge_kutta_4_step(f, t, y_rk, dt)
    t += dt
    res_euler.append(y_e)
    res_rk4.append(y_rk)

# テンソルに変換して可視化の準備
res_euler = torch.stack(res_euler)
res_rk4 = torch.stack(res_rk4)

# --- 可視化: 相空間プロット (theta vs omega) ---
plt.figure(figsize=(10, 6))
plt.plot(res_euler[:, 0].numpy(), res_euler[:, 1].numpy(), 'r--', label='Euler Method')
plt.plot(res_rk4[:, 0].numpy(), res_rk4[:, 1].numpy(), 'b-', label='Runge-Kutta 4')

# 理論的なエネルギー保存ライン (参考)
plt.title("Phase Space: Pendulum Swing (dt=0.2)")
plt.xlabel("Theta (Angle)")
plt.ylabel("Omega (Angular Velocity)")
plt.legend()
plt.grid(True)
plt.axis('equal')
plt.show()

コードの解説と結果の分析

  1. 関数の定義:
    f(t, state) では、現在の状態から微分値を計算して返します。PyTorchのテンソル計算を用いることで、多変量の方程式もシンプルに記述できます。
  2. 更新アルゴリズム:
    euler_step と runge_kutta_4_step を見比べてください。ルンゲ・クッタ法は4回の関数評価を行っていますが、その構造自体は非常に整理されており、Pythonで容易に実装できることがわかります。
  3. 相空間の軌跡: 実行結果のグラフを見てみましょう。本来、この振り子の軌跡は「閉じたループ（楕円に近い形）」を描き、同じところを回り続けるはずです。
    • オイラー法（赤の破線）: 軌跡が外側へどんどん広がっています。これは「数値的なエネルギー増大」が発生している証拠で、時間が経つと振り子が激しく回転し始めるという物理的に誤った結果になります。
    • ルンゲ・クッタ法（青の実線）: 同じ時間刻み \(dt=0.2\) でも、ほぼ完全に元のループを維持しています。非常に安定しており、信頼性が高いことが視覚的にわかります。

5. 流体シミュレーションへの応用

今回学んだ数値積分は、流体力学の核心部分でどのように使われるのでしょうか。

例えば、第4回で紹介したナビエ・ストークス方程式の時間発展項を考えてみましょう。
$$ \frac{\partial \mathbf{u}}{\partial t} = \text{ (移流項) } + \text{ (拡散項) } + \text{ (圧力勾配項) } $$

右辺をまとめて \(F(\mathbf{u})\) と置けば、これは \(\frac{d\mathbf{u}}{dt} = F(\mathbf{u})\) という巨大な常微分方程式とみなせます。

多くの単純な流体計算（例えば第13回で扱う「流れ関数-渦度法」など）では、計算コストを優先してオイラー法が使われることもあります。しかし、気象予測や航空宇宙などの高精度な解析では、今回学んだルンゲ・クッタ法や、さらに多段階の予測を行う「アダムス・バシュフォース法」などが活躍しています。

6. まとめと次回の予告

第8回では、シミュレーションの「時間」を制御するための強力なツールを手にしました。

今回の「初期値問題」は、時間に沿って変化を追いかける手法でした。しかし、流体の世界には「時間が経っても変化しない、落ち着いた状態」を知りたい場合や、「空間全体のバランスを一度に解きたい」場合もあります。

次回、第9回では、もう一つの重要な微分方程式の解き方、「境界値問題」について解説します。行列計算を駆使して空間の分布を一気に決定する手法を学び、定常的な物理現象の解法を習得していきましょう。

最後まで読んでいただきありがとうございます。

モバイルバージョンを終了