第 10 章：代数方程组的求解（Solving the System of Algebraic Equations）

10.1 引言（Introduction）

离散化的产物是一个线性方程组 \(\mathbf{A}\boldsymbol{\phi} = \mathbf{b}\)，未知量 \(\boldsymbol{\phi}\) 位于网格单元的形心位置。矩阵 \(\mathbf{A}\) 中的系数是线性化过程与网格几何的综合结果，向量 \(\mathbf{b}\) 则汇集了源项、常数、边界条件以及不可线性化的非定常分量。线性方程组的求解技术大致分为直接法与迭代法两大类，两类之下又各有若干子类。由于流动问题高度非线性，其线性化产生的系数一般依赖于解本身；同时在每一次非线性迭代中并不需要获得准确解——基于这两个原因，直接法在 CFD 中极少被采用。迭代法更受青睐，因为每次迭代计算量低、内存占用小，与 CFD 应用的契合度更高。本章先介绍几种适用于结构网格与非结构网格的直接法（Gauss 消去、LU 分解、三对角与五对角矩阵算法），为后续讨论迭代法铺垫基础；随后审视若干基本迭代法的性能与局限，包括带预条件与不带预条件的 Jacobi、Gauss-Seidel、不完全 LU 分解以及共轭梯度法；最后介绍通常与迭代求解器联用以弥补其局限性的多重网格方法。

任意线性求解器的起点都是离散化产生的方程组，将其写成 \(\mathbf{A}\boldsymbol{\phi} = \mathbf{b}\)（式 10.1）的形式，其中 \(\mathbf{A}\) 是以 \(a_{ij}\) 为元素的系数矩阵，\(\boldsymbol{\phi}\) 是以 \(\phi_i\) 为分量的未知向量，\(\mathbf{b}\) 是以 \(b_i\) 为分量的右端源向量。按矩阵编号展开的形式为式 10.2。一般而言，上述矩阵的每一行代表控制体上的一个方程；非零系数出现在与该单元相邻的那些列上。系数 \(a_{ij}\) 度量的是 \(\phi_i\) 在控制体形心处的值与其邻居之间联系的强度。由于每个单元只与少数几个邻居相连（具体数量取决于离散域中单元的连接关系），许多系数为零，所得到的矩阵 \(\mathbf{A}\) 总是稀疏的——非零系数仅占矩阵元素的极小比例。如果再额外采用结构网格，\(\mathbf{A}\) 会呈带状分布，所有非零元素沿少数几条对角线对齐。因此，构造高效的线性求解方法应当利用这种结构。

如前所述，求解代数方程组的技术大致分为直接法与迭代法两类。直接法对矩阵 \(\mathbf{A}\) 求逆，一步得到 \(\boldsymbol{\phi} = \mathbf{A}^{-1}\mathbf{b}\)。当矩阵规模较大时，求逆在计算量与内存上都代价高昂，因此在 CFD 中基本上不可行——CFD 涉及的方程组是非线性的，系数依赖于解本身，必须使用迭代过程。迭代法求解器的算法则按需反复施加，直到达到预设的收敛判据，而不需要每次迭代都获得完全收敛的解。

本章先介绍几种适用于结构网格与非结构网格的直接线性求解器；接着讨论利用结构网格系数矩阵带状结构的求解算法；本章的重点则放在一类特定的迭代线性代数求解器上——经验证明它们与 FVM 配合时效率高、代价低，并且几乎所有基于有限体积法的代码都将其作为线性求解器来唯一实现。

10.2 直接法与高斯消去法（Direct or Gauss Elimination Method）

虽然直接法在求解稀疏线性方程组时代价高昂，但其讨论将为下一节介绍高效迭代法铺垫基础。最简单的求式 10.1 所描述线性方程组解的直接法是高斯消去法，下文将首先描述。高斯消去法将系统变换为等价的上三角系统的过程，催生了下三角—上三角（LU）分解法，也将一并讨论。在这种方法中，矩阵 \(\mathbf{A}\) 被分解为两个矩阵的乘积 \(\mathbf{L}\) 与 \(\mathbf{U}\)，其中 \(\mathbf{L}\) 为下三角矩阵、\(\mathbf{U}\) 为上三角矩阵。这一过程又称 LU 分解。此外，还将讨论利用 \(\mathbf{A}\) 的带状结构、适用于结构网格的直接法。

10.2.1 高斯消去（Gauss Elimination）

描述高斯消去法最直接的方式是从一个简单例子开始。为此，考虑如下两个未知量 \(\phi_1\) 与 \(\phi_2\) 的线性方程组：

\[ a_{11}\phi_1 + a_{12}\phi_2 = b_1 \quad (10.3) \]

\[ a_{21}\phi_1 + a_{22}\phi_2 = b_2 \quad (10.4) \]

该方程组可以通过从其中一个方程里消去某一个变量来求解（例如从式 10.4 中消去 \(\phi_1\)）。具体做法是将式 10.3 乘以 \(a_{21}/a_{11}\)，再从式 10.4 中减去所得方程，得到

\[ \left(a_{22} - \frac{a_{21}}{a_{11}}a_{12}\right)\phi_2 = b_2 - \frac{a_{21}}{a_{11}}b_1 \quad (10.5) \]

该方程只含一个未知量，因此可直接解出 \(\phi_2\)：

\[ \phi_2 = \frac{b_2 - \dfrac{a_{21}}{a_{11}}b_1}{a_{22} - \dfrac{a_{12}a_{21}}{a_{11}}} \quad (10.6) \]

知道 \(\phi_2\) 后，将其代回式 10.3 即可解得 \(\phi_1\)。执行此步可得

\[ \phi_1 = \frac{b_1}{a_{11}} - \frac{a_{12}}{a_{11}} \cdot \frac{b_2 - \dfrac{a_{21}}{a_{11}}b_1}{a_{22} - \dfrac{a_{12}a_{21}}{a_{11}}} \quad (10.7) \]

上述过程由两步组成。第一步对方程进行操作以消去某个未知量，最终结果是一个只含一个未知量的方程。第二步直接求解该方程并将结果代回其中一个方程，求出剩余未知量。同样的过程可以推广到式 10.1 或式 10.2 所描述的 \(N\) 元方程组，详见下文。

10.2.2 前向消去（Forward Elimination）

在以下推导中，\(\mathbf{A}\) 的第一行对应 \(\phi_1\) 的离散方程，第二行对应 \(\phi_2\) 的方程，以此类推，第 \(i\) 行对应 \(\phi_i\) 的方程。过程从消去 \(\mathbf{A}\) 第一行以下所有方程中的 \(\phi_1\) 开始。为从第 \(i\) 行（\(i = 2, 3, \ldots, N\)）消去 \(\phi_1\)，将第一行的系数乘以 \(a_{i1}/a_{11}\)，再从第 \(i\) 行中减去所得方程。完成此步后方程组变为式 10.8 的形式。

然后消去修改后的 \(\mathbf{A}\) 第二行以下所有方程中的 \(\phi_2\)。为从第 \(i\) 行（\(i = 3, 4, \ldots, N\)）消去 \(\phi_2\)，将第二行的系数乘以 \(a'_{i2}/a'_{22}\)，再从第 \(i\) 行中减去所得方程。接着从第三行以下的修改系数矩阵的方程中消去 \(\phi_3\)，如此继续，直至从第 \(N\) 行消去 \(\phi_{N-1}\)，得到一个系数矩阵变为上三角矩阵的等价方程组（式 10.9）。所得算法描述如下。

10.2.3 前向消去算法（Forward Elimination Algorithm）

For k = 1 to N-1
  {
    For i = k + 1 to N
      {
        Ratio = a_ik / a_kk
        {
          For j = k + 1 to N
            a_ij = a_ij - Ratio * a_kj
        }
        b_i = b_i - Ratio * b_k
      }
  }

10.2.4 回代（Backward Substitution）

式 10.9 给出的修改方程组表明，第 \(N\) 个方程中只有一个未知量 \(\phi_N\)。因此该方程可用来求解 \(\phi_N\)：

\[ \phi_N = \frac{b^{N-1}_N}{a^{N-1}_{NN}} \quad (10.10) \]

第 \(N - 1\) 个方程是 \(\phi_{N-1}\) 与 \(\phi_N\) 的函数。在求出 \(\phi_N\) 后，可由该方程得到 \(\phi_{N-1}\)：

\[ \phi_{N-1} = \frac{b^{N-2}_{N-1} - a^{N-2}_{N-1,N}\phi_N}{a^{N-2}_{N-1,N-1}} \quad (10.11) \]

过程继续反向进行；到达第 \(i\) 个方程时， \(\phi_{i+1}, \phi_{i+2}, \phi_{i+3}, \ldots, \phi_{N-1}, \phi_N\) 已经全部获得，因此 \(\phi_i\) 可用下式计算：

\[ \phi_i = \frac{b^{i-1}_i - \sum\limits_{j=i+1}^{N} a^{i-1}_{ij}\phi_j}{a^{i-1}_{ii}} \quad (10.12) \]

该过程一直进行到 \(\phi_1\) 被求出。用算法表示如下。

10.2.5 回代算法（Back Substitution Algorithm）

phi_N = b_N / a_NN
For i = N-1 to 1
  {
    Term = 0
    {
      For j = i + 1 to N
        Term = Term + a_ij * phi_j
    }
    phi_i = (b_i - Term) / a_ii
  }

通过主元置换（交换行以选取最大的主元素）避免除零以及减少大系统中舍入误差的技术是可行的，但此处不展开讨论。有兴趣的读者可参阅这方面的专著 [1–4]。所呈现的算法表明该方法代价昂贵——求解 \(N\) 元线性系统所需的操作数与 \(N^3/3\) 成正比，其中回代部分仅需 \(N^2/2\) 次算术运算。这种高昂的计算代价促使研究者为稀疏矩阵系统寻找更高效、更专门的求解器。

10.2.6 LU 分解（LU Decomposition）

另一种直接求解线性代数方程组的方法是 LU 或更一般的 PLU 分解（其中 \(\mathbf{P}\) 指上文提及的主元置换过程）；本书只涉及 LU 分解法，它是高斯消去法的变体。这些方法相对于高斯消去法的优势在于：一旦完成了 (P)LU 分解，对应不同右端向量 \(\mathbf{b}\) 的线性系统就可以反复求解，无需再做额外的消去运算——而高斯消去法仍需进行消去。

根据上一节中的消去过程，式 10.1 被变换为式 10.9 给出的上三角矩阵形式，可写为

\[ \mathbf{U}\boldsymbol{\phi} - \mathbf{c} = 0 \quad (10.14) \]

设 \(\mathbf{L}\) 为单位下三角矩阵（对角元素设为 1 以使分解唯一），其形式如式 10.15 所给出，使得式 10.14 左乘 \(\mathbf{L}\) 后恢复式 10.1。若此条件成立，则可写出

\[ \mathbf{L}(\mathbf{U}\boldsymbol{\phi} - \mathbf{c}) = \mathbf{L}\mathbf{U}\boldsymbol{\phi} - \mathbf{L}\mathbf{c} = \mathbf{A}\boldsymbol{\phi} - \mathbf{b} \quad (10.16) \]

由矩阵性质可得

\[ \mathbf{L}\mathbf{U} = \mathbf{A} \quad (10.17) \]

\[ \mathbf{L}\mathbf{c} = \mathbf{b} \quad (10.18) \]

式 10.17 表明 \(\mathbf{A}\) 被写成一个下三角矩阵与一个上三角矩阵的乘积，即 LU 分解。

10.2.7 分解步骤（The Decomposition Step）

下文给出求 \(\mathbf{L}\) 和 \(\mathbf{U}\) 系数的高效过程，称为 Crout 分解 [1–4]。原始的 Crout 算法使用单位上三角矩阵，而此处假设单位下三角矩阵。该过程基于将 \(\mathbf{L}\) 与 \(\mathbf{U}\) 相乘得到 \(\mathbf{A}\)（式 10.19）。

系数计算首先将 \(\mathbf{L}\) 的第一行与 \(\mathbf{U}\) 的所有列相乘，并与 \(\mathbf{A}\) 对应系数相等，得到

\[ u_{1j} = a_{1j}, \quad j = 1, 2, 3, \ldots, N \quad (10.20) \]

然后将 \(\mathbf{L}\) 的第二到第 \(N\) 行与 \(\mathbf{U}\) 的第一列相乘，得到

\[ \ell_{i1}u_{11} = a_{i1} \Rightarrow \ell_{i1} = \frac{a_{i1}}{u_{11}}, \quad i = 2, 3, \ldots, N \quad (10.21) \]

将该过程重复：将 \(\mathbf{L}\) 的第二行与 \(\mathbf{U}\) 的第二到第 \(N\) 列相乘，得到

\[ u_{2j} = a_{2j} - \ell_{21}u_{1j}, \quad j = 2, 3, \ldots, N \quad (10.22) \]

随后将 \(\mathbf{L}\) 的第三到第 \(N\) 行与 \(\mathbf{U}\) 的第二列相乘，得到

\[ \ell_{i2}u_{22} + \ell_{i1}u_{12} = a_{i2} \Rightarrow \ell_{i2} = \frac{a_{i2} - \ell_{i1}u_{12}}{u_{22}}, \quad i = 3, 4, \ldots, N \quad (10.23) \]

一般而言，将 \(\mathbf{L}\) 的第 \(i\) 行与 \(\mathbf{U}\) 的第 \(i\) 到第 \(N\) 列相乘，得到

\[ u_{ij} = a_{ij} - \sum_{k=1}^{i-1}\ell_{ik}u_{kj}, \quad j = i, i+1, \ldots, N \quad (10.24) \]

再将 \(\mathbf{L}\) 的第 \(i+1\) 到第 \(N\) 行与 \(\mathbf{U}\) 的第 \(i\) 列相乘，得到

\[ \ell_{ki} = \frac{a_{ki} - \sum_{j=1}^{i-1}\ell_{kj}u_{ji}}{u_{ii}}, \quad k = i+1, i+2, \ldots, N \quad (10.25) \]

对 \(\mathbf{L}\) 的第 \(N\) 行，其系数与 \(\mathbf{U}\) 的第 \(N\) 列相乘，由此可得 \(u_{NN}\)：

\[ u_{NN} = a_{NN} - \sum_{k=1}^{N-1}\ell_{Nk}u_{kN} \quad (10.26) \]

LU 分解过程算法化总结如下。

10.2.8 LU 分解算法（LU Decomposition Algorithm）

u_1j = a_1j,  j = 1 to N
ℓ_i1 = a_i1 / u_11,  i = 2 to N
For i = 2 to N-1
  {
    u_ij = a_ij - sum_{k=1}^{i-1} ℓ_ik * u_kj,  j = i, i+1, ..., N
    ℓ_ki = (a_ki - sum_{j=1}^{i-1} ℓ_kj * u_ji) / u_ii,  k = i+1, i+2, ..., N
  }
u_NN = a_NN - sum_{i=1}^{N-1} ℓ_Ni * u_iN

10.2.9 代入步骤（The Substitution Step）

将原始矩阵 \(\mathbf{A}\) 分解为 \(\mathbf{L}\) 与 \(\mathbf{U}\) 后，方程组可以通过式 10.18 与式 10.14 进行两步求解。需要注意的是，该两步过程等价于求解两个线性方程组，但因为 \(\mathbf{L}\) 与 \(\mathbf{U}\) 分别是下三角与上三角矩阵，所以求解得到简化。

第一步通过前向代入由式 10.18 得到向量 \(\mathbf{c}\)。该过程描述为

\[ c_1 = b_1 \]

\[ c_i = b_i - \sum_{j=1}^{i-1}\ell_{ij}c_j, \quad i = 2, 3, \ldots, N \quad (10.27) \]

第二步通过回代由式 10.14 求出 \(\boldsymbol{\phi}\) 值。过程描述为

\[ \phi_N = \frac{c_N}{u_{NN}} \]

\[ \phi_i = \frac{c_i - \sum_{j=i+1}^{N} u_{ij}\phi_j}{u_{ii}}, \quad i = N-1, N-2, \ldots, 3, 2, 1 \quad (10.28) \]

如果 \(\mathbf{A}\) 不再需要，\(\mathbf{L}\) 与 \(\mathbf{U}\) 的元素可以直接存储在原矩阵 \(\mathbf{A}\) 中。这是因为 \(\mathbf{A}\) 的元素只在计算对应的 \(\mathbf{L}\) 或 \(\mathbf{U}\) 元素时才需要。对 \(N \times N\) 方阵执行 LU 分解所需的操作数为 \(2N^3/3\)，是用高斯消去法求解同一系统所需操作数的两倍。再次强调，使用 LU 分解的优势在于同一矩阵 \(\mathbf{A}\) 对应多个不同 \(\mathbf{b}\) 向量的情形。然而，引入 LU 分解的主要原因是它构成了发展某些更高效的线性代数方程组迭代求解器的基础，下一节将介绍这些迭代求解器。

10.2.10 LU 分解与高斯消去（LU Decomposition and Gauss Elimination）

这一点也许并不直观，但高斯消去法可以用来执行 LU 分解。前文已表明，前向消去步骤的结果是上三角矩阵 \(\mathbf{U}\)。然而在这一过程中，\(\mathbf{L}\) 实际上已经被产生。\(\mathbf{L}\) 的元素就是消去过程中用以乘行的那些因子（在高斯消去算法中记为 ratio）。下面的算法假设 \(\mathbf{L}\) 为单位下三角矩阵，通过高斯消去法对 \(\mathbf{A}\) 进行 LU 分解。

10.2.11 基于高斯消去的 LU 分解算法（LU Decomposition Algorithm by Gauss Elimination）

u_1j = a_1j,  j = 1 to N
For k = 1 to N-1
  {
    For i = k + 1 to N
      {
        ℓ_ik = a_ik / a_kk
        {
          For j = k + 1 to N
            u_ij = a_ij - ℓ_ik * a_kj
        }
      }
  }

例 1：用 LU 分解法求解下列线性代数方程组（具体矩阵与右端向量按式 10.19 所给系数逐元素代入计算）。依据上述过程，元素按如下方式计算：\(u_{1j} = a_{1j}\) 得 \(u_{11} = 3, u_{12} = -1, u_{13} = 0, u_{14} = 0\)；\(\ell_{i1} = a_{i1}/u_{11}\) 得 \(\ell_{21} = -2/3, \ell_{31} = 0, \ell_{41} = 0\)；\(u_{2j} = a_{2j} - \ell_{21}u_{1j}\) 得 \(u_{22} = 16/3, u_{23} = -1, u_{24} = 0\)；\(\ell_{i2} = (a_{i2} - \ell_{i1}u_{12})/u_{22}\) 得 \(\ell_{32} = -3/8, \ell_{42} = 0\)； \(u_{33} = a_{33} - \ell_{31}u_{13} - \ell_{32}u_{23} = 45/8, u_{34} = -1\) ； \(\ell_{43} = (a_{43} - \ell_{41}u_{13} - \ell_{42}u_{23})/u_{33} = -16/45\) ； \(u_{44} = a_{44} - \ell_{41}u_{14} - \ell_{42}u_{24} - \ell_{43}u_{34} = 299/45\) 。

因此 \(\mathbf{L}\) 与 \(\mathbf{U}\) 矩阵如本节公式所示；\(\mathbf{c}\) 向量应满足 \(\mathbf{L}\mathbf{c} = \mathbf{b}\)，由前向代入得 \(c_1 = 3, c_2 = 6, c_3 = 29/4, c_4 = -19/45\)。原方程的解通过求解 \(\mathbf{U}\boldsymbol{\phi} = \mathbf{c}\) 得到：由回代过程可得 \(\phi_4 = -19/299, \phi_3 = 382/299, \phi_2 = 408/299, \phi_1 = 435/299\) 。

10.2.12 带状稀疏矩阵的直接法（Direct Methods for Banded Sparse Matrices）

高斯消去法与 LU 分解法适用于任何方程组。具体而言，它们可以用于求解本书涉及的结构或非结构网格上守恒方程离散化所产生的方程组。当采用结构网格方法时，离散化产生的方程组的系数矩阵非零元素沿少数几条对角线对齐。根据所使用的离散模板与所求解问题的维数，可能会出现三对角或五对角矩阵；针对这两种情形已发展出高效的算法，下文予以描述。

10.2.13 三对角矩阵算法（TDMA）

三对角矩阵算法（TDMA），又称 Thomas 算法 [5, 6]，用于求解具有三对角系数矩阵的代数方程组，写为

\[ a_i\phi_i + b_i\phi_{i+1} + c_i\phi_{i-1} = d_i, \quad i = 1, 2, 3, \ldots, N; \quad c_1 = b_N = 0 \quad (10.29) \]

对本书所采用的网格排列，\(i\) 指向图 10.1 中所示的网格点位置。当 \(i = 1\) 时，该方程可用来以 \(\phi_2\) 表示 \(\phi_1\)：

\[ \phi_1 = -\frac{b_1}{a_1}\phi_2 + \frac{d_1}{a_1} \quad (10.30) \]

类似地，对 \(i = 2\)，式 10.29 在式 10.30 的帮助下可将 \(\phi_2\) 仅表示为 \(\phi_3\) 的函数：

\[ \phi_2 = -\frac{a_1 b_2}{a_1 a_2 - c_2 b_1}\phi_3 + \frac{d_2 a_1 - c_2 d_1}{a_1 a_2 - c_2 b_1} \quad (10.31) \]

对 \(\phi_3\) 直到 \(\phi_N\) 重复同样的过程，可以推断一般情况下 \(\phi_i\) 可以表达为 \(\phi_{i+1}\) 的函数：

\[ \phi_i = P_i\phi_{i+1} + Q_i, \quad i = 1, 2, 3, \ldots, N \quad (10.32) \]

将 \(i - 1\) 时的式 10.32 与式 10.29 联立，得到

\[ \phi_{i-1} = P_{i-1}\phi_i + Q_{i-1} \]

\[ a_i\phi_i + b_i\phi_{i+1} + c_i\phi_{i-1} = d_i \Rightarrow \phi_i = -\frac{b_i}{a_i + c_i P_{i-1}}\phi_{i+1} + \frac{d_i - c_i Q_{i-1}}{a_i + c_i P_{i-1}} \quad (10.33) \]

比较式 10.32 与式 10.33，可得 \(P_i\) 与 \(Q_i\) 的如下递推关系：

\[ P_i = -\frac{b_i}{a_i + c_i P_{i-1}}, \quad Q_i = \frac{d_i - c_i Q_{i-1}}{a_i + c_i P_{i-1}}, \quad i = 1, 2, \ldots, N \quad (10.34) \]

对 \(i = 1\)，由式 10.30 计算 \(P_1\) 与 \(Q_1\) 为

\[ P_1 = -\frac{b_1}{a_1}, \quad Q_1 = \frac{d_1}{a_1} \quad (10.35) \]

对 \(i = N\)，由于 \(b_N = 0\)，可得

\[ b_N = 0 \Rightarrow P_N = 0 \Rightarrow \phi_N = Q_N \quad (10.36) \]

TDMA 求解算法可总结如下：(1) 用式 10.35 计算 \(P_1\) 与 \(Q_1\) 的值；(2) 对 \(i = 2, 3, \ldots, N\)，用前向递推由式 10.34 计算 \(P_i\) 与 \(Q_i\) 的值；(3) 按式 10.36 令 \(\phi_N = Q_N\)；(4) 对 \(i = N-1, N-2, \ldots, 3, 2, 1\)，用反向递推由式 10.32 计算 \(\phi_i\) 的值。

10.2.14 五对角矩阵算法（PDMA）

五对角矩阵算法（PDMA）[7–10] 用于求解具有五对角系数矩阵的代数方程组，源于将网格点 \(i\) 处的 \(\phi_i\) 与其两个上游（\(i - 1\) 与 \(i - 2\)）及两个下游（\(i + 1\) 与 \(i + 2\)）邻居的 \(\phi\) 值相关联的离散格式。对图 10.1 所示记号，一般代数方程写为

\[ a_i\phi_i + b_i\phi_{i+2} + c_i\phi_{i+1} + d_i\phi_{i-1} + e_i\phi_{i-2} = f_i, \quad i = 1, 2, 3, \ldots, N \quad (10.37) \]

满足

\[ d_1 = e_1 = e_2 = 0 \]

\[ b_{N-1} = b_N = c_N = 0 \quad (10.38) \]

对 \(i = 1\)，式 10.37 给出

\[ \phi_1 = -\frac{b_1}{a_1}\phi_3 - \frac{c_1}{a_1}\phi_2 + \frac{f_1}{a_1} \quad (10.39) \]

而对 \(i = 2\)，\(\phi_2\) 的值为

\[ \phi_2 = -\frac{a_1 b_2}{a_1 a_2 - d_2 c_1}\phi_4 - \frac{a_1 c_2 - b_1 d_2}{a_1 a_2 - d_2 c_1}\phi_3 + \frac{a_1 f_2 - d_2 f_1}{a_1 a_2 - d_2 c_1} \quad (10.40) \]

该过程可对其他 \(i\) 值继续，一般情况下 \(\phi_i\) 可以表达为

\[ \phi_i = P_i\phi_{i+2} + Q_i\phi_{i+1} + R_i, \quad i = 1, 2, 3, \ldots, N \quad (10.41) \]

利用式 10.41 计算 \(\phi_{i-1}\) 与 \(\phi_{i-2}\)，并将它们代入式 10.37，推导出 \(\phi_i\) 的方程为

\[ \phi_i = -\frac{b_i}{a_i + e_i P_{i-2} + (d_i + e_i Q_{i-2})Q_{i-1}}\phi_{i+2} - \frac{c_i + (d_i + e_i Q_{i-2})P_{i-1}}{a_i + e_i P_{i-2} + (d_i + e_i Q_{i-2})Q_{i-1}}\phi_{i+1} + \frac{f_i - e_i R_{i-2} - (d_i + e_i Q_{i-2})R_{i-1}}{a_i + e_i P_{i-2} + (d_i + e_i Q_{i-2})Q_{i-1}} \quad (10.42) \]

比较式 10.41 与式 10.42，\(P_i, Q_i, R_i\) 为

\[ P_i = -\frac{b_i}{a_i + e_i P_{i-2} + (d_i + e_i Q_{i-2})Q_{i-1}} \]

\[ Q_i = -\frac{c_i + (d_i + e_i Q_{i-2})P_{i-1}}{a_i + e_i P_{i-2} + (d_i + e_i Q_{i-2})Q_{i-1}} \quad (10.43) \]

\[ R_i = \frac{f_i - e_i R_{i-2} - (d_i + e_i Q_{i-2})R_{i-1}}{a_i + e_i P_{i-2} + (d_i + e_i Q_{i-2})Q_{i-1}} \]

其中 \(i = 1\) 与 \(i = 2\) 时取值为

\[ P_1 = -\frac{b_1}{a_1}, \quad Q_1 = -\frac{c_1}{a_1}, \quad R_1 = \frac{f_1}{a_1} \]

\[ P_2 = -\frac{b_2}{a_2 + d_2 Q_1}, \quad Q_2 = -\frac{c_2 + d_2 P_1}{a_2 + d_2 Q_1}, \quad R_2 = \frac{f_2 - d_2 R_1}{a_2 + d_2 Q_1} \quad (10.44) \]

由于 \(b_{N-1} = b_N = c_N = 0\)，则 \(P_{N-1} = P_N = Q_N = 0\)。因此，\(\phi_{N-1}\) 与 \(\phi_N\) 的方程为

\[ \phi_N = R_N \]

\[ \phi_{N-1} = Q_{N-1}\phi_N + R_{N-1} \quad (10.45) \]

PDMA 求解算法可总结如下：(1) 用式 10.44 计算 \(P_1, Q_1, R_1, P_2, Q_2, R_2\) 的值；(2) 对 \(i = 3, 4, \ldots, N\)，用前向递推由式 10.43 计算 \(P_i, Q_i, R_i\) 的值；(3) 由式 10.45 计算 \(\phi_N\) 与 \(\phi_{N-1}\)；(4) 对 \(i = N-2, \ldots, 3, 2, 1\)，用反向递推由式 10.41 计算 \(\phi_i\) 的值。

10.3 迭代法（Iterative Methods）

直接法通常不适合求解大型方程组，尤其是在系数矩阵稀疏（非零元素少）时更是如此。当线性化的方程组是非线性的、系数依赖于解、或者涉及时间相关问题时，这一情况更为突出——这正是求解流体流动问题所遇到的方程类型。

相比之下，迭代法对这些问题更具吸引力，因为线性化系统的求解成为迭代求解过程的一部分。再加上这种方法相对于直接法在计算机存储与计算代价上的低要求，使其成为首选。本节首先介绍基本迭代方法，附带对其缺陷的常用补救手段——多重网格算法的评述。前面介绍的高斯消去法与 LU 分解直接法，仅为澄清理解迭代法所需的基本数值过程。

为统一表述这些方法，系数矩阵将写成如下形式：

\[ \mathbf{A} = \mathbf{D} + \mathbf{L} + \mathbf{U} \quad (10.46) \]

其中 \(\mathbf{D}, \mathbf{L}, \mathbf{U}\) 分别是对角矩阵、严格下三角矩阵与严格上三角矩阵。

对形如 \(\mathbf{A}\boldsymbol{\phi} = \mathbf{b}\) 的线性系统，迭代法计算一系列解 \(\boldsymbol{\phi}^{(n)}\)，在满足一定条件时收敛于精确解 \(\boldsymbol{\phi}\)。因此，为求解，首先选定一个起点（即选取 \(\boldsymbol{\phi}^{(0)}\) 作为初始条件或初始猜测），再发展一种迭代过程，由先前算得的 \(\boldsymbol{\phi}^{(n-1)}\) 场计算 \(\boldsymbol{\phi}^{(n)}\)。

"不动点"迭代总可以通过将矩阵 \(\mathbf{A}\) 分解为

\[ \mathbf{A} = \mathbf{M} - \mathbf{N} \quad (10.47) \]

与上述系统关联。利用该分解，式 10.1 可改写为

\[ (\mathbf{M} - \mathbf{N})\boldsymbol{\phi} = \mathbf{b} \quad (10.48) \]

采用不动点迭代的求解过程，式 10.48 变为

\[ \mathbf{M}\boldsymbol{\phi}^{(n)} = \mathbf{N}\boldsymbol{\phi}^{(n-1)} + \mathbf{b} \quad (10.49) \]

可进一步写成

\[ \boldsymbol{\phi}^{(n)} = \mathbf{B}\boldsymbol{\phi}^{(n-1)} + \mathbf{C}\mathbf{b}, \quad n = 1, 2, \ldots \quad (10.50) \]

其中 \(\mathbf{B} = \mathbf{M}^{-1}\mathbf{N}\)，\(\mathbf{C} = \mathbf{M}^{-1}\)。这些矩阵的不同选择定义不同的迭代方法。

在开始描述各种迭代方法之前，先给出保证收敛所需的一组最基本的迭代方法特性。

A. 迭代方程在收敛时可写为

\[ \boldsymbol{\phi} = \mathbf{B}\boldsymbol{\phi} + \mathbf{C}\mathbf{b} \quad (10.51) \]

整理后变为

\[ \mathbf{C}^{-1}(\mathbf{I} - \mathbf{B})\boldsymbol{\phi} = \mathbf{b} \quad (10.52) \]

将式 10.52 与式 10.1 比较，可得系数矩阵为

\[ \mathbf{A} = \mathbf{C}^{-1}(\mathbf{I} - \mathbf{B}) \quad (10.53) \]

或等价地写为

\[ \mathbf{B} + \mathbf{C}\mathbf{A} = \mathbf{I} \quad (10.54) \]

这一矩阵之间的关系保证：一旦达到精确解，后续迭代不会改变它。

B. 从某个猜测 \(\boldsymbol{\phi}^{(0)} \neq \boldsymbol{\phi}\) 出发，方法应保证 \(\boldsymbol{\phi}^{(n)}\) 在 \(n\) 增大时收敛到 \(\boldsymbol{\phi}\)。由于 \(\boldsymbol{\phi}^{(n)}\) 可以用 \(\boldsymbol{\phi}^{(0)}\) 表示为

\[ \boldsymbol{\phi}^{(n)} = \mathbf{B}^n\boldsymbol{\phi}^{(0)} + \sum_{i=0}^{n-1}\mathbf{B}^i\mathbf{C}\mathbf{b} \quad (10.55) \]

则上述条件成立要求 \(\mathbf{B}\) 满足

\[ \lim_{n \to \infty}\mathbf{B}^n = \lim_{n \to \infty}\underbrace{\mathbf{B} \cdot \mathbf{B} \cdot \mathbf{B} \cdots \mathbf{B}}_{n \text{ times}} = 0 \quad (10.56) \]

式 10.56 意味着 \(\mathbf{B}\) 的谱半径应小于 1，即

\[ \rho(\mathbf{B}) < 1 \quad (10.57) \]

该条件保证迭代方法是自校正的，即对解向量 \(\boldsymbol{\phi}\) 中任何不利地引入的误差都具有鲁棒性。

通过定义解的误差 \(e^{(n)}\) 为精确值与第 \(n\) 次迭代值之差，可对该条件获得更深入的认识：

\[ e^{(n)} = \boldsymbol{\phi}^{(n)} - \boldsymbol{\phi}, \quad e^{(n-1)} = \boldsymbol{\phi}^{(n-1)} - \boldsymbol{\phi} \quad (10.58) \]

从式 10.50 减去式 10.51，并利用式 10.58 的定义，可得第 \(n\) 次与第 \(n - 1\) 次迭代误差之间的关系：

\[ e^{(n)} = \mathbf{B}e^{(n-1)} \quad (10.59) \]

因此，方法收敛需满足

\[ \lim_{n \to \infty} e^{(n)} = 0 \quad (10.60) \]

为将式 10.60 翻译为有意义的条件，假设 \(\mathbf{B}\) 的特征向量是完备的、构成完整集合，即它们形成 \(\mathbb{R}^N\) 的一组基。这种情况下，\(e\) 可以表示为 \(\mathbf{B}\) 的 \(N\) 个特征向量 \(v_i\) 的线性组合：

\[ e = \sum_{i=1}^{N}\alpha_i v_i \quad (10.61) \]

其中每个特征向量满足

\[ \mathbf{B}v_i = \lambda_i v_i \quad (10.62) \]

\(\lambda_i\) 是特征向量 \(v_i\) 对应的特征值。从第一次迭代开始，式 10.59 给出

\[ e^{(1)} = \mathbf{B}e^{(0)} = \mathbf{B}\sum_{i=1}^{N}\alpha_i v_i = \sum_{i=1}^{N}\alpha_i(\mathbf{B}v_i) = \sum_{i=1}^{N}\alpha_i\lambda_i v_i \quad (10.63) \]

对第二次迭代，误差为

\[ e^{(2)} = \mathbf{B}e^{(1)} = \mathbf{B}\sum_{i=1}^{N}\alpha_i\lambda_i v_i = \sum_{i=1}^{N}\alpha_i\lambda_i(\mathbf{B}v_i) = \sum_{i=1}^{N}\alpha_i\lambda_i^2 v_i \quad (10.64) \]

继续该过程，由归纳法容易得到

\[ e^{(n)} = \sum_{i=1}^{N}\alpha_i\lambda_i^n v_i \quad (10.65) \]

因此，为使迭代过程在 \(n\) 趋于无穷时收敛，所有特征值的绝对值都应小于 1。若其中任何一个大于 1，则误差将趋于无穷。这解释了上面提到的矩阵 \(\mathbf{B}\) 的谱半径 \(\rho\) 的重要性，定义为

\[ \rho(\mathbf{B}) = \max_{i=1}^{N}(\lambda_i) \quad (10.66) \]

迭代法的收敛速度通过减小迭代矩阵的谱半径来加快。这是迭代技术的核心。

C. 迭代法需要某种形式的停机判据。常用的判据基于残差的某种范数变化，残差定义为

\[ r^{(n)} = \mathbf{A}\boldsymbol{\phi}^{(n)} - \mathbf{b} \quad (10.67) \]

一种判据是找出整个区域内最大的残差并要求其值小于某个阈值 \(\varepsilon\) 来宣布解已收敛：

\[ \max_{i=1}^{N}\left|b_i - \sum_{j=1}^{N}a_{ij}\phi_j^{(n)}\right| \leq \varepsilon \quad (10.68) \]

或要求均方根残差小于 \(\varepsilon\)：

\[ \frac{\sum_{i=1}^{N}\left[b_i - \sum_{j=1}^{N}a_{ij}\phi_j^{(n)}\right]^2}{N} \leq \varepsilon \quad (10.69) \]

另一种可能的判据是连续两次迭代之间最大归一化差值低于 \(\varepsilon\)。该条件写为

\[ \max_{i=1}^{N}\left|\frac{\phi_i^{(n)} - \phi_i^{(n-1)}}{\phi_i^{(n)}}\right| \cdot 100 \leq \varepsilon \quad (10.70) \]

10.3.1 Jacobi 方法（Jacobi Method）

Jacobi 方法也许是求解线性方程组最简单的迭代方法，在图 10.2 中以图形方式给出。

考虑式 10.1 描述的方程组，若对角元素非零，则第一个方程可用于求解 \(\phi_1\)，第二个方程用于求解 \(\phi_2\)，以此类推。求解过程首先为未知向量 \(\boldsymbol{\phi}\) 赋初值。这些猜测值用于计算新估计，从 \(\phi_1\) 开始，再到 \(\phi_2\)，计算继续进行，直至 \(\phi_N\) 的新估计被算出。这代表一次迭代。所得到的结果被视为下一次迭代的新猜测，求解过程重复。迭代一直进行，直到连续两次迭代之间预测值的变化降至某个微小值或满足预设的收敛判据。一旦满足，收敛解即被获得。在该方法中，给定某个当前估计 \(\boldsymbol{\phi}^{(n-1)}\)，按下式获得更新：

\[ \phi_j^{(n)} = \frac{1}{a_{ii}}\left(b_i - \sum_{\substack{j=1\\j \neq i}}^{N} a_{ij}\phi_j^{(n-1)}\right), \quad i = 1, 2, 3, \ldots, N \quad (10.71) \]

式 10.71 表明，一次迭代中得到的值在同一次迭代的后续计算中并不使用，而是保留给下一次迭代。用矩阵表示，式 10.71 的展开形式由式 10.72 给出。求解 \(\boldsymbol{\phi}^{(n)}\) 后，式 10.72 化为式 10.73。

利用式 10.46，式 10.73 可更简洁地写为

\[ \boldsymbol{\phi}^{(n)} = -\mathbf{D}^{-1}(\mathbf{L} + \mathbf{U})\boldsymbol{\phi}^{(n-1)} + \mathbf{D}^{-1}\mathbf{b} \quad (10.74) \]

Jacobi 方法在 \(\rho\left(-\mathbf{D}^{-1}(\mathbf{L} + \mathbf{U})\right) < 1\) 时收敛。这一条件对一大类矩阵满足，包括对角占优矩阵，其系数满足

\[ \sum_{\substack{j=1\\j \neq i}}^{N}|a_{ij}| \leq |a_{ii}|, \quad i = 1, 2, 3, \ldots, N \quad (10.75) \]

10.3.2 Gauss-Seidel 方法（Gauss-Seidel Method）

Gauss-Seidel 方法是 Jacobi 方法的一种更流行的变体，其收敛特性更好。它在内存方面开销略低，因为不需要将新估计存放在单独的数组中，而是使用 \(\boldsymbol{\phi}\) 的最新估计进行计算。Gauss-Seidel 方法的迭代公式如图 10.3 所示，给出为

\[ \phi_i^{(n)} = \frac{1}{a_{ii}}\left(b_i - \sum_{j=1}^{i-1}a_{ij}\phi_j^{(n)} - \sum_{j=i+1}^{N}a_{ij}\phi_j^{(n-1)}\right), \quad i = 1, 2, 3, \ldots, N \quad (10.76) \]

用矩阵形式，式 10.76 写为

\[ \boldsymbol{\phi}^{(n)} = -(\mathbf{D} + \mathbf{L})^{-1}\mathbf{U}\boldsymbol{\phi}^{(n-1)} + (\mathbf{D} + \mathbf{L})^{-1}\mathbf{b} \quad (10.77) \]

实际上，Gauss-Seidel 方法在迭代中使用最新值，具体而言是所有 \(j < i\) 的 \(\phi_j^{(n)}\) 值，因为当计算 \(\phi_i\) 时，\(\phi_1, \phi_2, \phi_3, \ldots, \phi_{i-1}\) 在当前迭代下已经算得。这种方法也节省内存，因为新值总是覆盖旧值。Gauss-Seidel 迭代在

\[ \rho\left(-(\mathbf{D} + \mathbf{L})^{-1}\mathbf{U}\right) < 1 \quad (10.78) \]

时收敛。虽然在某些情况下 Jacobi 方法收敛更快，Gauss-Seidel 是首选方法。

例 2：对例 1 中的方程组施加 5 次 Gauss-Seidel 与 Jacobi 迭代，并利用精确解计算每次迭代的误差。精确解为 \(\boldsymbol{\phi} = (435/299, 408/299, 382/299, -19/299)\) 。以 \(\boldsymbol{\phi}' = [0, 0, 0, 0]\) 作为初始猜测。

Jacobi 方法中，以前一次迭代的值加撇号表示，待求解的方程为 \(\phi_1 = \frac{1}{3}(\phi'_2 + 3)\)，

\[ \phi_2 = \frac{1}{6}(2\phi'_1 + \phi'_3 + 4) \]

，

\[ \phi_3 = \frac{1}{6}(2\phi'_2 + \phi'_4 + 5) \]

，\(\phi_4 = \frac{1}{7}(2\phi'_3 - 3)\)，误差为 \(e = |\phi_{\text{exact}} - \phi_{\text{computed}}|\)。第一次迭代的解为：\(\phi_1 = 1\)，\(e_1 = 0.4548\)；\(\phi_2 = 0.6667\)，\(e_2 = 0.6978\)；\(\phi_3 = 0.8333\)，\(e_3 = 0.4443\)；\(\phi_4 = -0.4286\)，\(e_4 = 0.3650\)。计算继续以相同方式进行，所得结果作为新的猜测。前五次迭代的结果列于表 10.1。

Gauss-Seidel 方法中，待求解的方程为 \(\phi_1 = \frac{1}{3}(\phi'_2 + 3)\)，

\[ \phi_2 = \frac{1}{6}(2\phi_1 + \phi'_3 + 4) \]

，

\[ \phi_3 = \frac{1}{6}(2\phi_2 + \phi'_4 + 5) \]

，\(\phi_4 = \frac{1}{7}(2\phi_3 - 3)\)。第一次迭代的解为：\(\phi_1 = 1\)，\(e_1 = 0.4548\)；\(\phi_2 = 1\)，\(e_2 = 0.3645\)；\(\phi_3 = 1.1667\)，\(e_3 = 0.1109\)；\(\phi_4 = -0.09523\)，\(e_4 = 0.03169\)。前五次迭代的结果列于表 10.2。

10.3.3 预条件与迭代方法（Preconditioning and Iterative Methods）

迭代法的收敛速度取决于迭代矩阵 \(\mathbf{B}\) 的谱性质，而后者又取决于系数矩阵。基于此，迭代法寻求将原方程组变换为一个等价但具有更好谱性质的方程组。在这些条件下，等价系统的特征值更为聚集，从而迭代求解可以比原系统更快地获得。一个预条件子（preconditioner）被定义为实现这种变换的矩阵。

预条件矩阵 \(\mathbf{P}\) 定义为使系统

\[ \mathbf{P}^{-1}\mathbf{A}\boldsymbol{\phi} = \mathbf{P}^{-1}\mathbf{b} \quad (10.79) \]

与原系统 \(\mathbf{A}\boldsymbol{\phi} = \mathbf{b}\) 具有相同的解，但其系数矩阵 \(\mathbf{P}^{-1}\mathbf{A}\) 的谱性质更有利于迭代收敛。在定义预条件子 \(\mathbf{P}\) 时，困难在于找到一个既近似 \(\mathbf{A}^{-1}\) 又易于求逆（即易于求 \(\mathbf{P}^{-1}\)）且代价合理的矩阵。

再次写出式 10.47，但用 \(\mathbf{P}\) 替换 \(\mathbf{M}\)（即 \(\mathbf{M} = \mathbf{P}\) 且 \(\mathbf{A} = \mathbf{P} - \mathbf{N}\)），则关联的不动点迭代系统为

\[ \boldsymbol{\phi}^{(n)} = \mathbf{B}\boldsymbol{\phi}^{(n-1)} + \mathbf{C}\mathbf{b} \]

\[ = \mathbf{P}^{-1}\mathbf{N}\boldsymbol{\phi}^{(n-1)} + \mathbf{P}^{-1}\mathbf{b} \quad (10.80) \]

\[ = \mathbf{P}^{-1}(\mathbf{P} - \mathbf{A})\boldsymbol{\phi}^{(n-1)} + \mathbf{P}^{-1}\mathbf{b} \]

\[ = \left(\mathbf{I} - \mathbf{P}^{-1}\mathbf{A}\right)\boldsymbol{\phi}^{(n-1)} + \mathbf{P}^{-1}\mathbf{b} \]

以残差形式可写为

\[ \boldsymbol{\phi}^{(n)} = \left(\mathbf{I} - \mathbf{P}^{-1}\mathbf{A}\right)\boldsymbol{\phi}^{(n-1)} + \mathbf{P}^{-1}\mathbf{b} \]

\[ = \boldsymbol{\phi}^{(n-1)} + \mathbf{P}^{-1}\left(\mathbf{b} - \mathbf{A}\boldsymbol{\phi}^{(n-1)}\right) \quad (10.81) \]

\[ = \boldsymbol{\phi}^{(n-1)} + \mathbf{P}^{-1}r^{(n-1)} \]

从上述两式可清楚地看出，迭代过程本质上是对分解 \(\mathbf{A} = \mathbf{P} - \mathbf{N}\) 所关联的预条件系统的不动点迭代，新的谱性质为

\[ \rho\left(\mathbf{I} - \mathbf{P}^{-1}\mathbf{A}\right) < 1 \quad (10.82) \]

作为对比，Jacobi（\(J\)）和 Gauss-Seidel（\(GS\)）方法的预条件矩阵就是

\[ \mathbf{P}_J = \mathbf{D}, \quad \mathbf{P}_{GS} = \mathbf{D} + \mathbf{L} \quad (10.83) \]

其中 \(\mathbf{D}\) 和 \(\mathbf{L}\) 分别是矩阵 \(\mathbf{A}\) 的对角部分与下三角部分。

因此，预条件是对原系统的一种操作以改善其谱性质，预条件矩阵 \(\mathbf{P}\) 用于关联的迭代过程。如下几节所述，可以构造系数定义更为复杂的更先进的预条件矩阵。

10.3.4 矩阵分解技术（Matrix Decomposition Techniques）

Gauss-Seidel 与 Jacobi 方法的低收敛率是发展更快迭代技术的主要动机。加速求解器收敛率和发展迭代方法的一种途径是使用更先进的预条件子。为此目的，一种简单而高效的方法是对原系数矩阵 \(\mathbf{A}\) 进行不完全分解。强调"不完全"是关键的，因为将 \(\mathbf{A}\) 完全分解为一个下三角矩阵 \(\mathbf{L}\) 与一个上三角矩阵 \(\mathbf{U}\) 等价于直接求解法，在内存需求（填充与稀疏性丧失）和计算代价方面都非常昂贵。

10.3.5 不完全 LU（ILU）分解（Incomplete LU (ILU) Decomposition）

如例 1 所示，\(\mathbf{L}\) 与 \(\mathbf{U}\) 矩阵的结果在原矩阵 \(\mathbf{A}\) 中为零的位置上出现非零元素（这被称为填充，fill-in）。因此，如果对 \(\mathbf{A}\) 进行不完全 LU（ILU）分解，使所得下三角矩阵 \(\mathbf{L}\) 与上三角矩阵 \(\mathbf{U}\) 的非零结构与 \(\mathbf{A}\) 的下三角与上三角部分相同，则有

\[ \mathbf{A} = \mathbf{L}\mathbf{U} + \mathbf{R} \quad (10.84) \]

其中 \(\mathbf{R}\) 是分解过程的残差。\(\mathbf{L}\) 与 \(\mathbf{U}\) 矩阵是稀疏的（与 \(\mathbf{A}\) 结构相同），比由完全分解得到的更易处理。然而，它们的乘积只是 \(\mathbf{A}\) 的近似，必须采用迭代求解过程来求解方程组。求解过程的第一步是将式 10.1 改写为

\[ \mathbf{A}\boldsymbol{\phi} = \mathbf{b} \Rightarrow 0 = \mathbf{b} - \mathbf{A}\boldsymbol{\phi} \Rightarrow (\mathbf{A} - \mathbf{R})\boldsymbol{\phi} = (\mathbf{A} - \mathbf{R})\boldsymbol{\phi} + (\mathbf{b} - \mathbf{A}\boldsymbol{\phi}) \quad (10.85) \]

用上标 \((n-1)\) 表示前一次迭代得到的值，\((n)\) 表示当前迭代得到的值，迭代过程通过将式 10.85 改写为以下形式得到：

\[ (\mathbf{A} - \mathbf{R})\boldsymbol{\phi}^{(n)} = (\mathbf{A} - \mathbf{R})\boldsymbol{\phi}^{(n-1)} + \mathbf{b} - \mathbf{A}\boldsymbol{\phi}^{(n-1)} \quad (10.86) \]

因此当前迭代的 \(\boldsymbol{\phi}^{(n)}\) 值可由前一次迭代得到的 \(\boldsymbol{\phi}^{(n-1)}\) 值求得。式 10.86 通常以残差形式求解——其中迭代 \((n)\) 的解 \(\boldsymbol{\phi}^{(n)}\) 表示为迭代 \((n-1)\) 的解 \(\boldsymbol{\phi}^{(n-1)}\) 加上一个修正量 \(\boldsymbol{\phi}'^{(n)}\)，即

\[ \boldsymbol{\phi}^{(n)} = \boldsymbol{\phi}^{(n-1)} + \boldsymbol{\phi}'^{(n)} \quad (10.87) \]

因此式 10.86 变为

\[ (\mathbf{A} - \mathbf{R})\boldsymbol{\phi}'^{(n)} = \mathbf{b} - \mathbf{A}\boldsymbol{\phi}^{(n-1)} \quad (10.88) \]

一旦求出 \(\boldsymbol{\phi}'^{(n)}\)，用式 10.87 在每次迭代中更新 \(\boldsymbol{\phi}\)。ILU 分解可以通过高斯消去法执行，同时丢弃预设位置上的某些非对角元素。决定元素被丢弃的位置便产生不同的 ILU 近似。

10.3.6 无填充的不完全 LU 分解 ILU(0)（Incomplete LU Factorization with no Fill-in ILU(0)）

ILU 分解技术存在多种变体，最简单的一种记为 ILU(0) [15–17]。在 ILU(0) 中，\(\mathbf{L}\) 与 \(\mathbf{U}\) 矩阵的组合零元素模式被认为恰好与原矩阵 \(\mathbf{A}\) 的零元素模式相同。使用高斯消去法，计算按完全 LU 分解的方式进行，但若过程中出现的任何新的非零元素（\(\ell_{ij}\) 与 \(u_{ij}\)）落在原矩阵 \(\mathbf{A}\) 中为零的位置上，则将其丢弃。因此，\(\mathbf{L}\) 与 \(\mathbf{U}\) 矩阵的组合所具有的非零元素数量与原矩阵 \(\mathbf{A}\) 相同。这种方法消除了分解稀疏矩阵时通常出现的填充问题（即在原矩阵为零的位置产生非零元素）。然而在过程中，精度会下降，从而需要更多的迭代次数才能收敛。为弥补这一不足，已经发展出精度更高、通常更高效更可靠的 ILU 分解方法。这些方法以所允许的填充级别区分，记为 ILU(\(p\))，其中 \(p\) 代表填充阶数。填充级别越高，ILU 分解步骤的代价就越大。此外，当在多重网格方法中使用时（将在后文解释），ILU(0) 方法作为光滑子已足够。因此更高级别的方法不在此介绍，对详细内容感兴趣的读者可参阅 Saad 的著作 [12]。

下面给出一种 ILU(0) 分解算法，该算法假设 \(\mathbf{L}\) 为单位下三角矩阵，并且同一矩阵 \(\mathbf{A}\) 被用来存储单位下三角矩阵 \(\mathbf{L}\) 与上三角矩阵 \(\mathbf{U}\) 的元素。

10.3.7 ILU(0) 分解算法（ILU(0) Factorization Algorithm）

For k = 1 to N-1
  {
    For i = k + 1 to N and if a_ik ≠ 0 Do:
      {
        a_ik = a_ik / a_kk  (ℓ values)
        {
          For j = k + 1 to N and if a_ij ≠ 0 Do:
            a_ij = a_ij - a_ik * a_kj  (u values)
        }
      }
  }

需要提及的是，对称正定矩阵的 ILU 分解记为不完全 Cholesky 分解。这种情况下分解只由下三角（或上三角）部分构成，对原矩阵的近似写为

\[ \mathbf{L}\mathbf{L}^T \approx \mathbf{A} \quad (10.89) \]

其中 \(\mathbf{L}\) 是分解得到的稀疏下三角矩阵（\(\mathbf{L}\) 的近似），预条件矩阵 \(\mathbf{P}\) 由下式给出

\[ \mathbf{P} = \mathbf{L}\mathbf{L}^T \approx \mathbf{A} \quad (10.90) \]

10.3.8 ILU 分解预条件子（ILU Factorization Preconditioners）

基于不完全分解的预条件子构成了一类非常流行的预条件子。在直接法的讨论中已经表明，将稀疏矩阵 \(\mathbf{A}\) 分解为下三角矩阵与上三角矩阵的乘积可能导致大量填充。因为预条件子只需作为 \(\mathbf{A}^{-1}\) 的近似即可，所以只需寻找 \(\mathbf{A}\) 的近似分解形式 \(\mathbf{A} \approx \mathbf{L}\mathbf{U}\)。选取 \(\mathbf{P} = \mathbf{L}\mathbf{U}\) 也带来预条件矩阵逆 \(\mathbf{P}^{-1}\) 的高效求取——因为该求逆可以通过前向代入和回代轻松实现（如前文所述），其中精确的 \(\mathbf{L}\) 与 \(\mathbf{U}\) 被替换为相应的近似 \(\mathbf{L}\) 与 \(\mathbf{U}\)。

对 ILU(0) 方法，不完全分解模仿原矩阵的非零元素稀疏模式，使得预条件子恰好具有原矩阵的大小。为减少所需的存储，Pommerell 引入了一种简化的 ILU 版本，称为对角 ILU（DILU）[18]。在 DILU 中，非对角元素的填充被消除（即矩阵的上、下三角部分保持不变），仅对角元素被修改。

在这种情况下，可以将预条件子写成

\[ \mathbf{P} = (\mathbf{D}^* + \mathbf{L})\mathbf{D}^{*-1}(\mathbf{D}^* + \mathbf{U}) \quad (10.91) \]

其中 \(\mathbf{L}\) 与 \(\mathbf{U}\) 是 \(\mathbf{A}\) 的下三角与上三角分解，\(\mathbf{D}^*\) 现在是一个真正的对角矩阵，不同于 \(\mathbf{A}\) 的对角部分。\(\mathbf{D}^*\) 矩阵的定义方式如下：使式 10.91 中矩阵乘积的对角线等于 \(\mathbf{A}\) 的对角线。

10.3.9 DILU 中 \(\mathbf{D}^*\) 的计算算法（Algorithm for the Calculation of \(\mathbf{D}^*\) in the DILU Method）

For i = 1 to N Do:
  {
    d_ii = a_ii
  }
For i = 1 to N Do:
  {
    For j = i + 1 to N and if a_ij ≠ 0, a_ji ≠ 0 Do:
      {
        d_jj = d_jj - (a_ji / d_ii) * a_ij
      }
  }

这种情况下预条件子的逆定义为

\[ \mathbf{P} = (\mathbf{D}^* + \mathbf{L})\mathbf{D}^{*-1}(\mathbf{D}^* + \mathbf{U}) = \mathbf{L}\mathbf{U}; \quad \mathbf{L} = (\mathbf{D}^* + \mathbf{L})\mathbf{D}^{*-1}, \quad \mathbf{U} = (\mathbf{D}^* + \mathbf{U}) \]

或

\[ \mathbf{P} = (\mathbf{D}^* + \mathbf{L})\left(\mathbf{I} + \mathbf{D}^{*-1}\mathbf{U}\right) = \mathbf{L}\mathbf{U}; \quad \mathbf{L} = (\mathbf{D}^* + \mathbf{L}), \quad \mathbf{U} = \left(\mathbf{I} + \mathbf{D}^{*-1}\mathbf{U}\right) \quad (10.92) \]

在求解 \(\mathbf{P}\boldsymbol{\phi}'^{(n+1)} = r^{(n)}\) 以求修正场 \(\boldsymbol{\phi}'^{(n+1)} = \mathbf{P}^{-1}r^{(n)}\) 时所需的预条件子之逆，可以通过以下前向与回代算法方便地计算。

10.3.10 DILU 方法的前向与回代算法（Forward and Backward Solution Algorithm with the DILU Method）

For i = 1 to N Do:
  {
    For j = 1 to i-1 Do:
      {
        t_i = d_ii^{-1} * (r_i - ℓ_ij * t_j)
      }
  }
For i = N to 1 Do:
  {
    For j = i + 1 to N Do:
      {
        φ'_i = t_i - d_ii^{-1} * (u_ij * t_j)
      }
  }

DILU 的显著优势在于（除其递推形式之外），它仅需要额外的一对对角线存储。

10.3.11 求解代数系统的梯度法（Gradient Methods for Solving Algebraic Systems）

求解线性代数方程组的另一类迭代过程是梯度法，包括最速下降法与共轭梯度法。它们最初是为系数矩阵 \(\mathbf{A}\) 对称正定（SPD）的情形而发展出来的，将问题重新表述为关于二次向量函数 \(Q(\boldsymbol{\phi})\) 的极小化问题：

\[ Q(\boldsymbol{\phi}) = \frac{1}{2}\boldsymbol{\phi}^T\mathbf{A}\boldsymbol{\phi} - \mathbf{b}^T\boldsymbol{\phi} + c \quad (10.93) \]

其中 \(c\) 是标量向量，其他变量如式 10.1 所定义。\(Q(\boldsymbol{\phi})\) 的极小值在 \(Q\) 关于 \(\boldsymbol{\phi}\) 的梯度为零时取得。向量场 \(Q(\boldsymbol{\phi})\) 在给定点 \(\boldsymbol{\phi}\) 处的梯度 \(Q'(\boldsymbol{\phi})\) 指向 \(Q(\boldsymbol{\phi})\) 增长最快的方向。经数学推导，梯度求得为

\[ Q'(\boldsymbol{\phi}) = \frac{1}{2}\mathbf{A}^T\boldsymbol{\phi} + \frac{1}{2}\mathbf{A}\boldsymbol{\phi} - \mathbf{b} \quad (10.94) \]

若 \(\mathbf{A}\) 对称，即 \(\mathbf{A} = \mathbf{A}^T\)，则式 10.94 意味着

\[ Q'(\boldsymbol{\phi}) = \mathbf{A}\boldsymbol{\phi} - \mathbf{b} \quad (10.95) \]

当 \(Q'(\boldsymbol{\phi}) = 0\) 时取得极小，由此得到

\[ Q'(\boldsymbol{\phi}) = 0 \Rightarrow \mathbf{A}\boldsymbol{\phi} = \mathbf{b} \quad (10.96) \]

因此极小化 \(Q(\boldsymbol{\phi})\) 等价于求解式 10.1，极小化问题的解就是线性方程组的解。

为使函数 \(Q(\boldsymbol{\phi})\) 具有全局极小，系数矩阵 \(\mathbf{A}\) 必须正定，即对所有 \(\boldsymbol{\phi} \neq 0\)，应满足 \(\boldsymbol{\phi}^T\mathbf{A}\boldsymbol{\phi} \geq 0\) 。这一要求可以通过考察精确解 \(\boldsymbol{\phi}\) 与其当前估计 \(\boldsymbol{\phi}^{(n)}\) 之间的关系来建立。若 \(e = \boldsymbol{\phi}^{(n)} - \boldsymbol{\phi}\) 表示精确解与当前估计之差，则式 10.93 给出

\[ Q(\boldsymbol{\phi} + e) = \frac{1}{2}(\boldsymbol{\phi} + e)^T\mathbf{A}(\boldsymbol{\phi} + e) - \mathbf{b}^T(\boldsymbol{\phi} + e) + c \]

\[ = \frac{1}{2}\boldsymbol{\phi}^T\mathbf{A}\boldsymbol{\phi} + \frac{1}{2}e^T\mathbf{A}\boldsymbol{\phi} + \frac{1}{2}\boldsymbol{\phi}^T\mathbf{A}e + \frac{1}{2}e^T\mathbf{A}e - \mathbf{b}^T\boldsymbol{\phi} - \mathbf{b}^T e + c \]

\[ = \underbrace{\frac{1}{2}\boldsymbol{\phi}^T\mathbf{A}\boldsymbol{\phi} - \mathbf{b}^T\boldsymbol{\phi} + c}_{Q(\boldsymbol{\phi})} + \underbrace{\left[e^T\mathbf{A}\boldsymbol{\phi} + \boldsymbol{\phi}^T\mathbf{A}e - \mathbf{b}^T e\right]}_{\mathbf{e}^T\mathbf{b} - \mathbf{b}^T\mathbf{e} = 0} + \frac{1}{2}e^T\mathbf{A}e \]

\[ = Q(\boldsymbol{\phi}) + \frac{1}{2}e^T\mathbf{A}e \quad (10.97) \]

这表明：若 \(\mathbf{A}\) 正定，则第二项除 \(e = 0\) 外始终为正，而 \(e = 0\) 时所需的解已经得到。此外，当 \(\mathbf{A}\) 正定时，其所有特征值为正，函数 \(Q(\boldsymbol{\phi})\) 具有唯一的极小值。

因此，对于对称正定矩阵，可以导出收敛的 \(\boldsymbol{\phi}^{(n)}\) 序列，使得

\[ \boldsymbol{\phi}^{(n+1)} = \boldsymbol{\phi}^{(n)} + \alpha^{(n)} \delta\boldsymbol{\phi}^{(n)} \quad (10.98) \]

其中 \(\alpha^{(n)}\) 是某个松弛因子，\(\delta\boldsymbol{\phi}^{(n)}\) 与每步迭代中使所述函数极小化所需的修正量相关。这可以通过多种方式实现，从而产生不同的方法。

10.3.12 最速下降法（The Method of Steepest Descent）

最速下降法用于求解式 10.1 给定形式的线性方程组，基于极小化式 10.93 给出的二次型。若 \(\boldsymbol{\phi}\) 是一个以标量 \(\phi\) 为分量的"一维"向量，则 \(Q(\boldsymbol{\phi})\) 表示一条抛物线。从某点 \(\phi^{(0)}\) 起迭代求抛物线函数的极小值，需要沿抛物线下行直至到达极小点。

同样的思路被用于 \(N\) 维情形。此时 \(Q(\boldsymbol{\phi})\) 可被视为一个抛物面，从初始位置 \(\boldsymbol{\phi}^{(0)}\) 起迭代下移抛物面直至到达极小点。为加快收敛，步长序列

\[ \boldsymbol{\phi}^{(0)}, \boldsymbol{\phi}^{(1)}, \boldsymbol{\phi}^{(2)}, \ldots \]

应选取为具有最快的下降速率，即沿 \(-Q'(\boldsymbol{\phi})\) 方向。根据式 10.95，该方向亦由

\[ -Q'(\boldsymbol{\phi}) = \mathbf{b} - \mathbf{A}\boldsymbol{\phi} \quad (10.99) \]

给出。精确解为 \(\boldsymbol{\phi}\)，第 \(n\) 步的误差与残差（分别记为 \(e^{(n)}\) 与 \(r^{(n)}\)）为

\[ e^{(n)} = \boldsymbol{\phi}^{(n)} - \boldsymbol{\phi} \Rightarrow r^{(n)} = -\mathbf{A}e^{(n)} \]

\[ r^{(n)} = \mathbf{b} - \mathbf{A}\boldsymbol{\phi}^{(n)} = -Q'(\boldsymbol{\phi}^{(n)}) \quad (10.100) \]

沿最速下降方向线性移动，第 \(n + 1\) 步的 \(\boldsymbol{\phi}\) 可以按第 \(n\) 步的 \(\boldsymbol{\phi}\) 表达为

\[ \boldsymbol{\phi}^{(n+1)} = \boldsymbol{\phi}^{(n)} + \alpha^{(n)} r^{(n)} \quad (10.101) \]

使 \(Q(\boldsymbol{\phi})\) 极小的 \(\alpha^{(n)}\) 值应满足

\[ \frac{d}{d\alpha^{(n)}}Q(\boldsymbol{\phi}^{(n+1)}) = 0 \quad (10.102) \]

这可展开为

\[ \frac{d}{d\alpha^{(n)}}Q(\boldsymbol{\phi}^{(n+1)}) = 0 \Rightarrow \frac{d}{d\boldsymbol{\phi}}Q(\boldsymbol{\phi}^{(n+1)})\frac{d\boldsymbol{\phi}^{(n+1)}}{d\alpha^{(n)}} = 0 \Rightarrow r^{(n+1)T}r^{(n)} = 0 \quad (10.103) \]

这表明新步长方向应与旧步长方向正交。\(\alpha^{(n)}\) 的值由式 10.103 计算如下：

\[ r^{(n+1)T}r^{(n)} = 0 \Rightarrow \left(\mathbf{b} - \mathbf{A}\boldsymbol{\phi}^{(n+1)}\right)^T r^{(n)} = 0 \]

\[ \Rightarrow \left(\mathbf{b} - \mathbf{A}\left(\boldsymbol{\phi}^{(n)} + \alpha^{(n)} r^{(n)}\right)\right)^T r^{(n)} = 0 \]

\[ \Rightarrow \left(\mathbf{b} - \mathbf{A}\boldsymbol{\phi}^{(n)}\right)^T r^{(n)} = \alpha^{(n)}(\mathbf{A}r^{(n)})^T r^{(n)} \quad (10.104) \]

\[ \Rightarrow r^{(n)T}r^{(n)} = \alpha^{(n)} r^{(n)T}\mathbf{A}r^{(n)} \]

\[ \Rightarrow \alpha^{(n)} = \frac{r^{(n)T}r^{(n)}}{r^{(n)T}\mathbf{A}r^{(n)}} \]

最速下降算法可总结如下：

r^(0) = b - Aφ^(0)   (choose residual as starting direction)
iterate starting at n until convergence
  r^(n) = b - Aφ^(n)   (Compute the residual vector)
  α^(n) = (r^(n))^T r^(n) / ((r^(n))^T A r^(n))   (Compute the factor in the orthogonal direction)
  φ^(n+1) = φ^(n) + α^(n) r^(n)   (Obtain new φ)

如上所述，算法在每次迭代中需要执行两次矩阵—向量乘法。其中一次可以通过对式 10.101 两边乘以 \(-\mathbf{A}\) 再加 \(\mathbf{b}\) 而消去：

\[ \boldsymbol{\phi}^{(n+1)} = \boldsymbol{\phi}^{(n)} + \alpha^{(n)} r^{(n)} \Rightarrow \mathbf{b} - \mathbf{A}\boldsymbol{\phi}^{(n+1)} \]

\[ = \mathbf{b} - \mathbf{A}\left(\boldsymbol{\phi}^{(n)} + \alpha^{(n)} r^{(n)}\right) \Rightarrow r^{(n+1)} = r^{(n)} - \alpha^{(n)}\mathbf{A}r^{(n)} \quad (10.105) \]

第 1 步中 \(r^{(n)}\) 的方程仅用于计算 \(r^{(0)}\)，之后可以使用式 10.105。采用这种形式后，无需计算 \(\mathbf{A}\boldsymbol{\phi}^{(n)}\)，因为已被 \(\mathbf{A}r^{(n)}\) 取代。然而该方法的不足之处在于 \(\boldsymbol{\phi}^{(n)}\) 的值对残差缺乏反馈，可能由于舍入误差的累积导致求解收敛到一个不同于精确解的值。这一不足可以通过周期性地按原方程计算残差来解决。

10.3.13 共轭梯度法（The Conjugate Gradient Method）

最速下降法虽然保证收敛，但收敛速度较慢。这种慢收敛源于围绕局部极小的振荡，迫使方法在同一方向上反复搜索。为避免这种不期望的行为，每一次新搜索的方向应与之前搜索的方向不同 [19]。这可以通过选择一组 A-正交的搜索方向 \(\mathbf{d}^{(0)}, \mathbf{d}^{(1)}, \mathbf{d}^{(2)}, \ldots, \mathbf{d}^{(N-1)}\) 来实现。若两个向量 \(\mathbf{d}^{(n)}\) 与 \(\mathbf{d}^{(m)}\) 满足如下条件，则称它们为 A-正交：

\[ \mathbf{d}^{(n)T}\mathbf{A}\mathbf{d}^{(m)} = 0 \quad (10.106) \]

若在每个搜索方向上取合适的步长，则解将在 \(N\) 步后被找到。第 \(n + 1\) 步的选取满足

\[ \boldsymbol{\phi}^{(n+1)} = \boldsymbol{\phi}^{(n)} + \alpha^{(n)}\mathbf{d}^{(n)} \quad (10.107) \]

从上式两边减去 \(\boldsymbol{\phi}\)，可得误差方程为

\[ e^{(n+1)} = e^{(n)} + \alpha^{(n)}\mathbf{d}^{(n)} \quad (10.108) \]

将式 10.100 与式 10.108 结合，可得残差方程为

\[ r^{(n+1)} = -\mathbf{A}e^{(n+1)} \]

\[ = -\mathbf{A}\left(e^{(n)} + \alpha^{(n)}\mathbf{d}^{(n)}\right) \quad (10.109) \]

\[ = r^{(n)} - \alpha^{(n)}\mathbf{A}\mathbf{d}^{(n)} \]

式 10.109 表明，新的残差 \(r^{(n+1)}\) 仅为前一次残差与 \(\mathbf{A}\mathbf{d}^{(n)}\) 的线性组合。还要求 \(e^{(n+1)}\) 与 \(\mathbf{d}^{(n)}\) A-正交。这一新条件等价于沿搜索方向 \(\mathbf{d}^{(n)}\) 找到极小点。利用 \(e^{(n+1)}\) 与 \(\mathbf{d}^{(n)}\) 之间的 A-正交性条件及式 10.108，可推导出 \(\alpha^{(n)}\) 的表达式为

\[ \mathbf{d}^{(n)T}\mathbf{A}e^{(n+1)} = 0 \Rightarrow \mathbf{d}^{(n)T}\mathbf{A}\left(e^{(n)} + \alpha^{(n)}\mathbf{d}^{(n)}\right) = 0 \]

\[ \Rightarrow \alpha^{(n)} = -\frac{\mathbf{d}^{(n)T}\mathbf{A}e^{(n)}}{\mathbf{d}^{(n)T}\mathbf{A}\mathbf{d}^{(n)}} = \frac{\mathbf{d}^{(n)T}r^{(n)}}{\mathbf{d}^{(n)T}\mathbf{A}\mathbf{d}^{(n)}} \quad (10.110) \]

以上要求也意味着

\[ \mathbf{d}^{(n)T}\mathbf{A}e^{(n+1)} = 0 \Rightarrow \mathbf{d}^{(n)T}r^{(n+1)} = 0 \quad (10.111) \]

若搜索方向已知，则可以计算 \(\alpha^{(n)}\)。为推导搜索方向，假设它由如下形式的方程控制：

\[ \mathbf{d}^{(n+1)} = r^{(n+1)} + \beta^{(n)}\mathbf{d}^{(n)} \quad (10.112) \]

\(\mathbf{d}\) 向量的 A-正交性要求

\[ \mathbf{d}^{(n+1)T}\mathbf{A}\mathbf{d}^{(n)} = 0 \quad (10.113) \]

将式 10.112 中 \(\mathbf{d}^{(n+1)}\) 的值代入式 10.113，得到

\[ \beta^{(n)} = -\frac{r^{(n+1)T}\mathbf{A}\mathbf{d}^{(n)}}{\mathbf{d}^{(n)T}\mathbf{A}\mathbf{d}^{(n)}} \quad (10.114) \]

由式 10.109 可得 \(\mathbf{A}\mathbf{d}^{(n)}\) 的表达式为

\[ \mathbf{A}\mathbf{d}^{(n)} = \frac{1}{\alpha^{(n)}}\left(r^{(n)} - r^{(n+1)}\right) \quad (10.115) \]

联立式 10.110、10.114 与 10.115，得到

\[ \beta^{(n)} = \frac{r^{(n+1)T}\left(r^{(n+1)} - r^{(n)}\right)}{\mathbf{d}^{(n)T}r^{(n)}} \]

\[ = \frac{r^{(n+1)T}r^{(n+1)} - \overbrace{r^{(n+1)T}r^{(n)}}^{=0}}{\mathbf{d}^{(n)T}r^{(n)}} \quad (10.116) \]

\[ = \frac{r^{(n+1)T}r^{(n+1)}}{\mathbf{d}^{(n)T}r^{(n)}} \]

上式的分母可以进一步表示为

\[ \mathbf{d}^{(n)T}r^{(n)} = \left(r^{(n)} + \beta^{(n-1)}\mathbf{d}^{(n-1)}\right)^T r^{(n)} \]

\[ = r^{(n)T}r^{(n)} + \beta^{(n-1)}\underbrace{\mathbf{d}^{(n-1)T}r^{(n)}}_{=0} \quad (10.117) \]

\[ = r^{(n)T}r^{(n)} \]

利用式 10.116 与 10.117，\(\beta^{(n)}\) 的最终表达式为

\[ \beta^{(n)} = \frac{r^{(n+1)T}r^{(n+1)}}{r^{(n)T}r^{(n)}} \quad (10.118) \]

共轭梯度算法变为

d^(0) = r^(0) = b - Aφ^(0)   (choose residual as starting direction)
iterate starting at n until convergence
  α^(n) = (d^(n))^T r^(n) / ((d^(n))^T A d^(n))   (Choose factor in d direction)
  φ^(n+1) = φ^(n) + α^(n) d^(n)   (Obtain new φ)
  r^(n+1) = r^(n) - α^(n) A d^(n)   (calculate new residual)
  β^(n) = (r^(n+1))^T r^(n+1) / ((r^(n))^T r^(n))   (Calculate coefficient to conjugate residual)
  d^(n+1) = r^(n+1) + β^(n) d^(n)   (obtain new conjugated search direction)

CG 方法的收敛速度可通过预条件来提高。这可以通过将原方程组乘以预条件矩阵 \(\mathbf{P}\) 的逆 \(\mathbf{P}^{-1}\)（其中 \(\mathbf{P}\) 是对称正定矩阵）来实现，得到式 10.79。问题在于即使 \(\mathbf{P}\) 与 \(\mathbf{A}\) 都是对称的，\(\mathbf{P}^{-1}\mathbf{A}\) 也不一定对称。为绕过这个问题，使用 Cholesky 分解将 \(\mathbf{P}\) 写成

\[ \mathbf{P} = \mathbf{L}\mathbf{L}^T \quad (10.119) \]

为保证对称性，方程组写为

\[ \mathbf{L}^{-1}\mathbf{A}\mathbf{L}^{-T}\mathbf{L}^T\boldsymbol{\phi} = \mathbf{L}^{-1}\mathbf{b} \quad (10.120) \]

其中 \(\mathbf{L}^{-1}\mathbf{A}\mathbf{L}^{-T}\) 对称且正定。CG 方法可用于求解 \(\mathbf{L}^T\boldsymbol{\phi}\)，再由此求出 \(\boldsymbol{\phi}\)。然而，通过变量替换，\(\mathbf{L}\) 可以从方程中消去而不破坏对称性或影响方法的有效性。完成此步并采用 CG 方法所用的术语，可以得到预条件 CG 方法中的各步。

预条件 CG 方法可总结如下：

r^(0) = b - Aφ^(0) and d^(0) = P^(-1) r^(0)   (choose starting direction)
iterate starting at n until convergence
  α^(n) = (r^(n))^T P^(-1) r^(n) / ((d^(n))^T A d^(n))   (Choose factor in d direction)
  φ^(n+1) = φ^(n) + α^(n) d^(n)   (Obtain new φ)
  r^(n+1) = r^(n) - α^(n) A d^(n)   (calculate new residual)
  β^(n) = (r^(n+1))^T P^(-1) r^(n+1) / ((r^(n))^T P^(-1) r^(n))   (Calculate coefficient to conjugate residual)
  d^(n+1) = P^(-1) r^(n+1) + β^(n+1) d^(n)   (obtain new conjugated search direction)

已经发展了许多预条件子，其复杂度范围广泛，从元素取自原矩阵 \(\mathbf{A}\) 对角元素的对角矩阵（Jacobi 预条件子）到使用不完全 Cholesky 分解的更复杂者。尽管如此，CG 方法在求解大型方程组时总是应与预条件子配合使用。

10.3.14 双共轭梯度法（BiCG）及其预条件形式（The Bi-conjugate Gradient Method (BiCG) and Preconditioned BICG）

第 8 章中介绍的对扩散方程进行离散化所得的系数矩阵以及一些其他方程（如第 15 章将介绍不可压缩压力或压力修正方程）的离散化矩阵是对称的，从而产生可以用前文讨论的 CG 方法求解的对称系统。然而，对 CFD 应用中出现的通用守恒方程进行离散化得到的矩阵 \(\mathbf{A}\) 是不对称的，所得的方程组也是不对称的。为能使用 CG 方法求解该系统，应将其变换为对称形式 [20]。一种做法是将式 10.1 改写为

\[ \begin{bmatrix} 0 & \mathbf{A} \\ \mathbf{A}^T & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol{\phi} \\ \tilde{\boldsymbol{\phi}} \end{bmatrix} = \begin{bmatrix} \mathbf{b} \\ 0 \end{bmatrix} \quad (10.121) \]

其中 \(\tilde{\boldsymbol{\phi}}\) 是为将原始不对称系统转换为可由 CG 方法求解的对称系统而引入的辅助变量。对该系统应用 CG 方法，会产生两组类 CG 向量序列：一组基于原系统（系数矩阵为 \(\mathbf{A}\)）的普通序列，用于计算 \(\boldsymbol{\phi}\)；另一组为基于系数矩阵 \(\mathbf{A}^T\) 的辅助系统（不需要求解）的"影子"序列，若需要可从中计算 \(\tilde{\boldsymbol{\phi}}\)。因为存在两组向量序列，该方法得名"双共轭梯度"（BiCG）。这里沿用 CG 方法中的同一术语。残差与搜索方向的普通序列分别记为 \(r\) 与 \(d\)，其影子形式分别记为 \(\tilde{r}\) 与 \(\tilde{d}\)。残差的双正交性通过如下构造得到保证：

\[ \tilde{r}^{(n)T} r^{(n)} = \tilde{r}^{(m)T} r^{(n)} = 0, \quad m \neq n \quad (10.122) \]

搜索方向的双共轭性通过要求

\[ \left(\tilde{d}^{(n)}\right)^T \mathbf{A} d^{(m)} = d^{(n)T}\mathbf{A}^T\tilde{d}^{(m)} = 0, \quad m \neq n \quad (10.123) \]

得到满足。

此外，残差与搜索方向的序列被构造成"普通形式"与"影子形式"相互正交。用数学语言写为

\[ \tilde{r}^{(n)T} d^{(m)} = r^{(n)T} \tilde{d}^{(m)} = 0, \quad m \neq n \quad (10.124) \]

该方法存在多种变体，且收敛不规则，有可能出现中断现象；下文给出的算法源自 Lanczos [21, 22]。

Lanczos 的 BiCG 算法可总结如下：

d^(0) = r^(0) = d̃^(0) = r̃^(0) = b - Aφ^(0)   (choose starting directions)
iterate starting at n until convergence
  α^(n) = (r^(n))^T r̃^(n) / ((d^(n))^T A d^(n))   (Choose factor in d direction)
  φ^(n+1) = φ^(n) + α^(n) d^(n)   (Obtain new φ)
  r^(n+1) = r^(n) - α^(n) A d^(n)   (calculate new r residual)
  r̃^(n+1) = r̃^(n) - α^(n) A^T d̃^(n)   (calculate new r̃ residual)
  β^(n) = (r^(n+1))^T r̃^(n+1) / ((r^(n))^T r̃^(n))   (Calculate coefficient to conjugate residual)
  d^(n+1) = r^(n+1) + β^(n+1) d^(n)   (obtain new search d direction)
  d̃^(n+1) = r̃^(n+1) + β^(n+1) d̃^(n)   (obtain new search d̃ direction)

BiCG 方法在每次迭代中需要与系数矩阵及其转置相乘，导致每次迭代的计算量约为 CG 方法的两倍。

预条件也可与 BiCG 方法配合使用。采用与预条件 CG 方法相同的术语，下文给出 Fletcher [23] 发展的该方法的一个鲁棒变体。

Fletcher 的 BiCG 预条件算法如下（\(\mathbf{P}\) 表示预条件矩阵）：

r^(0) = r̃^(0) = b - Aφ^(0),  d^(0) = P^(-1) r^(0),  d̃^(0) = P^T r̃^(0)   (choose starting directions)
iterate starting at n until convergence
  α^(n) = (r^(n))^T P^(-1) r̃^(n) / ((d^(n))^T A d^(n))   (Choose factor in d direction)
  φ^(n+1) = φ^(n) + α^(n) d^(n)   (Obtain new φ)
  r^(n+1) = r^(n) - α^(n) A d^(n)   (calculate new r residual)
  r̃^(n+1) = r̃^(n) - α^(n) A^T d̃^(n)   (calculate new r̃ residual)
  β^(n) = (r^(n+1))^T P^(-1) r̃^(n+1) / ((r^(n))^T P^(-1) r̃^(n))   (Calculate coefficient to conjugate residual)
  d^(n+1) = P^(-1) r^(n+1) + β^(n+1) d^(n)   (obtain new search d direction)
  d̃^(n+1) = P^T r̃^(n+1) + β^(n+1) d̃^(n)   (obtain new search d̃ direction)

更稳定、鲁棒的 BiCG 方法的其他变体已被报道，如 Sonneveld 的共轭梯度平方法（CGS）[24]、Van Der Vorst 的双共轭梯度稳定法（Bi-CGSTAB）[25] 以及广义最小残差法 GMRES [13, 26–29]。这些方法适用于求解 CFD 应用中产生的大型方程组，因为它们可应用于非对称矩阵以及结构与非结构网格。

10.4 多重网格方法（The Multigrid Approach）

迭代法的收敛速度随着代数系统规模的增大而显著下降，在中到大规模系统中甚至在初始误差被消除后仍能观察到收敛速度的下降。这已成为迭代求解器的一项严重局限。幸运的是，人们很快发现多重网格方法与迭代法的结合实际上可以弥补这一弱点。

多重网格方法的发展始于 Fedorenko [30]（几何多重网格）、Poussin [31]（代数多重网格）以及 Settari 和 Aziz [32] 的工作，并随着 Brandt [33] 的理论工作而获得更多关注。高频或振荡的误差容易被标准迭代求解器（Jacobi、Gauss-Seidel、ILU）消除，但这些方法难以消除光滑或低频误差分量 [34]。因此这些求解方法在多重网格的语境下被称为"光滑子"（smoothers）。误差频率的图示见图 10.4，其中绘制了一维问题不同误差频率模式的示意图。

图 10.4 中所示的误差模式从短波长的 \(\lambda_1\)（高频）到长波长的 \(\lambda_5\)（低频）不等，并在图的顶部集体绘出。一维域使用所示的一维网格进行离散化，各种模式分别绘制在同一网格上。可见，高频误差在一个单元上呈振荡状，易被迭代方法感知。随着误差频率降低或波长（\(\lambda\)）增大，误差在网格上变得越来越光滑——因为任何单元内只包含波长的一小部分。随着网格被进一步细化，方程数增多，这解释了为何收敛速度随系统规模增大而下降。

多重网格方法通过确保在任一网格层上施加光滑子所产生的低频误差被变换为粗网格层上的高频误差，从而提高迭代求解器的效率。通过使用粗网格层次（图 10.5），多重网格方法能够克服收敛速度的下降。

粗网格的构造通常可以采用两种方式：利用细网格的拓扑与几何（相当于在细网格上为每一粗层重新生成网格）；或直接对细网格单元进行聚类 [35–40]，后者又称代数多重网格法（AMG）。在 AMG 中，不直接需要或使用几何信息，聚合过程是纯代数的，每个粗层的方程通过聚合从细层重建。这种方法可用于为高度各向异性网格及/或方程系数大范围变化的问题构造高效且鲁棒的线性求解器。

无论采用哪种方法，都使用一种多重网格循环过程来引导对各网格层次的遍历。每一次从细网格到粗网格的遍历包含：(1) 限制过程，(2) 粗网格层方程组的建立或更新，(3) 若干次光滑子迭代。从粗网格到细网格的遍历包含：(1) 延拓过程，(2) 细网格层场值的修正，(3) 对限制阶段所建方程的若干次光滑子迭代。所需各步将在下文详细说明。

10.4.1 单元聚合/粗化（Element Agglomeration/Coarsening）

求解过程的第一步是通过聚合/粗化算法生成粗/细网格层。为此可以采用三种不同的方法。第一种方法是先生成粗网格，再通过细化得到细网格层 [41, 42]。这便于定义粗—细网格关系，在自适应网格设置中很有吸引力 [41–43]。然而其主要缺点在于细网格分布依赖于粗网格。第二种方法使用非嵌套网格 [44]，使得网格层间的信息传递代价高昂。此外，这两种方法都不能很好地解析复杂域。在第三种方法（此处推荐的方法）中，过程从生成最细网格开始。然后通过细网格单元的聚合 [45, 46]（如图 10.6 所示）建立粗网格层，聚合过程基于单元几何或基于邻居单元系数须满足的某种准则。下面的讨论与第三种方法相关。

粗网格层通过聚合算法融合细网格单元而生成。对每一粗网格层，算法反复施加，直至细层的所有网格单元都与粗网格单元相关联。

在这一启发式聚合过程中，逐一访问细网格点。选取一个单元作为种子单元，将满足预设准则的一定数量的相邻单元与之融合形成粗单元。一个粗单元所能融合的最大细单元数是预先设定的。若所选种子单元未能形成粗单元，则将其加入邻居中"最不饱满"的粗单元中。

一种高效的聚合算法是 Mavriplis 发展的方向性聚合（DA）算法 [47]。在 DA 中，聚合从某个种子单元出发，根据相邻细网格单元几何连接强度将它们与种子单元合并。该过程只需在求解开始时执行一次。

10.4.2 限制步骤与粗层系数（The Restriction Step and Coarse Level Coefficients）

求解从细网格层开始。经过几次迭代后，误差被传递或限制到粗一层，并在该层求解。然后在该层经过几次迭代后，误差再次被限制到更高层，过程序列重复，直至达到最高或最粗的网格层。设 \((k)\) 表示已通过求解如下校正形式方程组得到解的某一层：

\[ \mathbf{A}^{(k)}e^{(k)} = r^{(k)} \quad (10.125) \]

下一粗一层为 \((k + 1)\)，误差将被限制到该层。设 \(G_I\) 表示细网格层 \((k)\) 上聚合成粗网格层 \((k + 1)\) 的单元 \(I\) 的细网格单元 \(i\) 的集合。那么粗网格层上需要求解的方程组为

\[ \mathbf{A}^{(k+1)}e^{(k+1)} = r^{(k+1)} \quad (10.126) \]

式 10.126 右端的残差按下式计算：

\[ r^{(k+1)} = \mathbf{I}_k^{k+1} r^{(k)} \quad (10.127) \]

其中 \(\mathbf{I}_k^{k+1}\) 是从细网格到粗网格的限制算子（即插值矩阵），由聚合过程定义。在 AMG 中，限制算子以线性方式定义，对细网格残差作求和：

\[ r_I^{(k+1)} = \sum_{i \in G_I} r_i^{(k)} \quad (10.128) \]

此外，粗单元的系数由构成它的细单元的相应系数相加而构造。回顾离散化后线性方程的形式为

\[ a_C\phi_C + \sum_{F = NB(C)} a_F\phi_F = b_C \quad (10.129) \]

为当前目的，对细网格层改写为更合适的形式：

\[ a_i^{(k)}\phi_i^{(k)} + \sum_{j = NB(i)} a_{ij}^{(k)}\phi_j^{(k)} = b_i^{(k)} \quad (10.130) \]

其中 \(NB(i)\) 表示单元 \(i\) 的邻居。初始时式 10.130 并不满足，从而产生如下残差：

\[ r_i^{(k)} = b_i^{(k)} - \left(a_i^{(k)}\phi_i^{(k)} + \sum_{j = NB(i)} a_{ij}^{(k)}\phi_j^{(k)}\right) \quad (10.131) \]

记 \(\phi_I^{(k+1)}\) 为细网格单元 \(i\) 之父单元粗网格单元 \(I\) 上的解，则来自粗网格的细网格修正可写为

\[ \phi_i'^{(k)} = \phi_I^{(k+1)} - \phi_i^{(k)} \quad (10.132) \]

期望该修正能使粗网格单元 \(I\) 上的残差为零。这些新残差记为 \(\tilde{r}_i^{(k)}\)，其计算为

\[ \tilde{r}_i^{(k)} = b_i^{(k)} - \left[a_i^{(k)}\left(\phi_i^{(k)} + \phi_i'^{(k)}\right) + \sum_{j = NB(i)} a_{ij}^{(k)}\left(\phi_j^{(k)} + \phi_j'^{(k)}\right)\right] \quad (10.133) \]

或等价地写为

\[ \tilde{r}_i^{(k)} = \underbrace{b_i^{(k)} - \left(a_i^{(k)}\phi_i^{(k)} + \sum_{j = NB(i)} a_{ij}^{(k)}\phi_j^{(k)}\right)}_{r_i^{(k)}} - \left[a_i^{(k)}\phi_i'^{(k)} + \sum_{j = NB(i)} a_{ij}^{(k)}\phi_j'^{(k)}\right] \quad (10.134) \]

即

\[ \tilde{r}_i^{(k)} = r_i^{(k)} - \left[a_i^{(k)}\phi_i'^{(k)} + \sum_{j = NB(i)} a_{ij}^{(k)}\phi_j'^{(k)}\right] \]

。

强制 \(I\) 中的残差之和为零，即

\[ \sum_{i \in G_I} \tilde{r}_i^{(k)} = 0 \quad (10.135) \]

并将式 10.134 代入式 10.135，得到

\[ 0 = \sum_{i \in G_I} r_i^{(k)} - \sum_{i \in G_I}\left[a_i^{(k)}\phi_i'^{(k)} + \sum_{j = NB(i)} a_{ij}^{(k)}\phi_j'^{(k)}\right] \quad (10.136) \]

利用粗网格编号重写式 10.136，粗网格修正方程变为

\[ a_I^{(k+1)}\phi_I'^{(k+1)} + \sum_{J = NB(I)} a_{IJ}^{(k+1)}\phi_J'^{(k+1)} = r_I^{(k+1)} \quad (10.137) \]

其中 \(a_I^{(k+1)}, a_{IJ}^{(k+1)}, r_I^{(k+1)}\) 由细网格系数直接推导得到：

\[ a_I^{(k+1)} = \sum_{i \in G_I} a_i^{(k)} + \sum_{i \in G_I}\sum_{j \in G_I} a_{ij}^{(k)} \]

\[ a_{IJ}^{(k+1)} = \sum_{\substack{i \in G_I \\ j \notin G_I, j \in NB(I)}} a_{ij}^{(k)} \quad (10.138) \]

\[ r_I^{(k+1)} = \sum_{i \in G_I} r_i^{(k)} \]

这一过程如图 10.7 所示。

10.4.3 延拓步骤与细网格层修正（The Prolongation Step and Fine Grid Level Corrections）

延拓算子用于将修正从粗网格层传递到细网格层。可采用多种方式。如图 10.8 所示的一种方式是零阶延拓算子，其使细网格上的误差值相同——即粗网格单元上的误差将被其所有子单元在细网格层上继承。

修正基本上由粗网格上的方程组的解获得。向细网格层的插值或延拓记为

\[ e^{(k)} = \mathbf{I}_{k+1}^k e^{(k+1)} \quad (10.139) \]

其中 \(\mathbf{I}_{k+1}^k\) 是从粗网格到细网格的插值矩阵。最终细网格的解被修正为

\[ \boldsymbol{\phi}^{(k)} \leftarrow \boldsymbol{\phi}^{(k)} + e^{(k)} \quad (10.140) \]

所用网格层的数量取决于网格的规模。对更多的网格层，过程与图 10.8 中所勾勒的相同。

10.4.4 遍历策略与代数多重网格循环（Traversal Strategies and Algebraic Multigrid Cycles）

遍历策略是指在求解过程中访问粗网格的方式，又称多重网格循环 [48]。AMG 方法中常用的循环有 V 循环、W 循环与 F 循环 [35, 36, 49, 50]，如图 10.9 所示。

最简单的 AMG 循环如图 10.9a 所示，为 V 循环 [49, 50]，即仅访问每个网格层一次。通常的做法是在限制阶段执行几次迭代扫描，然后将残差注入粗网格直至达到最粗层。对非常 stiff 的系统，V 循环可能不足以加速收敛，因此需要在粗层上执行更多迭代。W 循环基于在每个被访问的粗网格层上施加较小的 V 循环。W 循环（图 10.9b）由嵌套的粗—细网格层扫描组成，其复杂性随 AMG 层数的增加而增加。F 循环是 W 循环的一种变体，可被视为 W 循环对半拆开，如图 10.9c 所示。F 循环所需的粗层扫描少于 W 循环但多于 V 循环。因此它介于 V 与 W 循环策略之间。

10.5.1 uFVM（uFVM）

uFVM 实现了两种线性代数求解器：逐次超松弛法（SOR）以及 ILU(0) 法。这些方法的实现遵循前述步骤。SOR 位于文件 cfdSORSolver.m 中，ILU(0) 实现位于 cfdILUSolver.m。

10.5.2 OpenFOAM®（OpenFOAM®）

OpenFOAM® 中迭代线性代数求解器的组织结构 [51] 遵循通常的做法。它首先定义基类，再由基类派生出每种类型的代数矩阵求解器。这些代数求解器分为三大类：求解器（solvers）、预条件子（preconditioners）以及光滑子（smoothers）。光滑子与预条件子通过将光滑子与不动点关系联系起来、并将其嵌入到预条件子框架中加以区分。回顾式 10.81，预条件子类实现乘积 \(\mathbf{P}^{-1}r\)，而光滑子类推进解。此外，求解器类汇集与共轭梯度及多重网格算法实现相关的必要信息。

线性代数求解器的源代码位于 …/src/OpenFOAM/matrices/lduMatrix/ 路径下的 lduMatrix 目录中，分为以下三个子目录：solvers、preconditioners、smoothers。每个子目录的名称反映其功能。solvers 目录包含 OpenFOAM® 实现的迭代求解器主体代码，包括：diagonalSolver（用于对称与非对称问题的对角求解器）；GAMG（几何聚合代数多重网格求解器，文档中也称为"广义几何—代数多重网格"）；ICC（不完全 Cholesky 预条件共轭梯度求解器）；PBiCG（用于非对称矩阵的预条件双共轭梯度求解器）；PCG（用于对称矩阵的预条件共轭梯度求解器）；smoothSolver（基于预条件子的对称与非对称矩阵光滑迭代求解器）。

preconditioners 目录包含各种对角 ILU 实现，包括：diagonalPreconditioner（对角预条件子）；DICPreconditioner、DILUPreconditioner（分别用于对称与非对称矩阵的对角不完全 Cholesky 预条件子）；FDICPreconditioner（DICPreconditioner 的更快版本，用于对称矩阵的对角不完全 Cholesky 预条件子，其中预条件对角元素的倒数以及除以对角后的上系数被预先计算并存储）；GAMGPreconditioner（几何聚合代数多重网格预条件子，使用多重网格循环作为预条件子以执行式 10.81 的第二部分）；noPreconditioner（用于对称与非对称矩阵的空预条件子）。

smoothers 目录包含以下内容：DIC、DILU（用于对称与非对称矩阵的对角不完全 Cholesky 光滑子）；DICGaussSeidel、DILUGaussSeidel（用于对称与非对称矩阵的 DIC/DILU 与 Gauss-Seidel 复合光滑子，其中 DIC/DILU 平滑之后接 Gauss-Seidel 以确保任何由 DIC/DILU 扫描产生的"尖峰"被抹平）；DILU（用于非对称矩阵的对角不完全 LU 光滑子）；GaussSeidel（用于对称与非对称矩阵的 Gauss-Seidel 方法）。

此外，OpenFOAM® 在 lduMatrix 类内定义了三个额外的基类，分别封装上述三个对应类别。因此 lduMatrix.H 文件读为（Listing 10.1）

class lduMatrix
{
    // private data

        //- LDU mesh reference
        const lduMesh& lduMesh_;

        //- Coefficients (not including interfaces)
        scalarField *lowerPtr_, *diagPtr_, *upperPtr_;

public:

    //- Abstract base-class for lduMatrix solvers
    class solver
    {
    protected:
        ...
    };

    class smoother
    {
    protected:
        ...
    };

    class preconditioner
    {
    protected:
        ...
    };
};

（Listing 10.1 由 lduMatrix 类定义的三个基类：solver、smoother 与 preconditioner。）

因此每个 smoother、solver 与 preconditioner 都必须从这三个基类派生。例如，DILU 预条件子按 Listing 10.2 所示声明：

class DILUPreconditioner
:
    public lduMatrix::preconditioner
{
    ...
};

（Listing 10.2 用于声明 DILU 预条件子的语法。）

而共轭梯度求解器则按 Listing 10.3 所示声明：

class PCG
:
    public lduMatrix::solver
{
    ...
};

（Listing 10.3 用于声明 PCG 求解器的语法。）

无论哪种情况，预条件子或求解器显然都是从 lduMatrix 类下定义的基类派生的。

在阐明 OpenFOAM® 中代数求解器的基本概念与组织结构后，下文给出一个详细说明预条件 CG 方法实现细节的示例。文件位于目录 $FOAM_SRC/OpenFOAM/matrices/lduMatrix/solvers/PCG。

从 lduMatrix::solver 类派生的类（如 Listing 10.3 所示）使用 Listing 10.4 中的脚本定义主成员函数 solve。

    // Member Functions

        //- Solve the matrix with this solver
        virtual solverPerformance solve
        (
            scalarField& psi,
            const scalarField& source,
            const direction cmpt=0
        ) const;

（Listing 10.4 用于定义成员函数 solve 的脚本。）

函数 solve 在文件 PCG.C 中实现所选线性代数求解器的求解算法。回顾预条件共轭梯度算法，其事件序列如下：

1. 计算 \(r^{(0)} = \mathbf{b} - \mathbf{A}\boldsymbol{\phi}^{(0)}\)
2. 计算 \(\mathbf{d}^{(0)} = \mathbf{P}^{-1}r^{(0)}\)
3. 计算

\[ \alpha^{(n)} = \dfrac{r^{(n)T}\mathbf{P}^{-1}r^{(n)}}{\mathbf{d}^{(n)T}\mathbf{A}\mathbf{d}^{(n)}} \]

1. 计算

\[ \boldsymbol{\phi}^{(n+1)} = \boldsymbol{\phi}^{(n)} + \alpha^{(n)}\mathbf{d}^{(n)} \]

1. 计算 \(r^{(n+1)} = r^{(n)} - \alpha^{(n)}\mathbf{A}\mathbf{d}^{(n)}\)
2. 若解已收敛则停止
3. 计算

\[ \beta^{(n+1)} = \dfrac{r^{(n+1)T}\mathbf{P}^{-1}r^{(n+1)}}{r^{(n)T}\mathbf{P}^{-1}r^{(n)}} \]

1. 计算 \(\mathbf{d}^{(n+1)} = \mathbf{P}^{-1}r^{(n+1)} + \beta^{(n+1)}\mathbf{d}^{(n)}\)
2. 跳至步骤 3

该算法在 PCG.C 中按相同过程直接实现。

在步骤 1 中（Listing 10.5），计算残差并存入变量 rA，而变量 wA 存储矩阵—解乘积 \(\mathbf{A}\boldsymbol{\phi}^{(0)}\)：

// --- Calculate A.psi
matrix_.Amul(wA, psi, interfaceBouCoeffs_, interfaces_, cmpt);
// --- Calculate initial residual field
scalarField rA(source - wA);

（Listing 10.5 用于计算残差的脚本。）

步骤 2 涉及预条件。因此，首先用对象 preconPtr 定义所使用的预条件子类型（Listing 10.6）：

    // --- Select and construct the preconditioner
    autoPtr<lduMatrix::preconditioner> preconPtr =
    lduMatrix::preconditioner::New
    (
        *this,
        controlDict_
    );

（Listing 10.6 定义所使用的预条件子类型的语法。）

所使用的构造函数是基于基类与 New 构造函数的通用构造函数。预条件子类型在运行时从字典中选取。然后，预条件操作 \(\mathbf{P}^{-1}r^{(n)}\) 被施加到残差 rA（依据算法步骤 2 中的方程）。结果被存入同一变量 wA 以减少内存使用，然后用 gSumProd 函数对 wA 与 rA 两个向量作内积即可求得 \(r^{(n)T}\mathbf{P}^{-1}r^{(n)}\)，如 Listing 10.7 所示。此外，前一次迭代 \((n-1)\) 的 wArA 旧值被存入变量 wArAold。

wArAold = wArA;

// --- Precondition residual
preconPtr->precondition(wA, rA, cmpt);

// --- Update search directions:
wArA = gSumProd(wA, rA, matrix().mesh().comm());

（Listing 10.7 用于计算 \(r^{(n)T}\mathbf{P}^{-1}r^{(n)}\) 并存储其旧值的语法。）

按预条件共轭梯度算法，步骤 3、4、5 如 Listing 10.8 所示执行：

// --- Update preconditioned residual
matrix_.Amul(wA, pA, interfaceBouCoeffs_, interfaces_, cmpt);

scalar wApA = gSumProd(wA, pA, matrix().mesh().comm());

// --- Update solution and residual:

scalar alpha = wArA/wApA;

for (register label cell=0; cell<nCells; cell++)
{
    psiPtr[cell] += alpha*pAPtr[cell];
    rAPtr[cell] -= alpha*wAPtr[cell];
}

（Listing 10.8 用于计算 \(\alpha\) 并更新因变量与残差值的脚本。）

现在变量 wA 存储乘积 \(\mathbf{A}\mathbf{d}^{(n)}\)，而 pA 表示向量 \(\mathbf{d}^{(n)}\)。再次使用 gSumProd 函数计算乘积 \(\mathbf{d}^{(n)T}\mathbf{A}\mathbf{d}^{(n)}\) 并存入变量 wApA。一旦 \(\alpha\) 被求出，步骤 4 与 5 中残差与解的更新在 for 循环中执行，Listing 10.8 中的变量 psiPtr 与 rAPtr 分别代表 \(\boldsymbol{\phi}\) 与 \(r\)。

通过执行 Listing 10.9 中的脚本完成算法中的步骤 7 与 8。

scalar beta = wArA/wArAold;

for (register label cell=0; cell<nCells; cell++)
{
    pAPtr[cell] = wAPtr[cell] + beta*pAPtr[cell];
}

（Listing 10.9 用于计算 \(\beta\) 与 \(\mathbf{d}\) 的新值以便在下次迭代中使用的脚本。）

在实际算例中，线性求解器的定义在 system 目录下的 fvSolution 文件中的 solvers {} 块内完成，如 Listing 10.10 所示。

solvers
{
    T
    {
        solver               PCG;
        preconditioner       DIC;
        tolerance            1e-06;
        relTol               0;
    }
}

（Listing 10.10 线性求解器的定义。）

Listing 10.10 中各项的含义为：solver 定义求解器类型（此处 PCG 为对称矩阵的预条件共轭梯度），可用选项包括 PCG（仅适用于对称矩阵的预条件共轭梯度）、PBiCG（仅适用于非对称矩阵的预条件双共轭梯度）、smoothSolver（仅用作光滑子以减少残差的求解器）、GAMG（广义几何—代数多重网格，应在大网格上用于压力方程）。preconditioner 定义所使用的预条件子类型，可选项包括 DIC（对称矩阵的对角不完全 Cholesky 预条件子）、DILU（非对称矩阵的对角不完全 LU 预条件子）。tolerance 为线性求解器停止迭代时所允许的最大绝对残差值。relTol 为初始残差与当前残差之比，线性求解器停止迭代时该比值应满足要求。其他示例如 Listing 10.11 所示。

solvers
{
    T
    {
        solver               PBiCG;
        preconditioner       DILU;
        tolerance            1e-06;
        relTol               0;
    }
}
solvers
{
    T
    {
        solver                smoothSolver;
        smoother              GaussSeidel;
        tolerance             1e-8;
        relTol                0.1;
        nSweeps               1;
    }
    T
    {
        solver           GAMG;
        tolerance        1e-7;
        relTol           0.01;
        smoother         GaussSeidel;
        nPreSweeps       0;
        nPostSweeps      2;
        cacheAgglomeration on;
        agglomerator     faceAreaPair;
        nCellsInCoarsestLevel 10;
        mergeLevels      1;
    }
}

（Listing 10.11 线性求解器定义的示例。）

10.6 小结（Closure）

本章介绍了直接法与迭代法两大类用于求解代数方程组的方法，并在每一类中描述了若干具体方法。代数多重网格技术也作了讨论。下一章将继续守恒方程的离散化，并详细描述对流项的离散。

10.7 练习（Exercises）

练习 1–14 包含若干 LU 分解、Gauss-Seidel/Jacobi 迭代（含与不含预条件）、ILU(0)、DILU 与预条件共轭梯度方法的练习题，以及针对 OpenFOAM® 中线性求解器、光滑子与多重网格 V 循环实现的查找任务；其中前几题要求对给定矩阵分别执行 LU 分解（包括前向与回代两步）、Gauss-Seidel 与 Jacobi 迭代若干次并比较误差、构造 Jacobi 与 Gauss-Seidel 预条件矩阵、按 ILU(0)/DILU 完成分解并计算误差矩阵 \(R = M - P\)、对示例矩阵运行 ILU(0) 与 DILU 三次迭代（零初值）、比较 Gauss-Seidel/Jacobi 与 ILU(0)/DILU 三次迭代误差、对称系统执行预条件共轭梯度两次迭代；后几题则要求列举 OpenFOAM® 提供的所有线性求解器、列在选择 smoothSolver 后可用的所有 smoother、查找 PBiCG 实现并与 Lanczos BiCG 算法比较、查找多重网格 V 循环实现（位于 $FOAM_SRC/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C）并与理论 V 循环比较、以及验证 Pommerell 对角 ILU(0) 的实现。

本章个人批注

本章是 Moukalled FVM 一书中工程味最重的章节之一，因为它把"求解线性方程组"这一抽象的数值代数话题与 CFD 的实际语境彻底结合了起来。从结构上看，本章有一个明显的递进：从 Gauss 消去到 LU 分解，再到稀疏结构（带状）的专门算法（TDMA、PDMA），接着跳到一般稀疏情形下的迭代法（Jacobi、Gauss-Seidel），最后到加速收敛的两类思路——预条件子（ILU、DILU）以及 Krylov 子空间方法（CG、BiCG 及其变体），并以多重网格作为对迭代法局限性的补充。这一安排本身就告诉读者：直接法只在结构网格带状情形下"几乎可行"，一旦遇到非结构网格或非线性耦合问题，迭代法成为唯一选项。

让我印象深刻的是第 10.3.5 节对 ILU 的论述。作者特意强调"incomplete"——完全 LU 分解等价于直接求解法且代价高昂，而不完全分解则保留稀疏结构但只近似 \(\mathbf{A}\)，于是不得不嵌入到一个外层迭代里。这种"分解 \(\to\) 近似 \(\to\) 迭代"的层次结构，正是 CFD 中代数求解器设计的核心思想。DILU（Pommerell 1992）作为 ILU 的简化版，只修改对角线而保持上/下三角部分不变——这种"只动对角、不动非对角"的简化思想在工程实现中非常重要，因为它把额外的存储开销压到一对角线，这正是 OpenFOAM 中 DILUPreconditioner 类的物理对应。

CG 与 BiCG 的对照则揭示了 CFD 中一类系统性的张力：扩散方程（第 8 章）以及压力/压力修正方程（第 15 章）所产生的矩阵是对称正定的，因此可以走 CG 这条"优雅、高效、可理论分析"的路线；但通用守恒方程（包括对流占优的情形）离散化得到的矩阵是非对称的，CG 不再适用。BiCG 通过引入影子系统把问题转化为对称形式，但代价是每次迭代要做两次矩阵—向量乘法（一次 \(\mathbf{A}\)、一次 \(\mathbf{A}^T\)）。作者进一步提到 CGS、Bi-CGSTAB、GMRES 等更稳定的变体——这暗示了一个工程现实：BiCG 的"不规则收敛"和"中断可能性"让它在生产代码中很少被直接采用；Bi-CGSTAB 与 GMRES 才是常见的稳健选择。但本章的章节安排（CGS、Bi-CGSTAB、GMRES 只在 BiCG 章节末尾一笔带过）也意味着 Moukalled 这本书是有取舍的——它讲了算法的"骨架"，但具体的实现细节和鲁棒性变体则留给 OpenFOAM/uFVM 的工程实践。

多重网格那一节也值得单独说几句。本章论述的多重网格并不是"几何多重网格"（GMG）的完整介绍，而是把多重网格作为"弥补迭代法低频误差的缺陷"的工具来讨论。具体做法是"细网格 → 粗网格"通过限制（restriction）算子传递残差，"粗网格 → 细网格"通过延拓（prolongation）算子传递修正。作者列出了 V、W、F 三种循环并说明了它们各自的代价—加速权衡。值得注意的是，作者对代数多重网格（AMG）特别强调——这是因为 AMG 不依赖几何信息，对于复杂域或大变形网格问题特别有用。这种"先讲几何理解、再讲代数实现"的写法是有意识的方法论选择。

第 10.5 节（Computational Pointers）从理论跳到 OpenFOAM/uFVM 的具体实现，是非常实用的桥梁。从这节中我学到的最重要一点是 OpenFOAM 把所有线性求解器都封装在 lduMatrix 类下的三个抽象基类（solver、smoother、preconditioner）里。这个设计反映了上述的理论分类：CG、BiCG 属于 solver；DIC、DILU 属于 preconditioner；Gauss-Seidel、DICGaussSeidel 属于 smoother。这种"基类—派生类"的 C++ 实现模式让 OpenFOAM 能在运行时通过字典（fvSolution 中的 solvers {} 块）切换不同的求解器与预条件子，无需重新编译——这是 OpenFOAM 高度可配置的根本原因。PCG.C 中的代码片段与本章第 10.3.13 节的算法步骤一一对应：步骤 1 对应 Listing 10.5 计算残差、步骤 2 对应 Listing 10.6 构造 preconditioner、步骤 3–5 对应 Listing 10.8 更新解与残差、步骤 7–8 对应 Listing 10.9 计算 \(\beta\) 与新的搜索方向。这种"理论算法 → 实际 C++ 代码"的对应关系，对我今后阅读 OpenFOAM 求解器源码会有直接的指导作用。

至于习题部分，本章的 14 道习题覆盖面非常广——从 LU 分解的手工演算（练习 1 给出的 5×5 矩阵对初学者是个不小的考验），到迭代法的收敛行为比较（练习 2、3 涉及 Jacobi、Gauss-Seidel 与其预条件版本的对比），到不完全分解的误差矩阵 \(R = M - P\) 的构造（练习 5、7），再到对称系统的预条件 CG（练习 10），最后到 OpenFOAM 实现的源代码考古（练习 11–14）。这种"从手工推导到代码考古"的梯度设计说明 Moukalled 这本书既重视理论推导也重视工程实现。

与上下章的衔接（一段话）

第 10 章在全书结构中处于一个明显的转折点：前 9 章（Ch.1–Ch.9）逐步构建了 FVM 的几何与离散化基础——从守恒方程的一般形式，到网格与控制体，再到对流项与扩散项的离散模板，最后到梯度计算——而从本章开始，后续章节将围绕"如何高效求解所得到的代数方程组"这一工程核心展开。具体来说，第 11 章将开启对流项离散格式的系统讨论（高阶格式 TVD 性质等），第 12 章讨论压力—速度耦合（SIMPLE、SIMPLEC、PISO），而这两类讨论最终都会回到本章所介绍的代数求解器：非结构网格下的对流项离散会产生高度不对称的矩阵，需要 BiCGSTAB 或 GMRES 类求解器；压力—速度耦合则会产生大型稀疏方程组，需要多重网格加速。本章放在这个位置正是因为它是连接"离散化"（前 9 章）与"求解实践"（后 11 章）之间的桥梁——前 9 章的所有努力最终都要落到本章所讨论的代数求解器上，而后 11 章的所有算法创新（如 SIMPLE 系列、压力修正方程）也都依赖本章所介绍的预条件子、共轭梯度类方法以及多重网格技术来获得工程上可接受的计算效率。