こんにちはJS2IIUです。Pythonで攻略する流体力学シミュレーション連載の第7回です。
前回、第6回では、微分という連続的な数学概念を、コンピュータが扱える離散的な「引き算と割り算」に翻訳する差分法の基礎を学びました。テイラー展開を用いて、前進差分、後退差分、そしてより精度の高い中心差分といった手法を導出し、格子間隔が計算精度にどう影響するかを実験しました。
しかし、微分を正しく離散化しただけでは、シミュレーションが成功するとは限りません。プログラムを走らせた瞬間、計算値が指数関数的に増大し、画面が「NaN（Not a Number）」で埋め尽くされる――そんな「計算の爆発（発散）」は、数値計算を学ぶ誰もが一度は通る道です。
今回は、シミュレーションが「正しく解に近づく」ための数学的な保証であるLAX（ラックス）の同値定理と、爆発を防ぐための絶対的なルールであるCFL条件について詳しく解説します。今回もよろしくお願いします。
1. 「正しい」シミュレーションとは何か？ ― 収束性
数値計算における「正しい」シミュレーションとは、計算によって得られる数値解が、現実の物理現象を記述する「真の解」にどれだけ近づくか、という観点で評価されます。ここで重要なのが収束性（Convergence）という概念です。
収束性とは、格子間隔（(\Delta x)）や時間ステップ（(\Delta t)）を限りなく小さくしていったとき、数値解が理論的な真の解に近づいていく性質を指します。つまり、計算を細かくすればするほど、シミュレーションの結果が現実に忠実になることが理想です。
しかし、実際のシミュレーションで「収束しているかどうか」を直接証明するのは非常に難しい問題です。なぜなら、真の解そのものが未知である場合が多く、また計算機の有限精度や様々な誤差が影響するためです。
この難問に対して、数値解析の巨匠ピーター・ラックス（Peter Lax）は、収束性をより現実的に判定するための2つの基準を提案しました。それが有名なLAXの同値定理です。
LAXの同値定理
LAXの同値定理は、特に線形な初期値問題に対して、次のように述べられます。
線形な初期値問題において、適切に構成された差分近似が収束するための必要十分条件は、その近似が適合性と安定性の両方を満たすことである。
$$
ext{収束性} \iff \text{適合性} + \text{安定性}
$$
この定理は、数値流体力学（CFD）をはじめとする数値シミュレーション分野で、計算手法の信頼性を判断するための最重要指針となっています。
それぞれの用語について、もう少し丁寧に説明します。
- 適合性（Consistency）:
- 差分法などで連続的な微分方程式を離散化したとき、その離散化方程式が、格子間隔 \(\Delta x\) や時間刻み \(\Delta t\) をゼロに近づけた極限で、元の微分方程式と一致する性質です。
- これは、前回学んだテイラー展開を使って、離散化による切断誤差（truncation error）が十分小さくなる、すなわちゼロに収束することを確認することで保証できます。
- 適合性があるということは、「計算式そのものが物理法則を正しく表現している」ことを意味します。
- 安定性（Stability）:
- 実際の計算では、丸め誤差や初期値のわずかなズレなど、様々な小さな誤差が必ず発生します。
- 安定性とは、こうした誤差が計算を進めるごとに増幅されて暴走することなく、抑えられている（または消えていく）性質を指します。
- 安定性がないと、どんなに正しい式を使っても、計算結果が発散（爆発）してしまい、物理的に意味のない値になってしまいます。
このように、
\(\text{適合性}\)（正しい式）と\(\text{安定性}\)（壊れない計算）の両方が揃って初めて、\(\text{収束性}\)（真の解に近づく）が保証される、というのがLAXの同値定理の本質です。
言い換えれば、「正しい数式を使っていても、計算が安定していなければ意味がない」「逆に、計算が安定でも、そもそも式が間違っていれば正しい解には辿り着けない」ということです。
この三位一体の考え方は、数値シミュレーションの信頼性を考える上で、常に立ち返るべき基本原則となります。
2. 計算の爆発を回避する ― CFL条件
安定性を考える上で絶対に避けて通れないのが、CFL条件（Courant-Friedrichs-Lewy Condition）です。CFL条件は、数値計算において「計算が爆発しない」ための物理的・数学的な制約であり、特に流体シミュレーションや波動現象の計算で極めて重要な役割を果たします。
物理的なイメージ
まず直感的なイメージから説明します。
例えば、ある場所で波が発生し、その波が速さ \(a\) で右方向に伝わっていくとします。計算機上では、空間を格子点（\(\Delta x\) 間隔）で区切り、時間も一定のステップ（\(\Delta t\)）で進めていきます。
このとき、1ステップ（\(\Delta t\)）の間に波がどれだけ進むかを考えます。もし波が1ステップで隣の格子点（\(\Delta x\)）に到達する前に、次の計算（時間更新）が行われれば、計算機は「情報の伝播」を正しく追いかけることができます。
しかし、もし波が1ステップで格子一つ分よりも遠くへ進んでしまう（つまり、\(a \Delta t > \Delta x\)）と、計算機はその情報の通過を捉えきれず、物理的な因果関係が壊れてしまいます。その結果、計算値が暴走し、現実とはかけ離れた「爆発」や「発散」が起こります。
この「情報の伝播速度」と「計算の進め方」のバランスを守るための条件が、CFL条件です。
数式での表現と意味
1次元の移流方程式（\(\frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0\)）を例にとると、CFL条件は次のような無次元数（クーラン数 \(C\)）で表されます。
$$
C = \frac{a \Delta t}{\Delta x}
$$
このクーラン数 \(C\) は、「1ステップで情報がどれだけ格子を進むか」を表す指標です。
- \(C \le 1\): 情報の伝わる速さに対して、十分に細かく時間を区切っている状態。計算は安定し、物理現象を正しく追従できます。
- \(C > 1\): 情報の伝わる速さが計算の網目（格子）を追い越してしまい、計算が不安定になり爆発・発散が起こります。
このCFL条件は、数値計算の安定性を守るための「安全装置」とも言えます。特に陽的（explicit）な時間積分法では、\(C\)が1を超えないように、\(\Delta t\)や\(\Delta x\)を慎重に選ぶ必要があります。
なぜCFL条件が必要なのか？
現実の物理現象では、情報や影響は有限の速度で伝わります。数値計算でも、その物理的な因果律を守る必要があります。CFL条件を破ると、計算機は「まだ伝わっていないはずの情報」を先回りして計算してしまい、現実とは異なる不自然な結果を生み出します。
例えば、波が1ステップで2つ先の格子点まで進んでしまうと、途中の格子点の情報を無視してしまうことになり、計算の論理が破綻します。これが「計算の爆発」の根本原因です。
実際のシミュレーション設計での使い方
シミュレーションを設計する際は、まず物理現象の伝播速度（\(a\)）を見積もり、空間格子幅（\(\Delta x\)）を決めた上で、CFL条件を満たすように時間ステップ（\(\Delta t\)）を設定します。
$$
\Delta t \leq \frac{\Delta x}{a}
$$
このガイドラインを守ることで、計算が爆発せず、物理現象を忠実に再現できる安定したシミュレーションが実現できます。
CFL条件は、数値流体力学や波動方程式の計算だけでなく、さまざまな分野のシミュレーションで「計算の信頼性」を支える最も基本的なルールの一つです。安定性を確保するために、必ず意識しておくべき指標です。
3. 精度（Accuracy）の再確認
安定して計算が回るようになったら、次に重要になるのが「どれだけ正確に真の解に近づいているか」、すなわち精度（Accuracy）です。精度は、数値計算の品質を左右する大きな要素であり、単に「爆発しない」だけでなく「どれだけ誤差が小さいか」を評価する指標です。
精度の「次数」とは？
数値計算における精度は、主に「次数（order）」という概念で表されます。これは、格子間隔 \(\Delta x\) や時間刻み \(\Delta t\) を小さくしたとき、誤差がどのように減少するかを示すものです。
- 1次精度（一次精度）: 誤差が \(\Delta x\) や \(\Delta t\) に比例して減少します。例えば、格子幅を半分にすると誤差も約半分になります。
- 2次精度（二次精度）: 誤差が \(\Delta x^2\) や \(\Delta t^2\) に比例して減少します。格子幅を半分にすると誤差は4分の1に減ります。
この「次数」が高いほど、同じ計算コストでより高精度な結果が得られるため、効率的なシミュレーションが可能になります。
誤差の種類と精度の意味
数値計算の誤差には主に2種類あります。
- 切断誤差（truncation error）: 微分方程式を差分方程式に置き換える際に生じる理論的な誤差。テイラー展開で評価できます。
- 丸め誤差（round-off error）: 計算機の有限精度による誤差。
精度の「次数」は主に切断誤差の減少速度を表します。高次のスキームを使うことで、同じ格子幅でもより小さな誤差で計算できます。
高次スキームの利点と注意点
2次精度以上の高次スキームは、格子を細かくしなくても高精度な結果が得られるため非常に強力です。例えば、2次精度なら格子幅を半分にするだけで誤差が4分の1に減ります。
しかし、高次スキームには「数値的な振動（ギブス現象）」が発生しやすい、境界条件の扱いが難しい、安定性が低下しやすいなどの注意点もあります。高精度と高安定性は必ずしも両立しないため、現実のシミュレーションでは「精度」と「安定性」のバランスを取ることが重要です。
精度と安定性のトレードオフ
高精度なスキームを使えば誤差は減りますが、安定性が損なわれる場合もあります。逆に、安定性を重視しすぎると精度が犠牲になることもあります。この「精度と安定性のトレードオフ」は、数値シミュレーション設計の永遠のテーマです。
今後の連載では、さまざまなスキームの精度・安定性の違いや、どのようにバランスを取るかについても詳しく解説していきます。
- 1次精度: 誤差が \(\Delta x\) や \(\Delta t\) に比例して減る。
- 2次精度: 誤差が \(\Delta x^2\) や \(\Delta t^2\) に比例して減る。
2次精度のスキームは、格子を半分にするだけで誤差が4分の1になるため非常に強力ですが、一方で「数値的な振動」を起こしやすいといった気難しい側面もあります。今後の連載では、この「安定性」と「精度」の絶妙なトレードオフを何度も扱うことになります。
4. Pythonによる実装：CFL条件と安定性の実験
それでは、PyTorchを用いて、CFL条件を満たす場合と満たさない場合で計算結果がどう変わるかを視覚化してみましょう。
題材は、最もシンプルな「1次元線形移流方程式」です。
$$ \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0 $$
この式は、ある分布 \(u\) が形を変えずに速さ \(a\) で右へ移動していく現象を表します。
import torch
import matplotlib.pyplot as plt
def advection_simulation(cfl, nx=100, nt=50):
"""
1次元移流方程式を計算し、安定性を検証する
cfl: 指定するクーラン数 (a * dt / dx)
"""
# 領域の設定
L = 2.0
dx = L / (nx - 1)
a = 1.0 # 移流速度
# CFL条件からdtを逆算
dt = cfl * dx / a
# 初期条件: 矩形波
x = torch.linspace(0, L, nx)
u = torch.zeros(nx)
u[(x >= 0.5) & (x <= 1.0)] = 1.0
u_init = u.clone()
# 時間発展 (1次後退差分スキーム)
# u[n+1, i] = u[n, i] - a * dt/dx * (u[n, i] - u[n, i-1])
for n in range(nt):
u_next = u.clone()
# ベクトル演算で一気に計算 (境界条件は周期境界とする)
u_next[1:] = u[1:] - cfl * (u[1:] - u[:-1])
u_next[0] = u[0] - cfl * (u[0] - u[-1])
u = u_next
# 数値が極端に大きくなったら爆発とみなして中断
if torch.any(torch.abs(u) > 1e5):
print(f"Calculation exploded at step {n} (CFL={cfl})")
return x, u_init, u, True
return x, u_init, u, False
# --- 実験: 安定なケースと不安定なケース ---
x_stable, u_init, u_stable, _ = advection_simulation(cfl=0.8) # CFL <= 1
x_unstable, _, u_unstable, exploded = advection_simulation(cfl=1.1) # CFL > 1
# --- 可視化 ---
plt.figure(figsize=(12, 5))
# 安定な計算 (CFL = 0.8)
plt.subplot(1, 2, 1)
plt.plot(x_stable.numpy(), u_init.numpy(), '--', label='Initial')
plt.plot(x_stable.numpy(), u_stable.numpy(), label='Final (Stable)')
plt.title(f"Stable Case (CFL = 0.8)")
plt.xlabel("x")
plt.ylabel("u")
plt.legend()
plt.grid(True)
# 不安定な計算 (CFL = 1.1)
plt.subplot(1, 2, 2)
plt.plot(x_unstable.numpy(), u_init.numpy(), '--', label='Initial')
if exploded:
plt.text(0.5, 0.5, "EXPLODED (NaN/Inf)", color='red', fontsize=12, ha='center')
else:
plt.plot(x_unstable.numpy(), u_unstable.numpy(), label='Final (Unstable)', color='red')
plt.title(f"Unstable Case (CFL = 1.1)")
plt.xlabel("x")
plt.ylabel("u")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()
コードの解説と出力結果の分析
- CFL数による時間刻みの決定:
コード内では、まず目標とするcflを決め、そこから物理法則に則ってdtを計算しています。これが、実際のシミュレーションを設計する際の標準的な手順です。 - 後退差分スキーム:
ここでは「情報の流れてくる方向（上流）」の点を使って微分を計算する「上流差分」を採用しています。 - 安定なケース (CFL=0.8):
左側のグラフを見てください。矩形波が崩れることなく（少しなまっていますが）、右へ移動しているのがわかります。誤差が増幅されることなく、計算が制御下にある状態です。 - 不安定なケース (CFL=1.1):
右側のグラフ（あるいは「EXPLODED」の警告）に注目してください。CFL条件をわずか0.1超えただけで、計算値は数ステップのうちに巨大な数値へと跳ね上がり、波の形を全く留めなくなります。これが「計算の爆発」の正体です。
5. 信頼性の高いシミュレーションへの第一歩
今回の実験で見た通り、差分方程式という「数学的な翻訳」が正しくても、パラメータ選びという「運用」を誤れば、シミュレーションは容易に崩壊します。
LAXの同値定理が教えてくれるのは、「計算が爆発していないか（安定性）」を確認することこそが、正しい答え（収束性）を得るための近道であるということです。
数値流体力学の世界では、この他にも「非線形性の問題」や「拡散項の影響」など、安定性を左右する要因が数多く存在します。しかし、今回の「情報の伝達速度を超えてはいけない（CFL条件）」という基本ルールは、そのすべての基礎となるものです。
6. まとめと次回の予告
第7回では、数値シミュレーションの信頼性を支える理論を学びました。
- 収束・適合・安定の三位一体を定義するLAXの同値定理。
- 情報の伝播速度と格子サイズの物理的制約であるCFL条件。
- 計算の爆発（発散）と、それを防ぐパラメータ設計の重要性。
- Pythonによる安定・不安定の比較実験。
これで、皆さんは「爆発しないシミュレーション」を設計するための理論的な武器を手に入れました。
しかし、ここまでの計算はすべて「時間ステップ」を一定のルールで進めるものでした。より複雑な時間変化を、より効率的に、より正確に追いかけるにはどうすればよいでしょうか？
次回、第8回からは、時間発展の計算手法をさらに深掘りします。常微分方程式の解法として知られる「オイラー法」や「ルンゲ・クッタ法」をPythonで実装し、流体シミュレーションの時間進展計算におけるその強力な性能を体験していきましょう。
最後まで読んでいただきありがとうございます。
コメント