FTCSが震えるとき、風上は鈍る — 1次移流スキーム、人工粘性、CFLの等価交換

中心差分が爆発し風上差分が生き残る理由、そしてその代償

1960年代初頭、John von Neumannの同僚たちは、最も単純なPDE — 1次元線形移流 — に最も自然な中心差分を組み込みました。結果は暴走でした。静かに流れるはずの矩形パルスが後ろの縁から振動を膨らませ、数ステップで格子の外へ飛び出します。この記事では、その暴走がどこから来るのか、なぜたった一格子の風上差分がそれを鎮める代わりに絵をぼやけさせるのか、そして格子幅と時間幅のあいだに引かれた見えない線(CFL条件)が何を語るのかを追っていきます。最後にPython 50行と対話デモで、3つのスキームを同じ格子の上で競争させます。

最も自然な差分が真っ先に崩れる#

1次元線形移流方程式はたった一行です。

\partial_t q + u\,\partial_x q = 0

ここで $q$ は運ばれる量、 $u>0$ は一定の速度です。解は単純で、初期形状が速度 $u$ でそのまま平行移動します。

時間に前進Euler、空間に中心(2次)差分を使うと、更新規則は次のように出てきます。

q_i^{n+1} = q_i^{n} - \tfrac{c}{2}\,\big(q_{i+1}^{n} - q_{i-1}^{n}\big)

$c \equiv u\,\Delta t / \Delta x$ はCourant数(1ステップで波が何セル進むか)です。空間2次・時間1次なのでLTE(局所打切り誤差)は $\mathcal{O}(\Delta t,\Delta x^2)$ 。それでもこのスキームは $c$ がどれほど小さくても発散します。

理由はvon Neumann安定解析の一行で見えます。 $q^n_i = g^n e^{i k x_i}$ を代入すると増幅因子は

g = 1 - i\,c\,\sin(k\Delta x)

となり、 $|g|^2 = 1 + c^2 \sin^2(k\Delta x) > 1$ 。すべてのフーリエモードが毎ステップ成長します。 $|g|>1$ のモードがひとつでもあれば全体が壊れます。

情報は片側からだけ流れてくる#

問題の本質は数式ではなく因果性です。 $u>0$ のとき、時刻 $t_{n+1}$ の $q_i$ が依存すべき情報は上流側 ( $x<x_i$ ) からのみ来るべきです。ところが中心差分は下流の格子値 $q_{i+1}^n$ を等しい重みで取り込みます。流れがまだ知らない未来から情報を先取りするのです。

風上差分は因果性を強制します。

q_i^{n+1} = q_i^{n} - c\,\big(q_i^{n} - q_{i-1}^{n}\big), \quad u>0

すると今度はvon Neumannが次を教えてくれます。

g = 1 - c + c\,e^{-i k \Delta x}

$|g|^2 = 1 - c(1-c)\cdot 4\sin^2(k\Delta x/2)$ 。 $0 \le c \le 1$ なら $|g|\le 1$ 。すべてのモードが減衰するか変わりません。

風上 = 中心差分 + 隠れた粘性#

風上が安定な理由をもっと直感的に見ましょう。1次風上差分を恒等式で分解します。

\frac{q_i - q_{i-1}}{\Delta x} = \frac{q_{i+1} - q_{i-1}}{2\Delta x} - \frac{\Delta x}{2}\,\frac{q_{i+1} - 2q_i + q_{i-1}}{\Delta x^2}

右辺第1項は中心差分、第2項はちょうど2階微分の離散形です。つまり1次風上は

\partial_t q + u\,\partial_x q = \frac{u\,\Delta x}{2}\,\partial_x^2 q

という修正方程式(modified equation)を解いていることになります。私たちは移流を解こうとしましたが、アルゴリズムは知らず知らず粘性 $\nu_{\text{num}} = u\Delta x/2$ の拡散をおまけで付けてくれたのです。これが人工粘性(artificial viscosity, または numerical diffusion)です。

この人工粘性が暴走を抑えます。ただし無料ではありません。矩形パルスの角は毎ステップごとに広がり、幅 $\sqrt{4\nu_{\text{num}} t}$ で滲みます。格子を細かく( $\Delta x$ 小)すれば粘性は減りますが、同じ最終時刻に到達するにはより多くのステップが必要です。両者がちょうど釣り合うので、滲みは決して完全には消えません。

CFL — 格子が追いつけなくなる瞬間#

風上が $0\le c \le 1$ でしか安定でないことには幾何学的な理由があります。時刻 $t_{n+1}$ の $q_i$ は、時刻 $t_n$ に位置 $x_i - u\Delta t$ にあった値です。その点が $x_{i-1}$ と $x_i$ の間にあれば、2つの格子値で補間できます。 $u\Delta t > \Delta x$ になると、その点は $x_{i-2}$ より向こうへ行きます — アルゴリズムのステンシルはそこまで見ていません。

これがCourant–Friedrichs–Lewy条件の核心です。

\Delta t \le \frac{\Delta x}{|u|}

数値的依存域(ステンシル)は物理的依存域(特性線)を必ず含まなければなりません。CFLは陽的スキームすべてに対する必要条件であって十分条件ではありません。FTCSはCFLを満たしていても爆発します — 原因はステンシル不足ではなく因果性違反だからです。

NumPy 50行で追う3スキーム#

同じ格子、同じ初期条件、同じCFLで3スキームをそのまま回します。コードは意図的に短くしています。

import numpy as np
 
N = 200
dx = 1.0 / N
u = 1.0
cfl = 0.8
dt = cfl * dx / u
n_steps = 400
 
x = (np.arange(N) + 0.5) * dx
q0 = ((x > 0.2) & (x < 0.4)).astype(float)  # 矩形パルス
 
def ftcs_advect(q, c):
    return q - 0.5 * c * (np.roll(q, -1) - np.roll(q, 1))
 
def upwind_advect(q, c):
    return q - c * (q - np.roll(q, 1))  # u > 0
 
def lax_wendroff_advect(q, c):
    flux_diff = 0.5 * c * (np.roll(q, -1) - np.roll(q, 1))
    visc = 0.5 * c * c * (np.roll(q, -1) - 2 * q + np.roll(q, 1))
    return q - flux_diff + visc
 
q_ftcs, q_up, q_lw = q0.copy(), q0.copy(), q0.copy()
for _ in range(n_steps):
    q_ftcs = ftcs_advect(q_ftcs, cfl)
    q_up = upwind_advect(q_up, cfl)
    q_lw = lax_wendroff_advect(q_lw, cfl)
 
# 厳密解: 1周して元の位置
q_exact = q0
print("FTCS  L1 err:", np.abs(q_ftcs - q_exact).sum() * dx)
print("Upwind L1 err:", np.abs(q_up - q_exact).sum() * dx)
print("LW    L1 err:", np.abs(q_lw - q_exact).sum() * dx)

cfl=0.8、400ステップ(=1周期)で厳密解は初期形状そのままです。典型的な結果は次の通り。

FTCS: 振幅が約 $1.78^{400}$ 倍に成長してNaNか inf に発散。
Upwind: 矩形パルスがガウス様に滲み、L1誤差が0.04付近に収まる。
Lax–Wendroff: 振幅は保たれるが、角の後ろに小波(Gibbs)が残る。

ひとつの小さなコードがこの記事のすべての主張 — 暴走、滲み、振動 — をそのまま再現します。

自分で触ってみる#

下のシミュレーションで直接操作してみましょう。CFLスライダーを0.05から1.5まで動かすと、3つのスキームが同じ格子の上でまったく違う運命をたどります。

CFL number c = u·dt / dx: 0.80

centered (FTCS)1st-order upwindLax–Wendroffanalytic shifttry CFL > 1 to break the safe zone

CFLをちょうど1.0にすると、風上(シアン)の線が解析解(点線)とほぼ完全に重なります — 1次風上は $c=1$ で厳密です。これは人工粘性係数 $u\Delta x(1-c)/2$ がそこで消えるからです。CFLを0.5まで下げると風上の滲みが最大になり、FTCSは依然発散。CFLを1.05に上げただけで風上もLax–Wendroffも発散します — ステンシルが特性線を追いつけないのです。

1次では足りないが、2次にはまた別の頭痛#

風上の人工粘性は衝撃波の近くでは歓迎されます — 非物理的振動を平らにします。しかし滑らかな領域でも同じように鈍るので、精度は2階微分レベルで止まります。実用CFDコードはたいてい次の2つのうちどちらかを取ります。

高次スキーム(Lax–Wendroff, MUSCL, WENO)を使い、衝撃近くで自動的に次数を落とすフラックスリミタやスロープリミタを組み込む。
因果構造を厳密に解くRiemann解法(Godunov, HLLC, Roe)をセル境界ごとに呼ぶ。

どちらも結局「どこで人工粘性を受け入れ、どこで拒むか」を決める仕組みです。1次風上はどこでも受け入れます — 最も単純で、最も安全で、最も不正確です。

覚えておく3行#

中心差分FTCSはCFLによらず発散する。 暴走の原因はステンシルではなく因果性違反です。
1次風上は安定の代償として人工粘性 $u\Delta x/2$ を払う。 $0\le c \le 1$ で安定、 $c=1$ で厳密。
CFL条件はすべての陽的スキームの必要条件であって十分条件ではない。 安定性は因果性とステンシル到達の両方が揃って初めて得られます。