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

N=200、c=0.8なら一周は250ステップです。コードの400ステップは1.6周で、初期パルスとは位置がずれます。

1960年代初頭、von Neumannの近くにいた人たちが1次元線形移流へ中心差分を入れました。矩形パルスは静かに流れるはずでした。後ろの縁で振動が膨らみ、数ステップで格子の外へ出ました。

原因は次数ではありません。下流の値を同じ重みで引くと因果が壊れます。一格子の風上差分はその暴走を鎮める代わりに角を殺します。その境界がCFLです。

CFLの目盛りを0.05から1.5まで動かすと、同じ矩形パルスでFTCSは震え続け、風上は1の近くでだけ破線に重なります。

CFL number c = u·dt / dx: 0.80

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

目盛りを0.5に置くと風上がいちばん広がります。1.05を超えると風上もLax–Wendroffも一緒に飛びます。ステンシルが特性線より短くなった場所です。

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

1次元の線形移流は一行です。

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

$q$ は流れる量、 $u>0$ は一定の速度です。解は初期形状が速度 $u$ で平行移動したものです。

時間に前進オイラー、空間に中心差分(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$ はクーラン数です。1ステップで波が横切る格子数です。空間2次、時間1次なのでLTE(局所打切り誤差)は $\mathcal{O}(\Delta t,\Delta x^2)$ です。 $c$ をどれだけ小さくしても、このスキームは暴走します。

$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を超えるモードが一つあれば場は壊れます。

コードの $c=0.8$ で $\sin=1$ のモードの大きさは $\sqrt{1+c^2}=\sqrt{1.64}\approx 1.28$ です。1ステップに1.28倍、400ステップはその累乗です。下に書いた目分量の $1.78^{400}$ はこのモードの率ではありません。式が返すのは1.28です。

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

$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

増幅率は

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$ です。

減衰項の係数 $c(1-c)$ は $c=0.5$ で0.25と最大です。 $c=0.8$ では0.16です。最短波長では $|g|^2 = 1 - 4\cdot 0.16 = 0.36$ 、つまり $|g|=0.6$ です。 $c=0.5$ ではそのモードの $|g|$ が0まで落ちます。 $c=1$ では引く項が消えて、すべてのモードが $|g|=1$ です。鈍りの頂点は0.5であり、安定域の端ではありません。

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

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}

前の項は中心差分です。後ろの項は2次微分の離散形です。だから風上が実際に積分している式は

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

です。移流を解こうとして粘性 $\nu_{\text{num}} = u\Delta x/2$ を差し込んだ修正方程式です。人工粘性(artificial viscosity)、あるいは numerical diffusion です。

この粘性が暴走を殺します。代償として矩形パルスの角が厚さ $\sqrt{4\nu_{\text{num}} t}$ に広がります。 $\Delta x$ を小さくすれば粘性も減りますが、同じ時刻までステップが増えて、ぼやけは残ります。

あとで出てくる係数 $u\Delta x(1-c)/2$ がこの式と違う理由は $c$ です。上の $u\Delta x/2$ は $c\to 0$ 側です。 $c=0.8$ なら $(1-c)=0.2$ なので、実際に残る係数はその5分の1、つまり $u\Delta x/10$ です。 $c=1$ で0になるのは $(1-c)$ が付いた側だけです。

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

時刻 $t_{n+1}$ の $q_i$ は、時刻 $t_n$ に位置 $x_i - u\Delta t$ にあった値です。その足が $x_{i-1}$ と $x_i$ の間にあれば、二つの格子値で補間できます。 $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です。

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ステップについて本文が残している目分量は次です。

FTCS: 振幅が約 $1.78^{400}$ 倍に増えて NaN または inf になります。
Upwind: 矩形パルスがガウス型に広がり、L1誤差は0.04付近に留まります。
Lax–Wendroff: 振幅は保たれますが、角の後ろにさざ波(Gibbs)が残ります。

c=0.5が風上をいちばん広げる#

$c(1-c)$ が0.5で最大だという事実が、画面のぼやけと一致します。 $c=1$ では $(1-c)$ が0なので、1次風上はその矩形を一格子ずつそのまま動かします。 $c=0.8$ はその間の0.16です。安定だからといって鈍らないわけではありません。0.8は0.5より粘性が小さいだけです。

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

風上の粘性は衝撃の近くの偽の振動を押さえてくれます。滑らかな領域でも同じ項がオンのままなので、誤差は2次微分の大きさに縛られます。実務のコードは次の二つへ進みます。

より高次(Lax–Wendroff、MUSCL、WENO)を使い、衝撃の近くでだけ次数を下げるフラックスリミッタまたはスロープリミッタを付けます。段差でLax–Wendroffが震える場所をリミッタで押さえる計算が、その選択です。
セル面でRiemannソルバー(Godunov、HLLC、Roe)を呼び、因果を直接積分します。

どちらも「粘性をどこで受け、どこで断るか」を決める装置です。1次風上はすべてのセルで受けます。単純で、安全で、いちばん広がります。

400ステップを一周と書くと失敗する#

$N=200$ 、 $u=1$ 、 $c=0.8$ なら $\Delta t = c/N = 0.004$ です。長さ1の領域を一周するステップ数は $1/\Delta t = 250$ 、つまり $N/c$ です。400ステップの移動距離は $400\times 0.004 = 1.6$ です。幅0.2のパルスが0.6だけずれると $[0.2,0.4]$ と $[0.8,1.0]$ は重なりません。鋭い二つの矩形のL1は0.4です。

本文の風上L1 0.04は、その比較にはなれません。0.04を一周分のぼやけとして読むなら、ステップは250である必要があります。1次元の風上が $c=1$ まで生きるこの上限は、運ぶものが界面になったとき、なぜずっと小さいCFLで止まるのでしょうか。その界面移流で0.05が残った場所を開くと、ステンシルの長さと界面の拘束がどこで分かれるかが見えます。