选择 CFD 线性求解器：从 Jacobi 到 BiCGSTAB 的收敛速度对比

矩阵保持 SPD 时，CG 会在 N 次内做完。对流一旦破坏对称，这个位置就不再属于 CG。

1952 年的会场上，Lanczos、Hestenes 和 Stiefel 各自带着一份稿子碰了面。那是把对称正定（SPD）矩阵在 $N$ 次迭代里做完的算法。名字后来才叫做共轭梯度法（CG）。七十年过去，压力方程仍在那个子空间里找解。

把网格 N 从 20 推到 80，四条残差曲线会整幅重画。把最大迭代上限抬高，谁先贴到地板的次序会变。

N 在 40 附近时，CG 折向下，Jacobi 几乎躺在对数轴上。把 N 推到 80，间距拉得更大。BiCGSTAB 的锯齿表示残差不会光滑地死去。

为什么 CFD 采用迭代法求解线性系统#

把 Navier–Stokes 离散以后，每一步都会掉下一个 $A\mathbf{x}=\mathbf{b}$ 。 $A$ 是系数矩阵， $\mathbf{x}$ 是单元中心未知量， $\mathbf{b}$ 是右端项。一百万个单元时， $A$ 是 $10^6 \times 10^6$ 的稀疏矩阵。LU 分解（拆成下三角与上三角之积的直接法）在这个尺寸上同时烧掉内存和时间。OpenFOAM、SU2、Gerris 走迭代法，原因在这里。

从初始猜测 $\mathbf{x}_0$ 出发，一直改 $\mathbf{x}$ ，直到残差 $\mathbf{r} = \mathbf{b} - A\mathbf{x}$ 落到容差以下。“一点一点”怎么定义，就是求解器的名字。

Krylov 子空间的概念#

用初始残差 $\mathbf{r}_0$ 反复去乘 $A$ ，这些向量张成一个空间。

\mathcal{K}_m(A, \mathbf{r}_0) = \text{span}\{\mathbf{r}_0,\, A\mathbf{r}_0,\, A^2 \mathbf{r}_0,\, \dots,\, A^{m-1} \mathbf{r}_0\}

$\mathcal{K}_m$ 是 Krylov 子空间， $m$ 是维数， $\text{span}$ 是线性组合填满的集合。更新限制在 $\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m$ 里，剩下的残差 $\mathbf{r}_m$ 要满足某个最优条件。

运算到矩阵-向量乘法就结束。稀疏图案不会被拆开。CG 在这个空间里最小化 A-范数误差，GMRES 最小化 2-范数残差。后者发生在 Arnoldi 把向量接起来的子空间里面。

四种求解器的特性对比

求解器	矩阵要求	内存占用	优点	缺点
Jacobi	推荐对角占优	极低	易于并行化，10 行代码实现	慢（谱半径 $1-O(1/N^2)$ ）
Gauss-Seidel	任意	极低	比 Jacobi 快约 2 倍	本质上是串行的，不利于 GPU
CG	必须为 SPD	4 个矢量	对 SPD 最优，三项递推，轻量级	不支持非对称矩阵
BiCGSTAB	任意	7 个矢量	支持非对称，无需转置	收敛过程有波动，可能出现停滞 (break-down)

压力泊松方程 $-\nabla^2 p = f$ 是 SPD，CG 对得上。掺了对流的动量方程不对称，就落到 BiCGSTAB 或 GMRES。OpenFOAM 的 PCG/PBiCGStab 和 SU2 的 BCGSTAB 就是这个分工。

Python 演示收敛速度#

一维泊松 $-u'' = 1$ ，边界 $u(0)=u(1)=0$ ，有限差分之后是三对角阵。在 $N=60$ 上，Jacobi 和 CG 用同一个容差 $10^{-10}$ 。

import numpy as np
 
def build_poisson_1d(n):
    A = np.zeros((n, n))
    for i in range(n):
        A[i, i] = 2.0
        if i > 0:     A[i, i - 1] = -1.0
        if i < n - 1: A[i, i + 1] = -1.0
    return A, np.ones(n)
 
def jacobi_sweep(A, b, tol=1e-10, max_iter=2000):
    D  = np.diag(A)
    LU = A - np.diag(D)
    x  = np.zeros_like(b)
    b0 = np.linalg.norm(b)
    history = []
    for _ in range(max_iter):
        x = (b - LU @ x) / D
        res = np.linalg.norm(b - A @ x) / b0
        history.append(res)
        if res < tol:
            break
    return x, history
 
def cg_sweep(A, b, tol=1e-10, max_iter=2000):
    x  = np.zeros_like(b)
    r  = b - A @ x
    p  = r.copy()
    rr = r @ r
    b0 = np.linalg.norm(b)
    history = [np.sqrt(rr) / b0]
    for _ in range(max_iter):
        Ap    = A @ p
        alpha = rr / (p @ Ap)
        x    += alpha * p
        r    -= alpha * Ap
        rr_new = r @ r
        history.append(np.sqrt(rr_new) / b0)
        if history[-1] < tol:
            break
        p  = r + (rr_new / rr) * p
        rr = rr_new
    return x, history
 
A, b = build_poisson_1d(n=60)
_, hist_j  = jacobi_sweep(A, b)
_, hist_cg = cg_sweep(A, b)
print(f"Jacobi: {len(hist_j):4d} 次, 最终残差 {hist_j[-1]:.2e}")
print(f"CG:     {len(hist_cg):4d} 次, 最终残差 {hist_cg[-1]:.2e}")

输出是这样印出来的。

Jacobi: 2000 次, 最终残差 2.11e-02
CG:       60 次, 最终残差 1.87e-13

把对数条件数滑条往上推，同样的迭代预算里四条曲线停在不同的地方。迭代上限决定曲线在哪里被切断。

condition number κ = 10^3.0iters80

조건수 κ를 키우면 (그림에서 오른쪽으로) Jacobi는 거의 멈추고 CG/GMRES는 √κ 한계를 따른다. BiCGSTAB는 비대칭에 강하지만 잔차가 출렁인다.

默认 $\kappa=10^3$ 、大约 80 次迭代时，Jacobi（红）几乎水平，CG（青）更陡。把 $\kappa$ 推向 $10^6$ ，Jacobi 即使给到 200 次也看不见地板。上面 $N=60$ 的输出是这张图上的一个点。Jacobi 用 2000 次停在 2.11e-2，CG 用 60 次到 1.87e-13。

把打印出来的残差折回每一步，数量级就出来了。初始残差取 1，2000 次之后的 2.11e-2 相当于平均衰减 0.998，也就是 0.0211 的 2000 次方根。这和表里的谱半径 $1-O(1/N^2)$ 是同一个量级。 $N=60$ 时 $1/N^2$ 是 1/3600，大约 0.00028。0.998 读成这个量级前面再乘上几倍常数，就对上了。

这个衰减若保持下去，降到 $10^{-10}$ 大约要 1.2 万次。已经到了 2.11e-2 之后，仍还要大约 1 万次。 $N$ 加倍， $1/N^2$ 变成四分之一，同样残差所需的次数也变成四倍。格子还要加密的计算里，若把 Jacobi 留作正式迭代，预算就跟着单元数的平方走。

CG 的 60 次和 $N$ 严格相等。它把“最多 $N$ 次”的上限用完了，不是大约。1.87e-13 比 2.11e-2 小约 $1.1\times 10^{11}$ 倍。双精度里，能加到 1 附近的最小数在 $10^{-16}$ 上下，60 次正交化磨掉了两三位数。次数本身停在 $N$ 。

一百万个单元乘 8 字节，一条向量是 8 MB。CG 的 4 条向量是 32 MB，BiCGSTAB 的 7 条是 56 MB。放在稀疏矩阵旁边，存残差是便宜的一侧。贵的是那一万次的尾巴。

关掉预条件之后的曲线

前面的曲线是预条件关掉的裸迭代。预条件矩阵一旦进去，横轴上一格的含义就变了。BiCGSTAB 垂直跳起的位置，是分母里的内积贴到了 0，也就是 break-down。只加迭代次数回不来。要么换预条件，要么改走 GMRES。重启周期 $m$ 的习惯值 30，是在内存和这次跳跃之间取的折中。

实际应用选择指南

压力泊松，以及只有扩散的传热，用 CG 加 ILU(0)。对称性还在的时候另选求解器，只是多存几条向量。

不对称的动量方程以 BiCGSTAB 为底。残差晃动就换成 GMRES( $m$ )。 $m$ 是重启周期。

并行时 Gauss-Seidel 是串行的，到了 GPU 上就使不上力。近十年更接近标准的是 Jacobi 光滑子加上代数多重网格（AMG）。

SU2 的 LU-SGS（Lower-Upper Symmetric Gauss-Seidel）是 Krylov 旁边的定常迭代。接到可压缩流动的隐式时间推进上，内存脚印更小。它没有 Krylov 那么稳。

改预条件比改求解器名字更能拉开迭代次数。PCG 加上 AMG，常常比裸的 CG 快 100 倍。差的就是上面那一万次的尾巴切不切得掉。

把残差 2.11e-2 当成收敛交出去的失误#

2000 次之后仍把 2.11e-2 交给下一个时间步，这个压力误差会走进下一轮非线性迭代。SPD 还在，这个位置属于 CG。对流把矩阵扭成不对称，CG 就当场失去资格，交给 BiCGSTAB。先挑名字、后量条件数，这个次序就是失误。

还没回答的问题是，要不要再堆一万次迭代。在粗网格上扫误差的 V 循环能不能切断这条 $1-O(1/N^2)$ 的尾巴，还是会在细网格上再次调用 Jacobi。下一次选择看的是这一下。