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

1952 年，Cornelius Lanczos、Magnus Hestenes 和 Eduard Stiefel 在一次学术会议上，偶然且独立地提出了同样的结果。那是一个仅需 $N$ 次迭代即可求解对称正定 (SPD) 矩阵的算法。这就是我们今天所称的共轭梯度法 (Conjugate Gradient, CG)。有趣的是，70 年后的今天，CFD 求解器的核心依然鲜活地保留着这个想法——在 Krylov 子空间中寻找解。本文将通过 Python 和交互式模拟，对比 CFD 中最常用的四种求解器（Jacobi、Gauss-Seidel、CG 和 BiCGSTAB）的收敛速度。读完本文，你将能直接判断在实际工作中“针对该问题应选择哪种求解器”。

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

对纳维-斯托克斯 (Navier–Stokes) 方程进行离散化后，每一步都会产生 $A\mathbf{x}=\mathbf{b}$ 形式的大型线性系统。其中 $A$ 是系数矩阵，$\mathbf{x}$ 是单元中心变量矢量，$\mathbf{b}$ 是右侧源项。对于 100 万个单元的 3D 网格，$A$ 将成为 $10^6 \times 10^6$ 大小的稀疏矩阵。在这种规模下，LU 分解（将矩阵分解为下三角和上三角乘积的直接法）会同时耗尽内存和 CPU。这就是 OpenFOAM、SU2、Gerris 等开源 CFD 代码无一例外使用迭代法 (iterative method) 的原因。

迭代法的基本思想很简单：从初始猜测值 $\mathbf{x}_0$ 出发，逐渐修正 $\mathbf{x}$，直到残差 $\mathbf{r} = \mathbf{b} - A\mathbf{x}$ 降至容差以下。问题在于如何定义“逐渐”。正是这个定义区分了不同求解器的名称。

Krylov 子空间的概念#

初始残差 $\mathbf{r}_0$ 与矩阵 $A$ 反复相乘后得到的矢量所构成的空间如下：

$\mathcal{K}_m$ 即为 Krylov 子空间，$m$ 为维数，$\text{span}$ 为通过线性组合填充的空间。Krylov 系列求解器在 $\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m$ 中寻找解，并使余下的残差 $\mathbf{r}_m$ 满足某种最优性条件。

将解限定在子空间内有两个优点：第一，除了矩阵-矢量乘法外没有其他运算，从而完整保留了稀疏结构；第二，理论上可以在该子空间内提取出最优解。CG 在 Krylov 空间内精确地使 A-范数误差最小化，而 GMRES 则使 2-范数残差最小化。

四种求解器的特性对比#

对于压力泊松方程 $-\nabla^2 p = f$，由于其具有 SPD 特性，CG 是标准选择。而包含平流-扩散的动量方程是非对称的，因此通常选择 BiCGSTAB 或 GMRES。OpenFOAM 中的 PCG/PBiCGStab 以及 SU2 中的 BCGSTAB 正是这种分工。

Python 演示收敛速度#

以 1D 泊松问题为例：将 $-u'' = 1$ 且边界条件为 $u(0)=u(1)=0$ 进行有限差分离散，会得到一个三对角矩阵。让我们直接对比 Jacobi 和 CG。

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

当 Jacobi 运行 2000 次仍无法突破 $10^{-2}$ 时，CG 仅在 $N=60$ 次后就达到了机器精度。理论上，CG 在最多 $N$ 次迭代内即可达到精确解（忽略浮点误差）。面对 SPD，CG 简直像拥有超能力。

观察收敛竞赛#

在下面的模拟中亲自尝试一下。通过滑块更改网格大小 $N$ 和最大迭代次数，四种求解器的残差曲线将以对数刻度重新绘制。

当 $N$ 增大到 80 时，Jacobi 和 Gauss-Seidel 的曲线几乎变为水平。而 CG 则呈阶梯状下降，并在 $N$ 附近迅速收敛。BiCGSTAB 在波动中下降，这种“锯齿状”收敛是潜在停滞（break-down，分子或分母爆炸为 0 导致数值崩溃）风险的可视化信号。在实战中，如果 BiCGSTAB 突然发散，通常的做法是强化预条件子或切换到 GMRES。

实际应用选择指南#

压力泊松方程、仅含扩散的热传递：CG + ILU(0) 预条件子。如果能保证对称性，选择其他方案都是浪费。

非对称动量方程、平流-扩散混合：首选 BiCGSTAB。如果收敛波动大或停滞，切换到 GMRES($m$)。$m$ 是重启周期，决定了内存与速度的权衡。行业惯例通常在 $m=30$ 左右。

并行 CFD：由于其串行特性，Gauss-Seidel 在 GPU 上表现不佳。过去十年的标准做法是基于 Jacobi 的平滑器 (smoother) 与代数多网格 (AMG) 的组合。

SU2 中的 LU-SGS：除 Krylov 外，SU2 还搭载了 LU-SGS (Lower-Upper Symmetric Gauss-Seidel) 这种固定迭代法。在压缩流动中结合隐式时间推进时，它的内存占用要小得多，但其稳健性不如 Krylov。

提示：在实际应用中，“预条件子的选择”对收敛速度的影响远大于“求解器的选择”。PCG (Preconditioned CG) + AMG 的组合往往比单纯使用 CG 快 100 倍。

核心总结#

• Jacobi 和 Gauss-Seidel 虽然简单，但在实际 CFD 问题中收敛率仅为 $1-O(1/N^2)$，表现堪忧。问题规模越大，迭代次数呈平方级增长。
• Krylov 系列（CG、BiCGSTAB、GMRES）仅需矩阵-矢量乘法即可在 $N$ 次内达到机器精度。SPD 选 CG，非对称选 BiCGSTAB 或 GMRES 是实际应用中的默认选择。
• 预条件子比求解器更重要。ILU(0)、AMG、Jacobi 平滑器——这三者的组合决定了 10 倍到 100 倍的速度差异。

下方调节 κ，立刻看到各求解器的极限。