隐式扩散与托马斯算法 — 用一次三对角求解换来无条件稳定

显式扩散的时间步被网格间距的平方锁死。在一维里拧开这把锁，代价只是一次三对角扫描。

晚上11点留下的计算，早上打开全是NaN。网格刚被切成一半，时间步没动。漏掉CFL(Courant–Friedrichs–Lewy)那一行，代价是集群一小时。

网格间距减半，显式扩散的时间步就变成四分之一。1024个点，也要为了不炸掉走上几十万步。出路很短。把右端改成未来时刻，再用托马斯法解那条三对角。

把r滑条从0.05拖到5，同一条高斯脉冲在显式和隐式上从哪里裂开，会叠在同一张图里。

CFL number r = D·dt / dx²: 0.45Explicit: stable (r ≤ 0.5)Implicit: unconditionally stable

0.5以下两条曲线几乎重合。一过0.5，橙色（显式）开始锯齿抖动；过了1就飞出画面，青色（隐式）的线还在。

CFL 锁住了细网格#

扩散方程如下。

\frac{\partial q}{\partial t} = D \frac{\partial^2 q}{\partial x^2}

其中 $q$ 是浓度（或温度）， $D$ 是扩散系数。最简单的显式格式（forward Euler + 中心差分）写成

q_i^{n+1} = q_i^{n} + r\,(q_{i+1}^{n} - 2 q_i^{n} + q_{i-1}^{n}),\quad r = \frac{D\,\Delta t}{\Delta x^2}

要稳定就需要 $r \le 1/2$ 。一行 von Neumann 分析：最短波长（波长 $= 2\Delta x$ ）的放大因子是 $|1 - 4r|$ ，要让它 $\le 1$ 就得 $r \le 1/2$ 。

从 $1024 \to 2048$ ， $\Delta t$ 变成 $1/4$ 。再切一刀到 $4096$ ，相对1024点的 $\Delta t$ 是 $1/16$ ，步数是16倍。点数加到8倍，走到同一时刻要64倍的步。账单来自这个平方，不来自单元个数。

$r = 0.6$ 时最短波的因子是 $|1 - 4r| = 1.4$ 。60步就是这个因子的60次方在吃模态。打印出来的 $1.33\times 10^{12}$ 里，还乘着初始脉冲投在这个模态上的那一份。 $r = 0.4$ 时因子是 $0.6$ ，显式峰值停在 $4.21\times 10^{-1}$ 。

把未来值搬到右边 — 隐式的简单把戏

改动只有一行。把右边的 $q^n$ 换成 $q^{n+1}$ 。

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

未来值在右边，就不能再一格一格往前算。所有单元同时绑在一起，变成矩阵方程。

A\,\mathbf{q}^{n+1} = \mathbf{q}^{n}

同样的分析给出放大因子 $1/(1 + 4r)$ ，正的 $r$ 下永远小于1。把 $\Delta t$ 加大也不会发散。这就是无条件稳定（unconditionally stable）。

代价是精度。过大的 $\Delta t$ 稳，但不准。即便如此，不炸掉这件事，对把作业放到过夜的人并不小。

同样的 $r = 0.4$ 、同样的60步，峰值是显式 $4.21\times 10^{-1}$ 、隐式 $0.4456$ 。都在稳定区，第三位已经分开。物理时刻相同，差的是截断误差。

三对角矩阵与托马斯的一次扫描

写出矩阵 $A$ 。一维里单元 $i$ 只和 $i-1$ 、 $i+1$ 耦合。

A = \begin{pmatrix} 1+r & -r & & & \\ -r & 1+2r & -r & & \\ & -r & 1+2r & -r & \\ & & \ddots & \ddots & \ddots \\ & & & -r & 1+r \end{pmatrix}

这是 三对角（tridiagonal） 矩阵。只有三条对角线非零。一般的高斯消元是 $\mathcal{O}(N^3)$ ，三对角则缩到 $\mathcal{O}(N)$ 。这个算法叫托马斯算法。1949年 Llewellyn Thomas 写在IBM内部备忘录里。

托马斯算法扫两遍。

前向消元 — 从上往下，一边消去次对角，一边算出新的对角 $b'$ 和右端 $d'$ 。
后向回代 — 从最后一行倒着代回去。

每遍每个单元4到5次浮点运算，存储只要三个数组。完整矩阵不必留在内存里。

只有边界对角是 $1+r$ ，内部是 $1+2r$ 。零通量把幽灵单元设成边界值，外侧邻居就消失了。对角上的 $2r$ 掉下一个 $r$ ，指向那个邻居的 $-r$ 也一起走。所以代码里是 b[0] = 1 + r。把内部那一行抄到边界上，总量会漏。

用 NumPy 比较显式与隐式#

import numpy as np
 
def explicit_diffusion_step(q, r):
    """前向欧拉 + 中心差分。"""
    next_q = q.copy()
    next_q[1:-1] = q[1:-1] + r * (q[2:] - 2 * q[1:-1] + q[:-2])
    return next_q
 
def thomas_solver(a, b, c, d):
    """以 O(N) 解三对角系统 A x = d。
    a: 次对角(a[0] 不使用)
    b: 主对角
    c: 上对角(c[-1] 不使用)
    d: 右端项
    """
    n = len(d)
    cp = np.zeros(n)
    dp = np.zeros(n)
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    for i in range(1, n):
        m = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / m
        dp[i] = (d[i] - a[i] * dp[i - 1]) / m
    x = np.zeros(n)
    x[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x
 
def implicit_diffusion_step(q, r):
    """后向欧拉 — 用托马斯法求解 (I - r·L) q^{n+1} = q^n。"""
    n = len(q)
    a = np.full(n, -r)
    b = np.full(n, 1 + 2 * r)
    c = np.full(n, -r)
    b[0] = 1 + r        # 零通量边界
    b[-1] = 1 + r
    return thomas_solver(a, b, c, q.copy())
 
# 用高斯脉冲对比
N = 80
x = np.linspace(0, 1, N)
q0 = np.exp(-((x - 0.5) ** 2) / 0.005)
D, dx = 1.0, 1.0 / N
r_values = [0.4, 0.6, 2.0]
for r in r_values:
    dt = r * dx**2 / D
    qe, qi = q0.copy(), q0.copy()
    for _ in range(60):
        qe = explicit_diffusion_step(qe, r)
        qi = implicit_diffusion_step(qi, r)
    print(f"r={r:.2f} | explicit max={np.max(np.abs(qe)):.2e} | implicit max={np.max(np.abs(qi)):.4f}")
# r=0.40 | explicit max=4.21e-01 | implicit max=0.4456
# r=0.60 | explicit max=1.33e+12 | implicit max=0.3873
# r=2.00 | explicit max=2.78e+38 | implicit max=0.2517

$r=0.6$ 时显式方法60步就被推到 $10^{12}$ 。同样的时间步下，隐式只是平稳地摊开，峰值是 $0.3873$ 。 $r=2$ 时显式到达 $2.78\times 10^{38}$ ，隐式峰值掉到 $0.2517$ 。

打印的三行并不是同一时刻

$\Delta t = r\,\Delta x^2 / D$ ，且 $N = 80$ 、 $D = 1$ 时，60步的物理时间正比于 $r$ 。 $r = 2$ 那一行比 $r = 0.4$ 远五倍。把峰值 $0.4456$ 和 $0.2517$ 并排放，说隐式“死得更多”，是假比较。要对齐时刻，得让 $r$ 乘步数相同。 $r = 2$ 的12步和 $r = 0.4$ 的60步是同一时刻，打印结果并没有做这个对齐。

隐式不是万能解的三个原因

无条件稳定不等于白给。代价分成三路。

1. 精度 vs 稳定性 — 稳定只保证振幅不炸，不保证答案对。后向欧拉的时间精度是 $\mathcal{O}(\Delta t)$ 。大 $\Delta t$ 还要精度时，换成 Crank–Nicolson（ $\mathcal{O}(\Delta t^2)$ ，仍然无条件稳定）更合适。

Crank–Nicolson 最短波的因子是 $(1-2r)/(1+2r)$ 。 $r = 2$ 时是 $-0.6$ 。绝对值小于1，所以不炸，但符号为负，每一步相位都翻过来。 $r$ 再加大，这个值会贴到 $-1$ 。

2. 非线性 — 扩散系数依赖解（ $D = D(q)$ ）时，矩阵每步都变，托马斯上面还要一层非线性求解（Picard 或 Newton）。

3. 多维 — 2D/3D 里矩阵不再是三对角。五点、七点模板会在远离主对角的地方生出旁带。这时用 ADI（交替方向隐式）拆成两到三次一维三对角，或者改走从 Jacobi 到 BiCGSTAB 的迭代法。

第三点会在从可压缩走到不可压缩时原样再出现。压力泊松方程 $\nabla^2 P = -\rho \nabla\!\cdot(\mathbf{u}\!\cdot\!\nabla\mathbf{u})$ 是同样结构的大稀疏矩阵。解过隐式扩散的人，写压力步用的也是同一套工具。

把大r写成正确答案就会失败#

$r = 2$ 的隐式峰值 $0.2517$ 不是“更接近真解”的记录。物理时间长了五倍，时间误差正比于 $\Delta t$ 。留下曲线之前如果不写上这些步的物理时间，比较就是假的。

网格一到二维，这段三对角代码就不能原样再用。当不组矩阵、只做扫描看起来更省时，显式炸掉之后不组矩阵就做隐式推进的 LU-SGS 能把这份稳定保证带到多远？扫描方向一旦错了，这份保证会不会一起碎掉？这一问还开着。