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

1960年代初，von Neumann 旁边的人把中心差分塞进一维线性对流。方波本该安静地漂。后缘的振荡长起来，几步就冲出了网格。

原因不是阶数。用同样的权重去取下游值，因果就断了。一格迎风压住失控，代价是把棱角磨掉。那条界就是 CFL。

把 CFL 刻度从 0.05 拨到 1.5，同一条方波上 FTCS 一直在抖，迎风只在 1 附近才和虚线重合。

刻度放在 0.5 时，迎风铺得最开。一过 1.05，迎风和 Lax–Wendroff 一起飞出去。模板比特征线短了。

最自然的差分最先崩溃

一维线性对流只有一行。

$q$ 是被输运的量，$u>0$ 是常速度。解就是初始形状以速度 $u$ 平移。

时间用前向欧拉，空间用二阶中心差分，更新写成

$c \equiv u\,\Delta t / \Delta x$ 是库朗数，一步之内波跨过的格子数。空间二阶、时间一阶，所以 LTE（局部截断误差）是 $\mathcal{O}(\Delta t,\Delta x^2)$。$c$ 再小，这个格式还是会失控。

代入 $q^n_i = g^n e^{i k x_i}$，放大因子是

于是 $|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.28 倍，400 步是这个因子的 400 次方。下面那句粗估 $1.78^{400}$ 不是这个模态。公式给出的是 1.28。

信息只从一侧流来

$u>0$ 时，时刻 $t_{n+1}$ 的 $q_i$ 只该看上游，也就是 $x<x_i$。中心差分把下游值 $q_{i+1}^n$ 以同样的权重拿过来。它在读流动还没到的那一侧。

迎风把方向固定住。

放大因子是

且 $|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，不在稳定区间的边上。

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

把一阶迎风的斜率拆开，得到恒等式。

前一项是中心差分。后一项是二阶导数的离散。所以迎风实际在积分的是

想解的是对流，算法塞进了粘性 $\nu_{\text{num}} = u\Delta x/2$。这就是人工粘性（artificial viscosity），也叫数值扩散。

这份粘性杀死失控。代价是方波的棱角摊成厚度 $\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$，真正留下的系数是它的五分之一，即 $u\Delta x/10$。只有带 $(1-c)$ 的那一个，才会在 $c=1$ 变成 0。

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 条件就是这个几何。

数值依赖域（模板）必须包住物理依赖域（特征线）。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）。

c=0.5把迎风铺得最开#

$c(1-c)$ 在 0.5 最大，和画面上的抹平对得上。$c=1$ 时 $(1-c)$ 为 0，一阶迎风把方波原样逐格平移。$c=0.8$ 夹在中间，是 0.16。稳定不等于没那么钝。0.8 只是比 0.5 的粘性小。

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

迎风粘性在激波附近是有用的。它把非物理的抖动压平。光滑区里同一项也开着，误差就被钉在二阶导数的量级。实际代码走两条路之一。

• 用更高阶（Lax–Wendroff、MUSCL、WENO），再加一个只在激波附近降阶的通量限制器或斜率限制器。用限制器压住 Lax–Wendroff 在台阶上抖动的计算就是这个选择。
• 在单元面上调用 Riemann 求解器（Godunov、HLLC、Roe），直接积分因果。

两条都是在决定粘性哪里收下、哪里拒绝。一阶迎风每个格子都收。简单，安全，铺得最开。

把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。一维迎风能活到 $c=1$。被输运的东西一旦变成界面，上限为什么会掉得那么低？打开界面输运里剩下 CFL 0.05 的地方，就能看到模板长度和界面约束在哪里分开。