FTCS 颤抖时,迎风变钝 — 一阶对流格式、人工粘性与 CFL 的等价交换

1960 年代初期,John von Neumann 的同事们把世上最简单的 PDE — 一维线性对流 — 接上了最自然的中心差分。结果是失控:本该静静漂移的方波在尾沿长出震荡,几步之内就冲出网格。本文跟踪这场失控从何而来,为什么仅一格的迎风差分能压住它却以画面变钝为代价,以及网格步长与时间步长之间那条看不见的界线(CFL 条件)究竟说了什么。最后用 50 行 Python 与一个交互演示,在同一张网格上让三种格式赛跑。

最自然的差分最先崩溃#

一维线性对流方程只有一行。

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

其中 $q$ 是被输运的量, $u>0$ 为常速度。解很简单:初始形状以速度 $u$ 平移。

时间用前向 Euler、空间用中心(二阶)差分,更新规则如下。

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 数(波在一个时间步穿过的格点数)。空间二阶、时间一阶,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$ 。所有 Fourier 模态每步都增长。只要有一个 $|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$ 。每个模态要么衰减要么不变。

迎风 = 中心差分 + 隐藏的粘性#

更直观地看为什么迎风稳定。把一阶迎风差分按恒等式拆开:

\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}

右边第一项是中心差分,第二项恰是二阶导数的离散形式。所以一阶迎风实际上在求解修正方程

\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$ )会减小粘性,但要走到同一最终时刻就得多走几步。两者刚好抵消,所以钝化永远不会消失。

CFL — 网格追不上特征线的瞬间#

迎风仅在 $0\le c \le 1$ 内稳定,有一个几何上的理由。 $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 行追三种格式#

同一格网、同一初始条件、同一 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)
 
# 解析解: 转一圈回到原位
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 涟漪。

一段简短代码就把本文所有主张 — 失控、钝化、震荡 — 一并复现出来。

自己上手#

下面的模拟可以亲手操作。把 CFL 滑块从 0.05 拖到 1.5,三种格式在同一张网格上呈现完全不同的命运。

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,迎风曲线(青色)几乎完全压在解析平移(虚线)上 — 一阶迎风在 $c=1$ 处精确。原因是人工粘性系数 $u\Delta x(1-c)/2$ 在那里恰好为零。把 CFL 降到 0.5,迎风的钝化最严重,而 FTCS 仍然发散。把 CFL 推到 1.05,迎风和 Lax–Wendroff 都发散 — 模板追不上特征线。

一阶不够,但二阶又有新麻烦#

迎风的人工粘性在激波附近受欢迎 — 它把非物理震荡抹平。但平滑区它同样钝化,精度卡在二阶导数级别。生产 CFD 代码通常二选一。

用更高阶格式(Lax–Wendroff, MUSCL, WENO),并在激波附近自动降阶的通量限制器或斜率限制器上叠加。
在每个单元面调用Riemann 求解器(Godunov, HLLC, Roe),精确求解因果结构。

两条路本质都在决定:在哪里接受人工粘性,在哪里拒绝。一阶迎风到处接受 — 最简单、最安全、最不准确。

三行要点#

中心差分 FTCS 与 CFL 无关地发散。 失控源自因果性,不是模板。
一阶迎风以人工粘性 $u\Delta x/2$ 换稳定。 仅在 $0\le c \le 1$ 稳定; $c=1$ 时精确。
CFL 是所有显式格式的必要条件,而非充分条件。 稳定性需要因果性与模板宽度同时到位。