FVM in CFD 学习笔记_第10章_求解代数方程组系统

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 10 Solving the System of Algebraic Equations

离散过程将会生成线性方程组系统，$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$，其中未知量$\mathbf{\varphi }$位于网格单元形心上，是待求量。该系统中，未知变量的系数构成了矩阵$\mathbf{A}$，其由线性化过程和网格几何量导出，而向量$\mathbf{b}$则包含了源项、常数项、边界条件和非线性分量。

求解线性方程组系统的技术通常被分为直接方法和迭代方法，每类又有很多具体方法。由于流动问题是高度非线性的，从其线性化过程得到的系数通常是依赖于解的，因此，在每步的迭代过程中并不需要得到非常精确的解，所以直接解法在CFD的应用中很少使用。最常用的还是迭代解法，因为它们更适合求解该类问题，其在每步迭代中所需内存更小，计算消耗更小。

本章首先讲解在结构和非结构网格上的一些直接解法（Gauss消元、LU分解、三对角和五对角矩阵算法），以便为在CFD应用中更加广泛使用的迭代方法提供基础。然后回顾一些基本的迭代解法（含预处理和不含预处理）的特性和局限性，包括Jacobi、Gauss-Seidel、不完全LU分解、以及共轭梯度（CG）方法。最后，简要讲讲多重网格方法，它通常是和迭代方法联合使用，以克服这些迭代方法的局限性。

1 引言

线性求解器求解的是如下形式的代数方程组
$$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$$
其中$\mathbf{A}$是单元的系数矩阵$a_{ij}$，$\mathbf{\varphi }$是未知变量$\varphi$的向量，而$\mathbf{b}$则是源项$b_i$的矢量。使用矩阵编号，该方程的展开形式为
$$\left[\begin{array}{ccccc}{a}_{11}& a_{12}& ...& a_{1N-1}& a_{1N}\\ a_{21}& a_{22}& ...& a_{2N-1}& a_{2N}\\ ...& ...& ...& ...& ...\\ a_{N1}& a_{N2}& ...& a_{NN-1}& a_{NN}\end{array}\right]\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ ...\\ {\varphi }_{N-1}\\ {\varphi }_N\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_{N-1}\\ b_N\end{array}\right]$$
一般来说呢，矩阵$\mathbf{A}$中的每一行代表计算域1个单元上定义的方程，非0系数是与该单元紧邻的邻居单元的，系数$a_{ij}$反映的是存储在$i$单元控制体形心的${\varphi }_i$与其邻近单元形心的${\varphi }_j$的关系紧密程度。由于一个单元只和少数几个单元相邻，且其邻接关系是与离散区域上单元的剖分编码方式相关的，所以大部分的系数实际上是0，这就使得矩阵$\mathbf{A}$总是个稀疏矩阵（即，非零元素的个数很少，在整体矩阵中占比很小，绝大多数的元素都是0）。更进一步，如果用的是结构化的网格，那么矩阵$\mathbf{A}$将会是带状矩阵，即，其仅含几条非零的对角线元素。因此，求解该类系统时，可针对其特性选用对应的高效方法。

如上所述，求解代数方程组系统的方法可以大体上分为两大类，即直接和迭代方法。在直接方法中，矩阵$\mathbf{A}$求逆，解$\varphi$通过$\varphi ={\mathbf{A}}^{-1}\mathbf{b}$来求出。当矩阵$\mathbf{A}$非常大的时候，应用直接线性解法来求解CFD问题是不切实际的，因为这些CFD问题通常包含非线性系统方程， 即它们的系数是与解相关联的，那么使用迭代方法是切合实际的。

另一方面，在迭代代数解法中，求解算法将重复多次直至达到预期的收敛水平，而不需要在每次迭代的过程中都获得完全收敛的解。

接下来，首先展示一些应用于结构和非结构网格的直接解法，紧接着是结构网格上带状矩阵的高效解法，该章的重点是，迭代线性代数求解器，它们已经被广为证明是在FVM中最有效的和最经济的方法，而且已经被包含进几乎所有有限体积代码的线性求解器中。

2 直接或Gauss消元方法

尽管直接方法并非求解线性代数方程组中稀疏系统的高效方法（它们的计算消耗实在是大的离谱），然而对于它们的讨论可为后续高效的迭代方法的引入铺平道路，所以还是讲一下为好。最简单的直接方法是Gauss消元方法，将被首先介绍。使用Gauss消元法时，将系统转化成等效的上三角系统，这激发了下三角上三角分解方法（LU分解方法），将随后介绍。该方法中，矩阵$\mathbf{A}$将被分解成两个矩阵$\mathbf{L}$和$\mathbf{U}$的乘积，其中$\mathbf{L}$是下三角矩阵，而$\mathbf{U}$是上三角矩阵。该过程也被称为LU因子分解。此外，还将讨论应用于从结构网格所推出的带状矩阵$\mathbf{A}$的直接方法。

2.1 Gauss消元

假设有如下2变量${\varphi }_1$和${\varphi }_2$的线性方程组
$$\begin{gathered} a_{11}{\varphi }_1+a_{12}{\varphi }_2=b_1 \\ {a}_{21}{\varphi }_1+a_{22}{\varphi }_2=b_2 \end{gathered}$$
为了消去${\varphi }_1$，将第1个式子乘上$a_{21}/a_{11}$，并减去第2个式子，可得
$$\left(a_{22}-\frac{a_{21}}{a_{11}}{a}_{12}\right){\varphi }_2=b_2-\frac{a_{21}}{a_{11}}{b}_1$$
从而可以直接求得${\varphi }_2$
$${\varphi }_2=\frac{b_2-\frac{a_{21}}{a_{11}}{b}_1}{a_{22}-\frac{a_{21}}{a_{11}}{a}_{12}}$$
再把这个求得的${\varphi }_2$代入到方程最初的第1个式子里，可以求得${\varphi }_1$的值
$${\varphi }_1=\frac{b_1}{a_{11}}-\frac{a_{12}}{a_{11}}\frac{b_2-\frac{a_{21}}{a_{11}}{b}_1}{a_{22}-\frac{a_{21}}{a_{11}}{a}_{12}}$$
上述过程实际上分成了2步，第1步是消去1个未知量，得到仅剩1个未知量的方程，求解该方程，可得该未知量的值；第2步是将该未知量回代到原方程中，求出剩余未知量的值。即，消元-回代，同样的思路也可用于N个变量的线性代数方程组。

2.2&2.3 前向消元Forward Elimination

直接给出算法流程

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
	}
}

2.4&2.5 反向回代Backward Substitution

直接给出算法流程

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
}

为防止被0除，消元的时候可选用最大行主元，并互换两行的位置。实际上，Gauss消元回代的方法计算量是非常大的，对于一个N个方程的线性系统，其计算量与$N^3/3$成正比，其中回代过程仅仅占了$N^2/2$，这样高的计算消耗迫使人们针对稀疏矩阵系统寻求更加高效的解法。

2.6 LU分解

求解线性代数方程组的另一个直接解法是LU或PLU（P即前面所讲的行主元选择过程），咱这里就讲讲LU就好了。LU实际上是Gauss方法的变种，LU方法比Gauss消元方法的优势在于，一旦LU分解执行之后，对于右端项$\mathbf{b}$不同的线性系统就想求解多少次就求解多少次，而不需要再做消元处理了（相当于把第1次的消元处理做了LU分解），然而在Gauss消元方法中，消元是始终需要进行的。

基于前面消元处理过程，是将原本的矩阵$\mathbf{A}$转化为上三角矩阵，即
$$\left[\begin{array}{cccccc}{u}_{11}& u_{12}& u_{13}& ...& u_{1N-1}& u_{1N}\\ 0& u_{22}& u_{23}& ...& u_{2N-1}& u_{2N}\\ 0& 0& u_{33}& ...& u_{3N-1}& u_{3N}\\ ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& 0& u_{NN}\end{array}\right]\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ {\varphi }_3\\ ...\\ {\varphi }_N\end{array}\right]=\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\\ ...\\ c_N\end{array}\right]$$
使用缩写形式，即
$$\mathbf{U}\mathbf{\varphi }-\mathbf{c}=0$$
令$\mathbf{L}$为单位下三角矩阵（对角线上的元素是1），即
$$\mathbf{L}=\left[\begin{array}{cccccc}1& 0& 0& ...& 0& 0\\ l_{21}& 1& 0& ...& 0& 0\\ l_{31}& l_{32}& 1& ...& 0& 0\\ ...& ...& ...& ...& ...& ...\\ l_{N1}& l_{N2}& l_{N3}& ...& l_{NN-1}& 1\end{array}\right]$$
将方程$\mathbf{U}\mathbf{\varphi }-\mathbf{c}=0$左乘以$\mathbf{L}$，则变回了最初的方程
$$\mathbf{L}(\mathbf{U}\mathbf{\varphi }-\mathbf{c})=\mathbf{L}\mathbf{U}\mathbf{\varphi }-\mathbf{L}\mathbf{c}=\mathbf{A}\mathbf{\varphi }-\mathbf{b}$$
这样，变得到
$$\mathbf{L}\mathbf{U}=\mathbf{A}$$
和
$$\mathbf{L}\mathbf{c}=\mathbf{b}$$
即，矩阵$\mathbf{A}$写成了下三角和上三角矩阵的乘积形式，即$LU$分解。

2.7 分解步骤

将$\mathbf{L}\mathbf{U}=\mathbf{A}$写成展开形式
$$\left[\begin{array}{cccccc}1& 0& 0& ...& 0& 0\\ l_{21}& 1& 0& ...& 0& 0\\ l_{31}& l_{32}& 1& ...& 0& 0\\ ...& ...& ...& ...& ...& ...\\ l_{N1}& l_{N2}& l_{N3}& ...& l_{NN-1}& 1\end{array}\right]\left[\begin{array}{cccccc}{u}_{11}& u_{12}& u_{13}& ...& u_{1N-1}& u_{1N}\\ 0& u_{22}& u_{23}& ...& u_{2N-1}& u_{2N}\\ 0& 0& u_{33}& ...& u_{3N-1}& u_{3N}\\ ...& ...& ...& ...& ...& ...\\ 0& 0& 0& ...& 0& u_{NN}\end{array}\right]=\left[\begin{array}{cccccc}{a}_{11}& a_{12}& a_{13}& ...& a_{1N-1}& a_{1N}\\ a_{21}& a_{22}& a_{23}& ...& a_{2N-1}& a_{2N}\\ a_{31}& a_{32}& a_{33}& ...& a_{3N-1}& a_{3N}\\ ...& ...& ...& ...& ...& ...\\ a_{N1}& a_{N2}& a_{N3}& ...& a_{NN-1}& a_{NN}\end{array}\right]$$
$\mathbf{L}$的第1行和$\mathbf{U}$的每一列（第1到N列）相乘，得到
$$u_{1j}=a_{1j}j=1,2,3,...,N$$
接下来，让$\mathbf{L}$的第2到第N行与$\mathbf{U}$的第1列相乘，得到
$$l_{i1}{u}_{11}=a_{i1}\Rightarrow l_{i1}=\frac{a_{i1}}{u_{11}}i=2,3,...,N$$
重复该过程，让$\mathbf{L}$的第2行与$\mathbf{U}$的第2到第N列相乘，得到
$$l_{21}{u}_{1j}+u_{2j}=a_{2j}\Rightarrow u_{2j}=a_{2j}-l_{21}{u}_{1j}j=2,3,...,N$$
然后呢，再把$\mathbf{L}$的第3到N行与$\mathbf{U}$的第2列相乘，得到
$$l_{i1}{u}_{12}+l_{i2}{u}_{22}=ai2\Rightarrow l_{i2}=\frac{ai2-l_{i1}{u}_{12}}{u_{22}}i=3,4,...,N$$
如此不断重复，直至解出所有$\mathbf{L}$和$\mathbf{U}$的元素值。

一般来说，$\mathbf{L}$的第i行与$\mathbf{U}$的第i到第N列相乘，可得
$$u_{ij}=a_{ij}-\sum _{k=1}^{i-1}{l}_{ik}{u}_{kj}j=i,i+1,...,N$$
而$\mathbf{L}$的第i+1到第N行与$\mathbf{U}$的第i列相乘，则可得
$$l_{ki}=\frac{a_{ki}-{\sum }_{j=1}^{i-1}{l}_{kj}{u}_{ji}}{u_{ii}}k=i+1,i+2,...,N$$
对于$\mathbf{L}$的第N行，其系数与$\mathbf{U}$的第N列相乘后，得到
$$u_{NN}=a_{NN}-\sum _{k=1}^{N-1}{l}_{Nk}{u}_{kN}$$
LU分解的算法汇总如下

2.8 LU分解算法

for j = 1 to N
	u_1j = a_1j
for i = 2 to N
	l_i1 = a_i1 / u_11
for i = 2 to N-1
{
	for j = i to N
	{
		sum = 0
		for k = 1 to i-1
			sum += l_ik * u_kj
		u_ij = a_ij - sum
	}
	for k = i+1 to N
	{
		sum = 0
		for j = 1 to i-1
			sum += l_kj * u_ji
		l_ki = (a_ki - sum) / u_ii
	}
}
sum = 0
for i = 1 to N-1
	sum += l_Ni * u_iN
u_NN = a_NN - sum

2.9 回代步

把原矩阵$\mathbf{A}$分解为$\mathbf{L}$和$\mathbf{U}$之后，方程系统可分两步来求解（$\mathbf{U}\mathbf{\varphi }=\mathbf{c},\mathbf{L}\mathbf{c}=\mathbf{b}$），两步方法等效于求解两个线性系统方程，然而却由于$\mathbf{L}$和$\mathbf{U}$的上下三角形式而极大简化了求解过程。

第1步是求解$\mathbf{L}\mathbf{c}=\mathbf{b}$，通过前向回代，可以直接得出
$$\begin{array}{rl}& c_1=b_1\\ & c_i=b_i-\sum _{j=1}^{i-1}{l}_{ij}{c}_j,i=2,3,...,N\end{array}$$

第2步是求解$\mathbf{U}\mathbf{\varphi }=\mathbf{c}$，采用反向回代，可以直接得到
$$\begin{array}{rl}& {\varphi }_N=\frac{C_N}{u_{NN}}\\ & {\varphi }_i=\frac{c_i-\sum _{j=i+1}^N{u}_{ij}{\varphi }_j}{u_{ii}},i=N-1,N-2,...,3,2,1\end{array}$$

实际上，$\mathbf{L}$和$\mathbf{U}$的元素可以直接存放在最初的矩阵$\mathbf{A}$中，如果$\mathbf{A}$不再需要的话（当LU分解完成后，$\mathbf{A}$的确是没啥用了，求解的时候只用$\mathbf{L}$和$\mathbf{U}$就行了）。对于$N\times N$维矩阵执行LU分解需要的计算量是$2N^3/3$，是用Gauss消元法求解线性代数方程组的两倍！然而，LU分解的优势就在于，只用做一次分解，就可以把同样的矩阵$\mathbf{A}$应用于不同的右端项$\mathbf{b}$来求解各种问题。然而，这里LU分解引入的另一个目的，是在于用它来引出其它更加高效的迭代解法。

2.10 LU分解和Gauss消元

也许并不明显，但是在Guass消元的过程中，悄然完成了LU分解。前向消元过程生成了上三角矩阵$\mathbf{U}$，在此过程中，下三角矩阵$\mathbf{L}$也悄悄算好了，因为$\mathbf{L}$的元素实际上就是在消元过程中各行乘上的那个因子。下面的算法，将在Gauss消元的基础上，同时完成LU分解。那么用这个算法显然更方便些哈。

2.11 由Gauss消元来做LU分解的算法

如果把$\mathbf{L}$和$\mathbf{U}$放在原来的矩阵$\mathbf{A}$中，那么直接这么做就好了，相当于在Gauss消元的算法中，把系数存放到了下三角位置上，注意，不需要对$\mathbf{b}$做处理，因为LU分解就是为了求解不同的右端项$\mathbf{b}$的

for k = 1 to N-1
{
	for i = k+1 to N
	{
		a_ik = a_ik / a_kk
		for j = k+1 to N
			a_ij = a_ij - a_ik * a_kj
	}
}

如果$\mathbf{L}$和$\mathbf{U}$另外找地方放，还要保留矩阵$\mathbf{A}$的话，那么直接把矩阵$\mathbf{A}$复制一份出来，还用上面的算法就行了。

---

例1 用LU分解方法求解如下线性代数方程组
$$\left[\begin{array}{cccc}3& -1& 0& 0\\ -2& 6& -1& 0\\ 0& -2& 6& -1\\ 0& 0& -2& 7\end{array}\right]\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ {\varphi }_3\\ {\varphi }_4\end{array}\right]=\left[\begin{array}{c}3\\ 4\\ 5\\ -3\end{array}\right]$$

解
求解过程从略，答案直接给出
$$\mathbf{L}=\left[\begin{array}{cccc}1& 0& 0& 0\\ -\frac{2}{3}& 1& 0& 0\\ 0& -\frac{3}{8}& 1& 0\\ 0& 0& -\frac{16}{45}& 1\end{array}\right],\mathbf{U}=\left[\begin{array}{cccc}3& -1& 0& 0\\ 0& \frac{16}{3}& -1& 0\\ 0& 0& \frac{45}{8}& -1\\ 0& 0& 0& \frac{299}{45}\end{array}\right]$$
$$\mathbf{c}=\left[\begin{array}{c}3\\ 6\\ 29/4\\ -19/45\end{array}\right],\mathbf{\varphi }=\left[\begin{array}{c}435/299\\ 408/299\\ 382/299\\ -19/299\end{array}\right]$$

---

2.12 带状稀疏矩阵的直接解法

Gauss消元和LU分解适用于任何方程组系统，当然，它们可以用于求解从结构或非结构网格的守恒方程离散得出的控制方程组。当使用结构网格时，离散过程得到的系统方程的非零元素仅存在于少数几个对角线上，依据不同的离散框架，常出现三对角或五对角矩阵，对于它们的求解可以提出更加高效的方法。实际上相当于把Gauss消元法应用于这类特定的问题中而已。

2.13 三对角矩阵算法（TDMA）

一维问题导出的方程往往具有三对角形式，即
$$a_i{\varphi }_i+b_i{\varphi }_{i+1}+c_i{\varphi }_{i-1}=d_i,i=1,2,3,...,N,c_1=b_N=0$$
对于$i=1$，可直接推得${\varphi }_1$和${\varphi }_2$的关系
$$i=1\Rightarrow a_1{\varphi }_1=-b_1{\varphi }_2+d_1\Rightarrow {\varphi }_1=-\frac{b_1}{a_1}{\varphi }_2+\frac{d_1}{a_2}$$
同样对于$i=2$，可导出${\varphi }_2$和${\varphi }_3$的关系
$$i=2\Rightarrow a_2{\varphi }_2=-b_2{\varphi }_3-c_2{\varphi }_1+d_2\Rightarrow {\varphi }_2=-\frac{a_1{b}_2}{a_1{a}_2-c_2{b}_1}{\varphi }_3+\frac{d_2{a}_1-c_2{d}_1}{a_1{a}_2-c_2{b}_1}$$
同样的操作可用于${\varphi }_3$到${\varphi }_N$，假设${\varphi }_i$可用${\varphi }_{i+1}$表示为
$${\varphi }_i=P_i{\varphi }_{i+1}+Q_i,i=1,2,3,...,N$$
假设对$i-1$上式已知，那么对$i$的推导如下
$$\begin{array}{c}{\varphi }_{i-1}=P_{i-1}{\varphi }_i+Q_{i-1}\\ a_i{\varphi }_i+b_i{\varphi }_{i+1}+c_i{\varphi }_{i-1}=d_i\end{array}\}\Rightarrow {\varphi }_i=-\frac{b_i}{a_i+c_i{P}_{i-1}}{\varphi }_{i+1}+\frac{d_i-c_i{Q}_{i-1}}{a_i+c_i{P}_{i-1}}$$
从而推得$P_i$和$Q_i$为
$$P_i=-\frac{b_i}{a_i+c_i{P}_{i-1}}{Q}_i=\frac{d_i-c_i{Q}_{i-1}}{a_i+c_i{P}_{i-1}}i=2,3,4,...,N$$
对于$i=1$，可直接算得$P_1$和$Q_1$
$$P_1=-\frac{b_1}{a_1}{Q}_1=\frac{d_1}{a_2}$$
而对于$i=N$，由于$b_N=0$，有
$$b_N=0\Rightarrow P_N=0\Rightarrow {\varphi }_N=Q_N$$
TDMA算法总结如下

1. 计算$P_1$和$Q_1$，用
   $$P_1=-\frac{b_1}{a_1}{Q}_1=\frac{d_1}{a_2}$$
2. 计算$P_i$和$Q_i$，$i=2,3,...,N$，用
   $$P_i=-\frac{b_i}{a_i+c_i{P}_{i-1}}{Q}_i=\frac{d_i-c_i{Q}_{i-1}}{a_i+c_i{P}_{i-1}}$$
3. 求出${\varphi }_N$，用
   $${\varphi }_N=Q_N$$
4. 反向依次求出${\varphi }_i$，$i=N-1,N-2,...,3,2,1$，用
   $${\varphi }_i=P_i{\varphi }_{i+1}+Q_i,i=N-1,N-2,...,3,2,1$$

2.14 五对角矩阵算法（PDMA）

还是1维问题，只不过离散格式用的精度更高，让i节点与i+1，i+2，i-1，i-2都有关系，这样子出来的离散矩阵便是5对角形式了，即
$$a_i{\varphi }_i+b_i{\varphi }_{i+2}+c_i{\varphi }_{i+1}+d_i{\varphi }_{i-1}+e_i{\varphi }_{i-2}=f_{i}i=1,2,3,...,N$$
对于边界的几个节点，有
$$\begin{gathered} d_1=e_1=e_2=0 \\ {b}_{N-1}=b_N=c_N=0 \end{gathered}$$
对于$i=1$，有
$${\varphi }_1=-\frac{b_1}{a_1}{\varphi }_3-\frac{c_1}{a_1}{\varphi }_2+\frac{f_1}{a_1}$$
对于$i=2$，有
$${\varphi }_2=-\frac{a_1{b}_2}{a_1{a}_2-d_2{c}_1}{\varphi }_4-\frac{a_1{c}_2-b_1{d}_2}{a_1{a}_2-d_2{c}_1}{\varphi }_3+\frac{a_1{f}_2-d_2{f}_1}{a_1{a}_2-d_2{c}_1}$$
该过程可重复下去，获取第$i$个变量${\varphi }_i$的通有形式
$${\varphi }_i=P_i{\varphi }_{i+2}+Q_i{\varphi }_{i+1}+R_{i}i=1,2,3,...,N$$
在${\varphi }_{i-1}$和${\varphi }_{i-2}$已知的情况下，可求得${\varphi }_i$的形式为
$$\begin{array}{rl}{\varphi }_i=& -\frac{b_i}{a_i+e_i{P}_{i-2}+(d_i+e_i{Q}_{i-2})Q_{i-1}}{\varphi }_{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}}{\varphi }_{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}}\end{array}$$
比较可得
$$\begin{gathered} 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}} \\ {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}} \end{gathered}$$
前两个节点，即$i=1,2$的值为
$$\begin{array}{rl}& P_1=-\frac{b_1}{a_1},Q_1=-\frac{c_1}{a_1},R_1=\frac{f_1}{a_1}\\ & P_2=-\frac{b_2}{a_2+d_2{Q}_1},Q_2=-\frac{c_2+d_2{P}_1}{a_2+d_2{Q}_1},R_2=\frac{f_2-d_2{R}_1}{a_2+d_2{Q}_1}\end{array}$$
由于$b_{N-1}=b_N=c_N=0$，所以$P_{N-1}=P_N=Q_N=0$，那么${\varphi }_N$和${\varphi }_{N-1}$的方程为
$$\begin{array}{rl}& {\varphi }_N=R_N\\ & {\varphi }_{N-1}=Q_{N-1}{\varphi }_N+R_{N-1}\end{array}$$
将PDMA算法的流程汇总如下

1. 计算 $P_1,Q_1,R_1,P_2,Q_2,R_2$，使用
   $$\begin{array}{rl}& P_1=-\frac{b_1}{a_1},Q_1=-\frac{c_1}{a_1},R_1=\frac{f_1}{a_1}\\ & P_2=-\frac{b_2}{a_2+d_2{Q}_1},Q_2=-\frac{c_2+d_2{P}_1}{a_2+d_2{Q}_1},R_2=\frac{f_2-d_2{R}_1}{a_2+d_2{Q}_1}\end{array}$$
2. 对于$i=3,4,5,...,N$计算$P_i,Q_i,R_i$，用
   $$\begin{gathered} 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}} \\ {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}} \end{gathered}$$
3. 计算${\varphi }_N$和${\varphi }_{N-1}$，用
   $$\begin{array}{rl}& {\varphi }_N=R_N\\ & {\varphi }_{N-1}=Q_{N-1}{\varphi }_N+R_{N-1}\end{array}$$
4. 对于$i=N-2,N-3,...,3,2,1$，计算${\varphi }_i$，用
   $${\varphi }_i=P_i{\varphi }_{i+2}+Q_i{\varphi }_{i+1}+R_i$$

3 迭代方法

直接解法通常并不适用于求系数疏矩阵为稀疏情况下大型的方程组，即，大多数矩阵元素是0的情况。这种情况经常会出现，比如线性化的系统方程反映的是非线性的问题（系数是依赖于解的），或是处理与时间相关的系统时。这些也正是求解流动问题时经常碰到的方程类型。

相对而言，迭代方法对于这些问题更具吸引性，因为线性化系统的解变成了迭代求解过程的一部分。而且迭代方法占用内存小、计算消耗低，这都是它比直接解法有优势的地方。

迭代方法有许多系列，对它们的详细回顾可参考相关文献书籍。本章主要介绍下基本的迭代方法，并介绍多重网格算法，其可有效克服迭代方法中的缺陷。Gauss消元和LU分解直接方法的讲解只是为了阐明数值过程中的一些基本概念，以便更加深入地理解迭代方法。

为了统一这些方法的表达形式，系数矩阵将统一写成如下形式
$$\mathbf{A}=\mathbf{D}+\mathbf{L}+\mathbf{U}$$
其中$\mathbf{D}$、$\mathbf{L}$、$\mathbf{U}$分别是对角、严格下三角、严格上三角矩阵。即，相当于把原来系数矩阵$\mathbf{A}$的元素剥离成对角矩阵、下三角矩阵和上三角矩阵（注意，这里的$\mathbf{L}$和$\mathbf{U}$并非是前面的LU分解出来的矩阵，而只是把原始矩阵$\mathbf{A}$的对角线元素、下三角元素、上三角元素拿出来构成的矩阵而已）。

迭代方法在求解线性系统$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$时，计算一系列的解${\mathbf{\varphi }}^{(n)}$，这些解当满足特定的条件时，收敛到精确解$\mathbf{\varphi }$。这样，对于特定解法而言，需要选择一个初值${\mathbf{\varphi }}^{(0)}$（选择为初始条件或初始猜测值），还需要构造一个从${\mathbf{\varphi }}^{(n-1)}$计算到${\mathbf{\varphi }}^{(n)}$的迭代过程。

定点迭代可通过将矩阵$\mathbf{A}$分解为
$$\mathbf{A}=\mathbf{M}-\mathbf{N}$$
这样，原方程$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$转化为
$$(\mathbf{M}-\mathbf{N})\mathbf{\varphi }=\mathbf{b}$$
使用定点迭代过程，上式变为
$$\mathbf{M}{\mathbf{\varphi }}^{(n)}=\mathbf{N}{\mathbf{\varphi }}^{(n-1)}+\mathbf{b}$$
可改写成如下形式
$${\mathbf{\varphi }}^{(n)}=\mathbf{B}{\mathbf{\varphi }}^{(n-1)}+\mathbf{C}\mathbf{b}n=1,2,...$$
其中$\mathbf{B}={\mathbf{M}}^{-1}\mathbf{N}$，而$\mathbf{C}={\mathbf{M}}^{-1}$。这些矩阵的选择不同将产生不同的迭代方法。

在详细描述各种迭代方法之前，先讲下迭代方法需要满足什么样的一些特性，才能让其得到收敛的解。

A. 在收敛的时候，有${\mathbf{\varphi }}^{(n)}={\mathbf{\varphi }}^{(n-1)}=\mathbf{\varphi }$，则迭代方程可以写成
$$\mathbf{\varphi }=\mathbf{B}\mathbf{\varphi }+\mathbf{C}\mathbf{b}$$
即
$${\mathbf{C}}^{-1}(\mathbf{I}-\mathbf{B})\mathbf{\varphi }=\mathbf{b}$$
与原方程$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$相比较，可得
$$\mathbf{A}={\mathbf{C}}^{-1}(\mathbf{I}-\mathbf{B})$$
或
$$\mathbf{B}+\mathbf{C}\mathbf{A}=\mathbf{I}$$
不同矩阵间的关系保证了当达到精确解的时候，再继续迭代不会更改所得值。

B. 从${\mathbf{\varphi }}^0\ne \mathbf{\varphi }$开始，该方法应保证${\mathbf{\varphi }}^n$随着迭代次数$n$的增加最终收敛到$\mathbf{\varphi }$。因为${\mathbf{\varphi }}^n$可以展开成${\mathbf{\varphi }}^0$的如下形式
$${\mathbf{\varphi }}^n={\mathbf{B}}^n{\mathbf{\varphi }}^0+\sum _{i=0}^{n-1}{\mathbf{B}}^i\mathbf{C}\mathbf{b}$$
那么，若要收敛成立，则需要让$\mathbf{B}$满足（即让${\mathbf{\varphi }}^0$的影响完全消失掉）
$$\underset{n\to \infty }{\lim }{\mathbf{B}}^n=\underset{n\to \infty }{\lim }\underset{ntimes}{\underbrace{\mathbf{B}\ast \mathbf{B}\ast \mathbf{B}\ast ...\ast \mathbf{B}}}=0$$
这意味着$\mathbf{B}$的谱半径要小于1，即
$$\rho (\mathbf{B})<1$$
该条件保证了迭代方法是有自我纠正能力的，即，其对于精确解的任何不利的错误扰动是robust（健壮的，鲁棒的）的。

关于该条件的深层理解，可以通过定义精确值与迭代值之间的错误矢量${\mathbf{e}}^{(n)}$，即
$$\begin{gathered} {\mathbf{e}}^{(n)}={\mathbf{\varphi }}^{(n)}-\mathbf{\varphi } \\ {\mathbf{e}}^{(n-1)}={\mathbf{\varphi }}^{(n-1)}-\mathbf{\varphi } \end{gathered}$$
将$\mathbf{\varphi }=\mathbf{B}\mathbf{\varphi }+\mathbf{C}\mathbf{b}$与${\mathbf{\varphi }}^{(n)}=\mathbf{B}{\mathbf{\varphi }}^{(n-1)}+\mathbf{C}\mathbf{b}$两式子互减，可得
$${\mathbf{e}}^{(n)}=\mathbf{B}{\mathbf{e}}^{(n-1)}$$
因此，若方法是收敛的，则应满足如下式子
$$\underset{n\to \infty }{\lim }{\mathbf{e}}^{(n)}=0$$
为使上式更有意义，假设$\mathbf{B}$的特征向量是完全的，且形成一个完备集，即它们是${\mathbf{R}}^N$的基本集。那么$\mathbf{e}$可以用$\mathbf{B}$的$N$个特征向量$\mathbf{v}$的线性组合表示，即
$$\mathbf{e}=\sum _{i=1}^N{\alpha }_i{\mathbf{v}}_i$$
每个特征向量满足
$$\mathbf{B}{\mathbf{v}}_i={\lambda }_i{\mathbf{v}}_i$$
其中${\lambda }_i$是与特征向量${\mathbf{v}}_i$相关的特征值，根据${\mathbf{e}}^{(n)}=\mathbf{B}{\mathbf{e}}^{(n-1)}$，从第1个迭代步开始，有
$${\mathbf{e}}^{(1)}=\mathbf{B}{\mathbf{e}}^{(0)}=\mathbf{B}\sum _{i=1}^N{\alpha }_i{\mathbf{v}}_i=\sum _{i=1}^N{\alpha }_i(\mathbf{B}{\mathbf{v}}_i)=\sum _{i=1}^N{\alpha }_i({\lambda }_i{\mathbf{v}}_i)$$
对于第2次迭代，误差变为
$${\mathbf{e}}^{(2)}=\mathbf{B}{\mathbf{e}}^{(1)}=\mathbf{B}\sum _{i=1}^N{\alpha }_i({\lambda }_i{\mathbf{v}}_i)=\sum _{i=1}^N{\alpha }_i{\lambda }_i(\mathbf{B}{\mathbf{v}}_i)=\sum _{i=1}^N{\alpha }_i{\lambda }_i^2{\mathbf{v}}_i$$
该过程可重复执行，有
$${\mathbf{e}}^{(n)}=\sum _{i=1}^N{\alpha }_i{\lambda }_i^n{\mathbf{v}}_i$$
因此，如果随着$n$趋近无穷大迭代过程要趋于收敛的话，那么所有的特征值必须小于1。如果特征值之中有任何大于1的，那么误差将趋于无穷大。这就很好地解释了为啥谱半径要小于1，因为谱半径实际上是最大特征值，即
$$\rho (\mathbf{B})=\underset{i=1}{\overset{N}{\max }}({\lambda }_i)<1$$
通过减小迭代矩阵的谱半径，也可以提高迭代方法的收敛速度，这也是迭代技术的核心。

C. 迭代方法中需要给出停止判据，应用较多的准则是给予残差的幅值变化，残差如下
$${\mathbf{r}}^{(n)}=\mathbf{A}{\mathbf{\varphi }}^{(n)}-\mathbf{b}$$
一个准则是找寻区域中的最大残差，让其值小于某个阈值$ϵ$定义为收敛条件，即
$$\underset{i=1}{\overset{N}{\max }}\left|b_i-\sum _{j=1}^N{a}_{ij}{\varphi }_j^{(n)}\right|\le ϵ$$
或者root mean square残差（RMS均方根残差）小于$ϵ$，即
$$\frac{\sum _{i=1}^N{\left(b_i-\sum _{j=1}^N{a}_{ij}{\varphi }_j^{(n)}\right)}^2}{N}\le ϵ$$
还有一种准则是让两个连续迭代过程中的最大单位化差值小于$ϵ$，即
$$\underset{i=1}{\overset{N}{\max }}\left|\frac{{\varphi }_i^{(n)}-{\varphi }_i^{(n-1)}}{{\varphi }_i^{(n)}}\right|\times 100\le ϵ$$
顺道说一句，我本人（梅冠华）最喜欢和最常用的判据，是计算两次迭代的差值向量的2范数（幅值），再用本次迭代结果向量的2范数（幅值）去单位化处理，然后让该值小于某个$ϵ\approx 10^{-6}$，个人感觉，这种方法兼顾了向量的各个分量，又抛掉了正负号的影响，还做了单位化处理，是比较合理的，即
$$\frac{\sqrt{{\sum }_{i=1}^N{\left({\varphi }_i^{(n)}-{\varphi }_i^{(n-1)}\right)}^2}}{\sqrt{{\sum }_{i=1}^N{\left({\varphi }_i^{(n)}\right)}^2}+eps}\le ϵ$$
其中的$eps$是一个非常小的数，比如$1e-12$，主要是为了防止分母是0，出现被0除的情况，比如你的初始速度场全是0值的情况。

3.1 Jacobi方法

这恐怕是求解线性方程组系统中最简单的迭代方法了，Jacobi方法的图示如下

对于原方程$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$，其迭代步骤为
$${\varphi }_i^{(n)}=\frac{1}{a_{ii}}\left(b_i-\sum _{j=1,j\ne i}^N{a}_{ij}{\varphi }_j^{(n-1)}\right)i=1,2,3,...,N$$
即，迭代中获得的新值并不用于本次迭代后续值的计算，而是在下次迭代时再去使用。使用矩阵，将上式的展开形式写出来，即
$$\left[\begin{array}{ccccc}{a}_{11}& 0& ...& 0& 0\\ 0& a_{22}& ...& 0& 0\\ ...& ...& ...& ...& ...\\ 0& 0& ...& a_{N-1,N-1}& 0\\ 0& 0& ...& 0& a_{NN}\end{array}\right]\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ ...\\ {\varphi }_{N-1}\\ {\varphi }_N\end{array}\right]+\left[\begin{array}{ccccc}0& a_{12}& ...& a_{1,N-1}& a_{1,N}\\ a_{21}& 0& ...& a_{2,N-1}& a_{2,N}\\ ...& ...& ...& ...& ...\\ a_{N-1,1}& a_{N-1,2}& ...& 0& a_{N-1,N}\\ a_{N,1}& a_{N,2}& ...& a_{N,N-1}& 0\end{array}\right]\left[\begin{array}{c}{\varphi }_1\\ {\varphi }_2\\ ...\\ {\varphi }_{N-1}\\ {\varphi }_N\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_{N-1}\\ b_N\end{array}\right]$$
求解${\mathbf{\varphi }}^{(n)}$，上式转化为
$$\left[\begin{array}{c}{\varphi }_1^{(n)}\\ {\varphi }_2^{(n)}\\ ...\\ {\varphi }_{N-1}^{(n)}\\ {\varphi }_N^{(n)}\end{array}\right]={\left[\begin{array}{ccccc}{a}_{11}& 0& ...& 0& 0\\ 0& a_{22}& ...& 0& 0\\ ...& ...& ...& ...& ...\\ 0& 0& ...& a_{N-1,N-1}& 0\\ 0& 0& ...& 0& a_{NN}\end{array}\right]}^{-1}\left(\left[\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_{N-1}\\ b_N\end{array}\right]-\left[\begin{array}{ccccc}0& a_{12}& ...& a_{1,N-1}& a_{1,N}\\ a_{21}& 0& ...& a_{2,N-1}& a_{2,N}\\ ...& ...& ...& ...& ...\\ a_{N-1,1}& a_{N-1,2}& ...& 0& a_{N-1,N}\\ a_{N,1}& a_{N,2}& ...& a_{N,N-1}& 0\end{array}\right]\left[\begin{array}{c}{\varphi }_1^{(n-1)}\\ {\varphi }_2^{(n-1)}\\ ...\\ {\varphi }_{N-1}^{(n-1)}\\ {\varphi }_N^{(n-1)}\end{array}\right]\right)$$
使用$\mathbf{A}=\mathbf{D}+\mathbf{L}+\mathbf{U}$，可将上式转化为
$${\mathbf{\varphi }}^{(n)}=-{\mathbf{D}}^{-1}(\mathbf{L}+\mathbf{U}){\mathbf{\varphi }}^{(n-1)}+{\mathbf{D}}^{-1}\mathbf{b}$$
那么只要$\rho (-{\mathbf{D}}^{-1}(\mathbf{L}+\mathbf{U}))<1$则Jacobi方法是收敛的。该条件对于很多类型的矩阵都是满足的，其中一种就是对角占优矩阵，即它们的系数满足
$$\sum _{j=1,j\ne i}^N|a_{ij}|\le |a_{ii}|i=1,2,3,...,N$$

3.2 Gauss-Seidel方法

从Jacobi方法改进而来的更好的方法是Gauss-Seidel方法，其收敛特性更好，而且占用内存较少，因为其不需要再单独存储计算出来的新值了，因为它用的都是最新算出来的值。Gauss-Seidel方法的示意图为

Guass-Seidel方法的迭代关系式为
$${\varphi }_i^{(n)}=\frac{1}{a_{ii}}\left(b_i-\sum _{j=1}^{i-1}{a}_{ij}{\varphi }_j^{(n)}-\sum _{j=i+1}^N{a}_{ij}{\varphi }_j^{(n-1)}\right)i=1,2,3,...,N$$
如果要写成矩阵形式，即
$${\mathbf{\varphi }}^{(n)}=-(\mathbf{D}+\mathbf{L}{)}^{-1}\mathbf{U}{\mathbf{\varphi }}^{(n-1)}+(\mathbf{D}+\mathbf{L}{)}^{-1}\mathbf{b}$$
Gauss-Seidel方法更为高效，因为其计算${\varphi }_j^{(n)}$的时候用的是最新算出的${\varphi }_i^{(n)},i=1,2,3,...,j-1$。该方法很节省内存，因为直接用新值覆盖旧值即可。Gauss-Seidel迭代收敛的要求是
$$\rho (-(\mathbf{D}+\mathbf{L}{)}^{-1}\mathbf{U})<1$$
尽管在某些情况下Jacobi方法收敛更快，然而Gauss-Seidel仍然是首选方法。

---

例2

对于例1给出的系统方程，采用Jacobi和Gauss-Seidel方法，迭代5次，并计算每次迭代解与真实解的误差，迭代初值选为
$${\mathbf{\varphi }}^{\ast }=[0000{]}^T$$
解
较为简单，略去不提，感兴趣的可自行计算。

---

3.3 预处理和迭代方法

迭代方法$${\mathbf{\varphi }}^{(n)}=\mathbf{B}{\mathbf{\varphi }}^{(n-1)}+\mathbf{C}\mathbf{b}n=1,2,...$$的收敛速度依赖于迭代矩阵$\mathbf{B}$的谱半径，这与矩阵的系数相关。尽管迭代方法是把原系统方程转化为有相同解的等效系统，但是拥有更好的谱半径特性。在这些条件下，等效系统的特征值更加聚集，使得其比原系统的迭代收敛速度更快。预处理器就是这样定义的一种矩阵，即
$${\mathbf{P}}^{-1}\mathbf{A}\mathbf{\varphi }={\mathbf{P}}^{-1}\mathbf{b}$$
其有着和原方程$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$同样的解，但是${\mathbf{P}}^{-1}\mathbf{A}$的谱半径特性更好。定义预处理器$\mathbf{P}$的时候，其难点在于找寻一个与${\mathbf{A}}^{-1}$近似的矩阵，而且还要很容易做求逆运算（为了计算${\mathbf{P}}^{-1}$）。

将之前的式子$\mathbf{A}=\mathbf{M}-\mathbf{N}$中的$\mathbf{M}$替换成$\mathbf{P}$，则$\mathbf{M}=\mathbf{P}$且$\mathbf{A}=\mathbf{P}-\mathbf{N}$，那么与之相关的迭代系统变为
$$\begin{array}{rl}{\mathbf{\varphi }}^{(n)}& =\mathbf{B}{\mathbf{\varphi }}^{(n-1)}+\mathbf{C}\mathbf{b}\\ & ={\mathbf{P}}^{-1}\mathbf{N}{\varphi }^{(n-1)}+{\mathbf{P}}^{-1}\mathbf{b}\\ & ={\mathbf{P}}^{-1}(\mathbf{P}-\mathbf{A}){\varphi }^{(n-1)}+{\mathbf{P}}^{-1}\mathbf{b}\\ & =(\mathbf{I}-{\mathbf{P}}^{-1}\mathbf{A}){\varphi }^{(n-1)}+{\mathbf{P}}^{-1}\mathbf{b}\end{array}$$
写成残差形式，为
$$\begin{array}{rl}{\mathbf{\varphi }}^{(n)}& =(\mathbf{I}-{\mathbf{P}}^{-1}\mathbf{A}){\varphi }^{(n-1)}+{\mathbf{P}}^{-1}\mathbf{b}\\ & ={\varphi }^{(n-1)}+{\mathbf{P}}^{-1}(\mathbf{b}-\mathbf{A}){\varphi }^{(n-1)}\\ & ={\varphi }^{(n-1)}+{\mathbf{P}}^{-1}{\mathbf{r}}^{(n-1)}\end{array}$$
从这两个方程中我们可以看出，迭代方法只不过是做了分解$\mathbf{A}=\mathbf{P}-\mathbf{N}$的预处理系统的定点迭代方法，其中的谱半径特性为
$$\rho (\mathbf{I}-{\mathbf{P}}^{-1}\mathbf{A})<1$$
作为比较，对于Jacobi和Gauss-Seidel的预处理矩阵J和JS分别为
$$\begin{gathered} {\mathbf{P}}_J=\mathbf{D} \\ {\mathbf{P}}_{GS}=\mathbf{D}+\mathbf{L} \end{gathered}$$
其中$\mathbf{D}$和$\mathbf{L}$分别是矩阵$\mathbf{A}$的对角和下三角部分。

这样，预处理是对原系统的一种提高谱半径特性的处理方式，其运用预处理矩阵$\mathbf{P}$来提高迭代效率。下面我们将会看到，可以提出更加优越的预处理矩阵，其系数的定义更加复杂。

3.4 矩阵分解技术

Gauss-Seidel和Jacobi方法的低收敛速度，激发人们寻求更加快速的迭代方法。一种加速迭代收敛速度的方法是通过使用更好的预处理器来开发迭代方法，一个简单和行之有效的方法是对原矩阵$\mathbf{A}$执行不完全分解，之所以是不完全分解，因为完全分解的话，把矩阵$\mathbf{A}$分解成下三角矩阵$\mathbf{L}$和上三角矩阵$\mathbf{U}$的工作量跟直接求解方程是相当的，而且计算耗费很大，内存消耗很大（丧失了稀疏特性）。

3.5 不完全LU分解（ILU）

如同在例1中看到的那样，$\mathbf{L}$和$\mathbf{U}$矩阵会出现原矩阵$\mathbf{A}$中$0$元素地方的非$0$元素（这是所谓的fill-in填补）。所以呢，如果是要做不完全LU（ILU）分解的话，那么得到的下三角$\mathbf{L}$和上三角$\mathbf{U}$矩阵应该和原矩阵$\mathbf{A}$中下三角和上三角的非零元素排布完全一致才好，即
$$\mathbf{A}=\mathbf{L}\mathbf{U}+\mathbf{R}$$
其中$\mathbf{R}$是分解过程产生的残差，矩阵$\mathbf{L}$和$\mathbf{U}$是稀疏矩阵（与矩阵$\mathbf{A}$的排布一样），这样处理起来非常方便，相比而言，如果是做完全分解的话，处理起来就麻烦多了。然而，$\mathbf{L}\mathbf{U}$只是矩阵$\mathbf{A}$的近似，就需要采用迭代求解过程来求解方程组系统了。求解过程得第1步是把方程组$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$重写成
$$\mathbf{A}\mathbf{\varphi }=\mathbf{b}\Rightarrow 0=\mathbf{b}-\mathbf{A}\mathbf{\varphi }\Rightarrow (\mathbf{A}-\mathbf{R})\mathbf{\varphi }=(\mathbf{A}-\mathbf{R})\mathbf{\varphi }+(\mathbf{b}-\mathbf{A}\mathbf{\varphi })$$
迭代形式可以写成
$$(\mathbf{A}-\mathbf{R}){\mathbf{\varphi }}^{(n)}=(\mathbf{A}-\mathbf{R}){\mathbf{\varphi }}^{(n-1)}+(\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(n-1)})$$
上式常写成修正形式，即定义
$${\mathbf{\varphi }}^{(n)}={\mathbf{\varphi }}^{(n-1)}+{\mathbf{\varphi }}^{\prime (n)}$$
这样迭代过程变为
$$(\mathbf{A}-\mathbf{R}){\mathbf{\varphi }}^{\prime (n)}=(\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(n-1)})$$
即每步迭代中，用上式算得${\mathbf{\varphi }}^{\prime (n)}$，然后用${\mathbf{\varphi }}^{(n)}={\mathbf{\varphi }}^{(n-1)}+{\mathbf{\varphi }}^{\prime (n)}$算得${\mathbf{\varphi }}^{(n)}$。

ILU分解可以通过Gauss消元法，并去掉一些在预设位置的非对角元素来实现，那么舍弃元素所处的位置不同便催生了不同的ILU近似。

3.6 不完全LU分解不含填补（ILU(0)）

有许多ILU分解技术的变种，然而最简单的莫过于ILU(0)了，在ILU(0)中，0元素在矩阵$\mathbf{L}$和$\mathbf{U}$中的位置是和原矩阵$\mathbf{A}$中的0元素位置严格对应的。使用Gauss消元，计算过程和完全LU分解是一样的，只是如果在原矩阵$\mathbf{A}$的0元素位置处出现了新的非0元素（$l_{ij}$和$u_{ij}$），那么直接把它们刨掉就是了。因此，矩阵$\mathbf{L}$和$\mathbf{U}$联合起来，所反映的就是原始矩阵$\mathbf{A}$中非0元素的位置，如此一来，在分解稀疏矩阵时经常出现的填补现象就被消除了。然而，精度却有所丧失导致达到收敛的迭代次数有所增加。为克服该缺陷，人们发展了更加精确的ILU分解方法，它们更加高效更加可靠。这些方法依据允许填补层次的不同，用$ILU(p)$来定义填补的阶数。填补的阶数越高，ILU分解步骤代价越大。此外，当使用多重网格方法的时候（后面会讲到），ILU(0)方法作为光滑器（smoother）是绰绰有余的。正因为如此，这里就不展示高阶方法了。

ILU(0)分解算法假设$\mathbf{L}$是个单位下三角矩阵，因此，可以把矩阵$\mathbf{L}$和$\mathbf{U}$还是存在矩阵$\mathbf{A}$原本的位置上，该算法如下

3.7 ILU(0)分解算法

for k = 1 to N-1
{
	for i = k+1 to N and a_ik != 0
	{
		a_ik = a_ik / a_kk	// 这个是l矩阵的元素
		for j = k+1 to N and a_ij != 0
			a_ij = a_ij - a_ik * a_kj	// 这个是u矩阵的元素
	}
}

其实跟Gauss消元推出的LU分解一样，不过是判断了原矩阵的原始元素是否为0，如果是0就不处理了，仅此而已。

值的注意的是，对于对称正定矩阵的ILU分解被定义成了不完全Cholesky分解，此时，分解只需做下三角（或上三角）部分，对原矩阵的近似可以写为
$$\overline{\mathbf{L}}{\overline{\mathbf{L}}}^{\text{T}}\approx \mathbf{A}$$
其中$\overline{\mathbf{L}}$是分解出来的稀疏下三角矩阵（矩阵$\mathbf{L}$的近似），预处理矩阵$\mathbf{P}$为
$$\mathbf{P}=\overline{\mathbf{L}}{\overline{\mathbf{L}}}^{\text{T}}\approx \mathbf{A}$$

3.8 ILU分解预处理器

基于ILU分解有很多好用的预处理器，前面提到直接解法中把稀疏矩阵$\mathbf{A}$分解成下三角和上三角矩阵乘积的时候会出现填补现象。由于预处理器仅需要是矩阵${\mathbf{A}}^{-1}$的近似，因此，寻找矩阵$\mathbf{A}$的近似分解就足够了，比如$\mathbf{A}=\overline{\mathbf{L}}\overline{\mathbf{U}}$。选择$\mathbf{P}=\overline{\mathbf{L}}\overline{\mathbf{U}}$将会使得${\mathbf{P}}^{-1}$的计算非常高效，因为该求逆可以用前向和后向回代来求出，如前所述，而此时精确的$\mathbf{L}$和$\mathbf{U}$已经分别被近似的$\overline{\mathbf{L}}$和$\overline{\mathbf{U}}$替代了。

对于ILU(0)方法，不完全分解减少了非零元素的稀疏度，预处理器与原矩阵的大小完全一样。为了减少存储量，Pommerell引入了ILU的简化版叫做对角ILU(DILU)，在DILU中非对角元的填补都被消去了（矩阵中的上下三角部分的元素保持不变），只有对角元素才被修改。

那么，预处理器可写成如下形式
$$\mathbf{P}=({\mathbf{D}}^{\ast }+\mathbf{L}){\mathbf{D}}^{\ast -1}({\mathbf{D}}^{\ast }+\mathbf{U})$$
其中$\mathbf{L}$和$\mathbf{U}$为矩阵$\mathbf{A}$的上下三角分量（注意，不是LU分解，这里的L和U就是A中原本的上下三角部分的元素值），而${\mathbf{D}}^{\ast }$为规则的对角矩阵，但是与$\mathbf{A}$的对角元素不同，矩阵${\mathbf{D}}^{\ast }$的定义是，让上式$({\mathbf{D}}^{\ast }+\mathbf{L}){\mathbf{D}}^{\ast -1}({\mathbf{D}}^{\ast }+\mathbf{U})$的对角元等于$\mathbf{A}$的对角元素。

3.9 计算DILU方法中D*的算法

for i = 1 to N
	d_ii = a_ii
for i = 1 to N
	for j = i+1 to N 且 a_ij != 0 且 a_ji != 0
		d_jj = d_jj - a_ji * a_ij / d_ii

为了求预处理器的逆，定义
$$\mathbf{P}=({\mathbf{D}}^{\ast }+\mathbf{L}){\mathbf{D}}^{\ast -1}({\mathbf{D}}^{\ast }+\mathbf{U})=\overline{\mathbf{L}}\overline{\mathbf{U}},\overline{\mathbf{L}}=({\mathbf{D}}^{\ast }+\mathbf{L}){\mathbf{D}}^{\ast -1},\overline{\mathbf{U}}=({\mathbf{D}}^{\ast }+\mathbf{U})$$
或
$$\mathbf{P}=({\mathbf{D}}^{\ast }+\mathbf{L})(\mathbf{I}+{\mathbf{D}}^{\ast -1}\mathbf{U})=\overline{\mathbf{L}}\overline{\mathbf{U}},\overline{\mathbf{L}}=({\mathbf{D}}^{\ast }+\mathbf{L}),\overline{\mathbf{U}}=(\mathbf{I}+{\mathbf{D}}^{\ast -1}\mathbf{U})$$
这样在求解$\mathbf{P}{\mathbf{\varphi }}^{\prime (n+1)}={\mathbf{r}}^{(n)}$时，需要用到修正场${\mathbf{\varphi }}^{\prime (n+1)}={\mathbf{P}}^{-1}{\mathbf{r}}^{(n)}$，就可以很方便地用前向和后向回代算法求出了。

3.10 DILU方法中的前向和后向回代算法

for i = 1 to N
	for j = 1 to i-1
		t_i = d_ii^-1 (r_i - l_ij * t_j)
for i = N to 1
	for j = i+1 to N
		phi'_i = t_i - d_ii^-1 (u_ij * t_j)

DILU的显著优势是，其仅需要存储一个额外的对角矩阵就行了。

3.11 求解代数系统的梯度方法

另一组求解线性代数方程组系统的迭代方法是梯度方法（Gradient Method），包含最速下降法（Steepest Descent）和共轭梯度方法（Conjugate Gradient）。它们最初是用于求解对称正定（SPD）矩阵$\mathbf{A}$为系数的二次矢量函数的最小值问题，即
$$\mathbf{Q}(\mathbf{\varphi })=\frac{1}{2}{\mathbf{\varphi }}^T\mathbf{A}\mathbf{\varphi }-{\mathbf{b}}^T\mathbf{\varphi }+\mathbf{c}$$
其中$\mathbf{c}$为标量矢量，$\mathbf{Q}(\mathbf{\varphi })$的最小值对应的是其对$\mathbf{\varphi }$的梯度值是零。而矢量场$\mathbf{Q}(\mathbf{\varphi })$在$\mathbf{\varphi }$处的梯度${\mathbf{Q}}^{\prime }(\mathbf{\varphi })$，指向的正是$\mathbf{Q}(\mathbf{\varphi })$的最大增长方向。梯度可算得
$${\mathbf{Q}}^{\prime }(\mathbf{\varphi })=\frac{1}{2}{\mathbf{A}}^T\mathbf{\varphi }+\frac{1}{2}\mathbf{A}\mathbf{\varphi }-\mathbf{b}$$
如果$\mathbf{A}$是对称的（$\mathbf{A}={\mathbf{A}}^T$），则
$${\mathbf{Q}}^{\prime }(\mathbf{\varphi })=\mathbf{A}\mathbf{\varphi }-\mathbf{b}$$
最小值是在${\mathbf{Q}}^{\prime }(\mathbf{\varphi })=0$处获得的，有
$${\mathbf{Q}}^{\prime }(\mathbf{\varphi })=0\Rightarrow \mathbf{A}\mathbf{\varphi }=\mathbf{b}$$
因此，$\mathbf{Q}(\mathbf{\varphi })$的最小值问题对应的就是求解$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$的问题，即该最小值问题的解和线性代数方程组的解是一样的。

还有个问题，要保证函数$\mathbf{Q}(\mathbf{\varphi })$的最小值是全局最小值，需要系数矩阵$\mathbf{A}$是正定的（positive definite），即，其应该对所有的$\mathbf{\varphi }\ne 0$满足${\mathbf{\varphi }}^T\mathbf{A}\mathbf{\varphi }>0$。这一要求可通过考察精确解$\mathbf{\varphi }$与其当前预估值${\mathbf{\varphi }}^{(n)}$的关系来阐明，令$\mathbf{e}={\mathbf{\varphi }}^{(n)}-\mathbf{\varphi }$代表精确值和当前值的差别，则
$$\begin{array}{rl}\mathbf{Q}(\mathbf{\varphi }+\mathbf{e})& =\frac{1}{2}(\mathbf{\varphi }+\mathbf{e}{)}^T\mathbf{A}(\mathbf{\varphi }+\mathbf{e})-{\mathbf{b}}^T(\mathbf{\varphi }+\mathbf{e})+\mathbf{c}\\ & =\frac{1}{2}{\mathbf{\varphi }}^T\mathbf{A}\mathbf{\varphi }+\frac{1}{2}{\mathbf{e}}^T\mathbf{A}\mathbf{\varphi }+\frac{1}{2}{\mathbf{\varphi }}^T\mathbf{A}\mathbf{e}+\frac{1}{2}{\mathbf{e}}^T\mathbf{A}\mathbf{e}-{\mathbf{b}}^T\mathbf{\varphi }-{\mathbf{b}}^T\mathbf{e}+\mathbf{c}\\ & =\underset{\mathbf{Q}(\mathbf{\varphi })}{\underbrace{\frac{1}{2}{\mathbf{\varphi }}^T\mathbf{A}\mathbf{\varphi }-{\mathbf{b}}^T\mathbf{\varphi }+\mathbf{c}}}+\frac{1}{2}({\mathbf{e}}^T\mathbf{A}\mathbf{\varphi }+{\mathbf{\varphi }}^T\mathbf{A}\mathbf{e})-{\mathbf{b}}^T\mathbf{e}+\frac{1}{2}{\mathbf{e}}^T\mathbf{A}\mathbf{e}\\ & =\mathbf{Q}(\mathbf{\varphi })+\frac{1}{2}{\mathbf{e}}^T\mathbf{A}\mathbf{e}\end{array}$$
之所以消去了一些项，是因为
$$\begin{array}{rl}\frac{1}{2}({\mathbf{e}}^T\mathbf{A}\mathbf{\varphi }+{\mathbf{\varphi }}^T\mathbf{A}\mathbf{e})-{\mathbf{b}}^T\mathbf{e}& =\frac{1}{2}((\mathbf{A}\mathbf{\varphi }{)}^T\mathbf{e}+{\mathbf{\varphi }}^T\mathbf{A}\mathbf{e})-{\mathbf{b}}^T\mathbf{e}\\ & =\frac{1}{2}({\mathbf{\varphi }}^T{\mathbf{A}}^T\mathbf{e}+{\mathbf{\varphi }}^T\mathbf{A}\mathbf{e})-{\mathbf{b}}^T\mathbf{e}\\ & =\frac{1}{2}({\mathbf{\varphi }}^T\mathbf{A}\mathbf{e}+{\mathbf{\varphi }}^T\mathbf{A}\mathbf{e})-{\mathbf{b}}^T\mathbf{e}\\ & ={\mathbf{\varphi }}^T\mathbf{A}\mathbf{e}-{\mathbf{b}}^T\mathbf{e}\\ & ={\mathbf{\varphi }}^T{\mathbf{A}}^T\mathbf{e}-{\mathbf{b}}^T\mathbf{e}\\ & =({\mathbf{\varphi }}^T{\mathbf{A}}^T-{\mathbf{b}}^T)\mathbf{e}\\ & =(\mathbf{A}\mathbf{\varphi }-\mathbf{b}{)}^T\mathbf{e}\\ & =0\end{array}$$
由$\mathbf{Q}(\mathbf{\varphi }+\mathbf{e})=\mathbf{Q}(\mathbf{\varphi })+\frac{1}{2}{\mathbf{e}}^T\mathbf{A}\mathbf{e}$可见，如果$\mathbf{A}$是正定的，那么第2项将永远是正值，除非$\mathbf{e}=0$，即获得了所要的全局最小值。此外，如果$\mathbf{A}$是正定的，它的所有特征值都是正的，且函数$\mathbf{Q}(\mathbf{\varphi })$有唯一的最小值。

这样对于一个对称正定矩阵，可获得其收敛序列${\mathbf{\varphi }}^{(n)}$如下
$${\mathbf{\varphi }}^{(n+1)}={\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}(\delta {\mathbf{\varphi }}^{(n)})$$
其中${\alpha }^{(n)}$是松弛因子，而$\delta {\mathbf{\varphi }}^{(n)}$是与每步迭代过程中最小化函数$\mathbf{Q}(\mathbf{\varphi })$的相关修正，当然这个式子的实现可以有很多方式，从而产生了不同的方法。

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

求解线性方程组系统的最速下降方法，是基于最小化式子$\mathbf{Q}(\mathbf{\varphi })=\frac{1}{2}{\mathbf{\varphi }}^T\mathbf{A}\mathbf{\varphi }-{\mathbf{b}}^T\mathbf{\varphi }+\mathbf{c}$中平方项的考虑。如果$\mathbf{\varphi }$是一维矢量，其分量为标量$\varphi$，那么函数$\mathbf{Q}(\varphi )$代表了一个抛物线，那么从某个初始点${\varphi }_0$开始，找到该抛物线函数的最小值的方法，可以沿着该抛物线向下运动，直到碰到其底部最小值。

同样的思想可以应用于$N$维系统，此时函数$\mathbf{Q}({\mathbf{\varphi }}^{(n)})$描述的可能是一个抛物面，而从一个初始位置${\mathbf{\varphi }}^{(0)}$开始的迭代求解过程，就可以是沿着该抛物面不断向下运动直到达到底部的最小值。为了获取更快的收敛序列，${\mathbf{\varphi }}^{(0)},{\mathbf{\varphi }}^{(1)},{\mathbf{\varphi }}^{(2)},...$应该选择最快速的下降方向，即，-${\mathbf{Q}}^{\prime }(\mathbf{\varphi })$。这个方向就是
$$-{\mathbf{Q}}^{\prime }(\mathbf{\varphi })=\mathbf{b}-\mathbf{A}\mathbf{\varphi }$$
精确解是$\mathbf{\varphi }$，则在第$n$步的误差${\mathbf{e}}^{(n)}$和残差${\mathbf{r}}^{(n)}$，为
$$\begin{array}{c}{\mathbf{e}}^{(n)}={\mathbf{\varphi }}^{(n)}-\mathbf{\varphi }\\ {\mathbf{r}}^{(n)}=\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(n)}=-{\mathbf{Q}}^{\prime }({\mathbf{\varphi }}^{(n)})\end{array}\}\Rightarrow {\mathbf{r}}^{(n)}=\mathbf{A}\mathbf{\varphi }-\mathbf{A}{\mathbf{\varphi }}^{(n)}=-\mathbf{A}{\mathbf{e}}^{(n)}$$
那么，沿着最速下降方向，${\mathbf{\varphi }}^{(n+1)}$的值可由${\mathbf{\varphi }}^{(n)}$的值写作
$${\mathbf{\varphi }}^{(n+1)}={\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{r}}^{(n)}$$
${\alpha }^{(n)}$的值是要让$\mathbf{Q}(\mathbf{\varphi })$最小，所以
$$\frac{d\left(\mathbf{Q}({\mathbf{\varphi }}^{(n+1)})\right)}{d{\alpha }^{(n)}}=0$$
即
$$\frac{d\left(\mathbf{Q}({\mathbf{\varphi }}^{(n+1)})\right)}{d{\alpha }^{(n)}}=0\Rightarrow {\left[\frac{d\left(\mathbf{Q}({\mathbf{\varphi }}^{(n+1)})\right)}{d{\mathbf{\varphi }}^{(n+1)}}\right]}^T\frac{d{\mathbf{\varphi }}^{(n+1)}}{d{\alpha }^{(n)}}\Rightarrow ({\mathbf{r}}^{(n+1)}{)}^T{\mathbf{r}}^{(n)}=0$$
表明新迭代步的方向应该与老迭代步的方向相互垂直，${\alpha }^{(n)}$的值可由上式计算
$$\begin{array}{rl}({\mathbf{r}}^{(n+1)}{)}^T{\mathbf{r}}^{(n)}=0& \Rightarrow (\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(n+1)}{)}^T{\mathbf{r}}^{(n)}=0\\ & \Rightarrow [\mathbf{b}-\mathbf{A}({\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{r}}^{(n)}){]}^T{\mathbf{r}}^{(n)}=0\\ & \Rightarrow (\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(n)}{)}^T{\mathbf{r}}^{(n)}={\alpha }^{(n)}(\mathbf{A}{\mathbf{r}}^{(n)}{)}^T{\mathbf{r}}^{(n)}\\ & \Rightarrow ({\mathbf{r}}^{(n)}{)}^T{\mathbf{r}}^{(n)}={\alpha }^{(n)}({\mathbf{r}}^{(n)}{)}^T\mathbf{A}{\mathbf{r}}^{(n)}\\ & \Rightarrow {\alpha }^{(n)}=\frac{({\mathbf{r}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}{({\mathbf{r}}^{(n)}{)}^T\mathbf{A}{\mathbf{r}}^{(n)}}\end{array}$$
将最速下降算法的流程总结如下

1. 选择残差作为起始方向 ${\mathbf{r}}^{(0)}=\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(0)}$
2. 第n步的迭代过程如下
   2.1. 计算残差矢量${\mathbf{r}}^{(n)}=\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(n)}$
   2.2. 计算正交方向的比率${\alpha }^{(n)}=\frac{({\mathbf{r}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}{({\mathbf{r}}^{(n)}{)}^T\mathbf{A}{\mathbf{r}}^{(n)}}$
   2.3. 获取新值${\mathbf{\varphi }}^{(n+1)}={\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{r}}^{(n)}$

算法中的每个迭代步需要执行两次矩阵-矢量相乘，其中一个可以消去，即
$$\begin{array}{rl}{\mathbf{\varphi }}^{(n+1)}={\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{r}}^{(n)}& \Rightarrow \mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(n+1)}=\mathbf{b}-\mathbf{A}({\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{r}}^{(n)})\\ & \Rightarrow {\mathbf{r}}^{(n+1)}=\mathbf{r}(n)-{\alpha }^{(n)}\mathbf{A}{\mathbf{r}}^{(n)}\end{array}$$
也就是说，${\mathbf{r}}^{(n)}$只需要计算第1步的${\mathbf{r}}^{(0)}$，用${\mathbf{r}}^{(n+1)}=\mathbf{r}(n)-{\alpha }^{(n)}\mathbf{A}{\mathbf{r}}^{(n)}$就可以把后面的算出来了，这样子就不需要再计算$\mathbf{A}{\mathbf{\varphi }}^{(n)}$，因为算$\mathbf{A}{\mathbf{r}}^{(n)}$就好了（残差矢量和正交方向比率都用这个来算，只用算一次就行了）。然而这样做的缺点在于缺乏从${\mathbf{\varphi }}^{(n)}$向残差的反馈，有可能会因为圆整误差导致解收敛到一个不同于精确解的值，该缺陷可通过使用最初的方程定期计算残差来修复。

3.13 共轭梯度方法（The Conjugate Gradient Method）

虽然最速下降方法保证了收敛性，然而其收敛速度非常低。之所以会这样，是由于局部最小值附近的震荡导致该方法不断地在同一个方向搜索引起的，为避免这类问题，需要让每一次搜索的方向都跟前一次不一样才好。这可以通过选择一系列搜索方向${\mathbf{d}}^{(0)},{\mathbf{d}}^{(1)},{\mathbf{d}}^{(2)},...,{\mathbf{d}}^{(N-1)}$，让其是$\mathbf{A}$正交的，两个矢量${\mathbf{d}}^{(m)}$和${\mathbf{d}}^{(n)}$如果满足如下关系则说它们是$\mathbf{A}$正交的
$$({\mathbf{d}}^{(n)}{)}^T\mathbf{A}{\mathbf{d}}^{(m)}=0$$
如果在每次搜索方向都采用了正确的步长，那么解将会在第$N$步收敛。

第n+1步的值为
$${\mathbf{\varphi }}^{(n+1)}={\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{d}}^{(n)}$$
将$\mathbf{\varphi }$从上式两端减去，可得误差方程
$${\mathbf{e}}^{(n+1)}={\mathbf{e}}^{(n)}+{\alpha }^{(n)}{\mathbf{d}}^{(n)}$$
由于在前面最速下降算法中提到${\mathbf{r}}^{(n)}=-\mathbf{A}{\mathbf{e}}^{(n)}$，故
$$\begin{array}{rl}{\mathbf{r}}^{(n+1)}& =-\mathbf{A}{\mathbf{e}}^{(n+1)}\\ & =-\mathbf{A}({\mathbf{e}}^{(n)}+{\alpha }^{(n)}{\mathbf{d}}^{(n)})\\ & ={\mathbf{r}}^{(n)}-{\alpha }^{(n)}\mathbf{A}{\mathbf{d}}^{(n)}\end{array}$$
上式表明，新残差${\mathbf{r}}^{(n+1)}$只是之前残差${\mathbf{r}}^{(n)}$和$\mathbf{A}{\mathbf{d}}^{(n)}$的线性组合。

接下来需要让${\mathbf{e}}^{(n+1)}$和${\mathbf{d}}^{(n)}$是$\mathbf{A}$正交的，这等效于沿着${\mathbf{d}}^{(n)}$方向搜索找到最小点，可推出${\alpha }^{(n)}$为
$$\begin{array}{rl}({\mathbf{d}}^{(n)}{)}^T\mathbf{A}{\mathbf{e}}^{(n+1)}=0& \Rightarrow ({\mathbf{d}}^{(n)}{)}^T\mathbf{A}({\mathbf{e}}^{(n)}+{\alpha }^{(n)}{\mathbf{d}}^{(n)})=0\\ & \Rightarrow {\alpha }^{(n)}=\frac{({\mathbf{d}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}{({\mathbf{d}}^{(n)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}}\end{array}$$
上式也意味着
$$({\mathbf{d}}^{(n)}{)}^T\mathbf{A}{\mathbf{e}}^{(n+1)}=0\Rightarrow ({\mathbf{d}}^{(n)}{)}^T{\mathbf{r}}^{(n+1)}=0$$
即，如果知道了搜索方向，那么可以算出${\alpha }^{(n)}$。

为了推导出搜索方向，假设其控制方程的形式为
$${\mathbf{d}}^{(n+1)}={\mathbf{r}}^{(n+1)}+{\beta }^{(n)}{\mathbf{d}}^{(n)}$$
$\mathbf{d}$矢量的$\mathbf{A}$正交需要其满足
$$({\mathbf{d}}^{(n+1)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}=0$$
结合上上式子，可得
$${\beta }^{(n)}=-\frac{({\mathbf{r}}^{(n+1)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}}{({\mathbf{d}}^{(n)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}}$$
由于${\mathbf{r}}^{(n+1)}={\mathbf{r}}^{(n)}-{\alpha }^{(n)}\mathbf{A}{\mathbf{d}}^{(n)}$，得到
$$\mathbf{A}{\mathbf{d}}^{(n)}=-\frac{1}{{\alpha }^{(n)}}({\mathbf{r}}^{(n+1)}-{\mathbf{r}}^{(n)})$$
结合上式、上上式，以及${\alpha }^{(n)}=\frac{({\mathbf{d}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}{({\mathbf{d}}^{(n)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}}$，可推得
$$\begin{array}{rl}{\beta }^{(n)}& =\frac{({\mathbf{r}}^{(n+1)}{)}^T({\mathbf{r}}^{(n+1)}-{\mathbf{r}}^{(n)})}{({\mathbf{d}}^{(n)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}{\alpha }^{(n)}}\\ & =\frac{({\mathbf{r}}^{(n+1)}{)}^T({\mathbf{r}}^{(n+1)}-{\mathbf{r}}^{(n)})}{({\mathbf{d}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}\\ & =\frac{({\mathbf{r}}^{(n+1)}{)}^T{\mathbf{r}}^{(n+1)}-\underset{=0}{\underbrace{{\mathbf{r}}^{(n+1){)}^T}{\mathbf{r}}^{(n)}}}}{({\mathbf{d}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}\\ & =\frac{({\mathbf{r}}^{(n+1)}{)}^T{\mathbf{r}}^{(n+1)}}{（{\mathbf{d}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}\end{array}$$
上式中的分母可进一步展开为
$$\begin{array}{rl}({\mathbf{d}}^{(n)}{)}^T{\mathbf{r}}^{(n)}& =({\mathbf{r}}^{(n)}+{\beta }^{(n-1)}{\mathbf{d}}^{(n-1)}{)}^T{\mathbf{r}}^{(n)}\\ & =({\mathbf{r}}^{(n)}{)}^T{\mathbf{r}}^{(n)}+{\beta }^{(n-1)}\underset{=0}{\underbrace{({\mathbf{d}}^{(n-1)}{)}^T{\mathbf{r}}^{(n)}}}\\ & =({\mathbf{r}}^{(n)}{)}^T{\mathbf{r}}^{(n)}\end{array}$$
因此，有
$${\beta }^{(n)}=\frac{({\mathbf{r}}^{(n+1)}{)}^T{\mathbf{r}}^{(n+1)}}{({\mathbf{r}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}$$
共轭梯度（CG）算法流程为

1. 选择残差作为初始方向 ${\mathbf{d}}^{(0)}={\mathbf{r}}^{(0)}=\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(0)}$
2. 第n步的迭代过程如下
   2.1. 选择$\mathbf{d}$方向的比率因子${\alpha }^{(n)}=\frac{({\mathbf{d}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}{({\mathbf{d}}^{(n)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}}$
   2.2. 获得新的$\mathbf{\varphi }$，用${\mathbf{\varphi }}^{(n+1)}={\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{d}}^{(n)}$
   2.3. 计算新残差${\mathbf{r}}^{(n+1)}={\mathbf{r}}^{(n)}-{\alpha }^{(n)}\mathbf{A}{\mathbf{d}}^{(n)}$
   2.4. 计算共轭残差系数${\beta }^{(n)}=\frac{({\mathbf{r}}^{(n+1)}{)}^T{\mathbf{r}}^{(n+1)}}{({\mathbf{r}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}$
   2.5. 获取新的共轭搜索方向${\mathbf{d}}^{(n+1)}={\mathbf{r}}^{(n+1)}+{\beta }^{(n)}{\mathbf{d}}^{(n)}$

共轭梯度（CG）方法的收敛速度可以通过预处理来提高，这可以通过将原系统方程乘上一个预处理矩阵的逆矩阵${\mathbf{P}}^{-1}$来实现，其中$\mathbf{P}$是对称正定矩阵，以便做预处理${\mathbf{P}}^{-1}\mathbf{A}\mathbf{\varphi }={\mathbf{P}}^{-1}\mathbf{b}$。问题在于，即便$\mathbf{P}$和$\mathbf{A}$都是对称的，${\mathbf{P}}^{-1}\mathbf{A}$也未必是对称的，为了克服该问题，把$\mathbf{P}$采用Cholesky分解成
$$\mathbf{P}=\mathbf{L}{\mathbf{L}}^T$$
为了保证对称性，将系统方程写为
$${\mathbf{L}}^{-1}\mathbf{A}{\mathbf{L}}^{-T}{\mathbf{L}}^T\mathbf{\varphi }={\mathbf{L}}^{-1}\mathbf{b}$$
其中${\mathbf{L}}^{-1}\mathbf{A}{\mathbf{L}}^{-T}$是对称和正定的。CG方法可用于求解${\mathbf{L}}^T\mathbf{\varphi }$，从而解出$\mathbf{\varphi }$。然而，通过替换，可以把$\mathbf{L}$从方程中消去，而且不会影响到对称性和该方法的可靠性。通过对CG方法做一定的处理，便可获得预处理CG方法。

预处理CG方法（PCG）算法流程如下

1. 选择初始方向 ${\mathbf{r}}^{(0)}=\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(0)}$和${\mathbf{d}}^{(0)}={\mathbf{P}}^{-1}{\mathbf{r}}^{(0)}$
2. 第n步的迭代过程如下
   2.1. 选择$\mathbf{d}$方向的比率因子${\alpha }^{(n)}=\frac{({\mathbf{r}}^{(n)}{)}^T{\mathbf{P}}^{-1}{\mathbf{r}}^{(n)}}{({\mathbf{d}}^{(n)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}}$
   2.2. 获得新的$\mathbf{\varphi }$，用${\mathbf{\varphi }}^{(n+1)}={\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{d}}^{(n)}$
   2.3. 计算新残差${\mathbf{r}}^{(n+1)}={\mathbf{r}}^{(n)}-{\alpha }^{(n)}\mathbf{A}{\mathbf{d}}^{(n)}$
   2.4. 计算共轭残差系数${\beta }^{(n+1)}=\frac{({\mathbf{r}}^{(n+1)}{)}^T{\mathbf{P}}^{-1}{\mathbf{r}}^{(n+1)}}{({\mathbf{r}}^{(n)}{)}^T{\mathbf{P}}^{-1}{\mathbf{r}}^{(n)}}$
   2.5. 获取新的共轭搜索方向${\mathbf{d}}^{(n+1)}={\mathbf{P}}^{-1}{\mathbf{r}}^{(n+1)}+{\beta }^{(n+1)}{\mathbf{d}}^{(n)}$

实际上，有很多预处理的方法，较为简单的方法是选用原矩阵$\mathbf{A}$的对角元素构建一个对角元矩阵（Jacobi pre-conditioner），较为复杂的方法是使用不完全Cholesky分解。无论如何，当CG方法用于大型矩阵求解时一定要采用预处理，毕竟能提高收敛效率不是。

3.14 双共轭梯度（BiCG）算法及其预处理器

一般来说，扩散方程、不可压缩的压力方程离散后的格式是对称的，可以用CG方法来求解，可是CFD的更加普遍的其它方程则不是对称的，就没法直接用CG方法了。但是CG方法的收敛实在是太快了，不用又怪可惜的。所以，要是能把不对称的矩阵转化成对称的矩阵，就能用CG方法了，有一种法子是这样的
$$\left[\begin{array}{cc}0& \mathbf{A}\\ {\mathbf{A}}^T& 0\end{array}\right]\left[\begin{array}{c}\hat{\mathbf{\varphi }}\\ \mathbf{\varphi }\end{array}\right]=\left[\begin{array}{c}\mathbf{b}\\ 0\end{array}\right]$$
其中$\hat{\mathbf{\varphi }}$是虚假/伪变量，其添加只是为了让原本不对称的系统转化成对称系统，从而可以用CG方法求解。那么，CG方法现在求解了两个矢量，其中一个就是原方程$\mathbf{A}\mathbf{\varphi }=\mathbf{b}$要求解的矢量$\mathbf{\varphi }$，另一个则是不需要的伪矢量$\hat{\mathbf{\varphi }}$（如果想算就把它算出来，如果不想算就不算它了）。因为要求解两个序列的矢量，所以这个方法便称之为双共轭梯度（BiCG）。

实际上基于这个思路提了很多BiCG方法，Lanczos提出的BiCG方法流程如下

1. 选择初始方向 ${\mathbf{d}}^{(0)}={\mathbf{r}}^{(0)}={\hat{\mathbf{d}}}^{(0)}={\hat{\mathbf{r}}}^{(0)}=\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(0)}$
2. 第n步的迭代过程如下
   2.1. 选择$\mathbf{d}$方向的比率因子${\alpha }^{(n)}=\frac{({\hat{\mathbf{r}}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}{({\hat{\mathbf{d}}}^{(n)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}}$
   2.2. 获得新的$\mathbf{\varphi }$，用${\mathbf{\varphi }}^{(n+1)}={\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{d}}^{(n)}$
   2.3. 计算新残差${\mathbf{r}}^{(n+1)}={\mathbf{r}}^{(n)}-{\alpha }^{(n)}\mathbf{A}{\mathbf{d}}^{(n)}$
   2.4. 计算新伪残差${\hat{\mathbf{r}}}^{(n+1)}={\hat{\mathbf{r}}}^{(n)}-{\alpha }^{(n)}{\mathbf{A}}^T{\hat{\mathbf{d}}}^{(n)}$
   2.5. 计算共轭残差系数${\beta }^{(n+1)}=\frac{({\hat{\mathbf{r}}}^{(n+1)}{)}^T{\mathbf{r}}^{(n+1)}}{({\hat{\mathbf{r}}}^{(n)}{)}^T{\mathbf{r}}^{(n)}}$
   2.6. 获取新的搜索方向${\mathbf{d}}^{(n+1)}={\mathbf{r}}^{(n+1)}+{\beta }^{(n+1)}{\mathbf{d}}^{(n)}$
   2.7. 获取新的伪搜索方向${\hat{\mathbf{d}}}^{(n+1)}={\hat{\mathbf{r}}}^{(n+1)}+{\beta }^{(n+1)}{\hat{\mathbf{d}}}^{(n)}$

BiCG方法在每步迭代中需要对矩阵和矩阵的转置做两次相乘运算，这使得其计算量是CG方法的两倍。

BiCG方法当然也可以做预处理，跟预处理CG方法一样，Fletcher的方法流程如下（$\mathbf{P}$代表预处理矩阵）：

1. 选择初始方向 ${\mathbf{r}}^{(0)}={\hat{\mathbf{r}}}^{(0)}=\mathbf{b}-\mathbf{A}{\mathbf{\varphi }}^{(0)}$，${\mathbf{d}}^{(0)}={\mathbf{P}}^{-1}{\mathbf{r}}^{(0)}$，${\hat{\mathbf{d}}}^{(0)}={\mathbf{P}}^{-T}{\hat{\mathbf{r}}}^{(0)}$
2. 第n步的迭代过程如下
   2.1. 选择$\mathbf{d}$方向的比率因子${\alpha }^{(n)}=\frac{({\hat{\mathbf{r}}}^{(n)}{)}^T{\mathbf{P}}^{-1}{\mathbf{r}}^{(n)}}{({\hat{\mathbf{d}}}^{(n)}{)}^T\mathbf{A}{\mathbf{d}}^{(n)}}$
   2.2. 获得新的$\mathbf{\varphi }$，用${\mathbf{\varphi }}^{(n+1)}={\mathbf{\varphi }}^{(n)}+{\alpha }^{(n)}{\mathbf{d}}^{(n)}$
   2.3. 计算新残差${\mathbf{r}}^{(n+1)}={\mathbf{r}}^{(n)}-{\alpha }^{(n)}\mathbf{A}{\mathbf{d}}^{(n)}$
   2.4. 计算新伪残差${\hat{\mathbf{r}}}^{(n+1)}={\hat{\mathbf{r}}}^{(n)}-{\alpha }^{(n)}{\mathbf{A}}^T{\hat{\mathbf{d}}}^{(n)}$
   2.5. 计算共轭残差系数${\beta }^{(n+1)}=\frac{({\hat{\mathbf{r}}}^{(n+1)}{)}^T{\mathbf{P}}^{-1}{\mathbf{r}}^{(n+1)}}{({\hat{\mathbf{r}}}^{(n)}{)}^T{\mathbf{P}}^{-1}{\mathbf{r}}^{(n)}}$
   2.6. 获取新的搜索方向${\mathbf{d}}^{(n+1)}={\mathbf{P}}^{-1}{\mathbf{r}}^{(n+1)}+{\beta }^{(n+1)}{\mathbf{d}}^{(n)}$
   2.7. 获取新的伪搜索方向${\hat{\mathbf{d}}}^{(n+1)}={\mathbf{P}}^{-1}{\hat{\mathbf{r}}}^{(n+1)}+{\beta }^{(n+1)}{\hat{\mathbf{d}}}^{(n)}$

BiCG方法的其它变种，更具稳定性和可靠性的，有Sonneveld的共轭梯度平方（CGS）、Van Der Vorst的双共轭梯度稳定（Bi-CGSTAB）、通用最小残差方法（GMRES），等。这些方法对于求解CFD当中出现的大型系统矩阵非常有用，可应用于结构和非结构网格所产生的非对称矩阵。

4 多重网格方法

迭代方法的收敛速度随着代数方程系统维数的增大而急剧恶化，这一情况即便是对于中等程度维数的系统也是存在的。这使得迭代方法的应用很受限制，幸运的是，人们很快发现把多重网格算法和迭代方法联合起来使用可以很好地克服该缺陷。

多重网格方法的提出源自于这些学者的贡献：Fedorenko（几何多重网格Geometric Multigrid）、Poussin（代数多重网格Algebraic Multigrid）、以及Settari和Azziz、还有Brandt在理论上的工作。尽管高频或振荡误差可以很容易地通过标准迭代方法（Jacobi、Gauss-Seidel、ILU）消除，然而这些求解方法无法轻易地去除光滑或低频误差分量。正因为如此，这些解法在多重网格方法中被定义成光滑器，下图展示了一维问题的误差频率模态的变化情况。

从上向下看，图中展示的误差模态从高频到低频，波长也从短${\lambda }_1$到长${\lambda }_5$，一维计算域采用一维网格离散，不同的模态分别绘制在同样的网格上。可见，高频误差的振荡发生在一个单元内部，其可以轻易地被迭代方法感受到。但是伴随着误差频率的降低，即，波长的增大，误差在网格上变得越来越光滑，在一个单元内部仅可感受到误差波长的一小部分，网格越是加密该问题反而越是严重，使得系统维数的增大导致方程数目的增多和收敛速率的恶化。

多重网格方法提高迭代求解器效率的法子是，确保在一个网格层级（level）上由于光滑器的应用所产生的低频误差，可以转化到，较粗糙网格层级上的高频误差。如图，通过使用多层次的粗细网格，多重网格方法能够克服收敛恶化的现象。

一般而言，粗网格可以用细网格的拓扑和几何关系来生成，这相当于在细网格的上面重新生成了一层粗网格；粗网格也可以直接通过细网格单元的聚合来生成，这样就是所谓的代数多重网格（Algebraic MutiGrid Method （AMG））。在AMG中，并不直接用到几何信息，聚合过程只是纯粹的代数处理，即通过聚合过程从细网格上的方程直接构建粗网格上的方程（并不需要重画网格啥的）。该方法可用于构造高效可靠的线性求解器，适用于高度各向异性网格或方程中系数变化很大的问题。

不管哪种方法，多重网格循环过程将用于遍历不同的网格层级。每次从细网格到粗网格的遍历包含：（i）限制（restriction）过程，（ii）设置或更新系统方程到粗网格，（iii）应用几次光滑迭代。每个从粗网格到细网格的遍历包含：（i）扩展（prolongation）过程，（ii）在细网格层级上的流场值修正，（iii）在限制过程中构造的方程上应用几次光滑迭代。下面将详细讲讲不同的步骤。

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

求解过程的首要任务是通过单元聚合/粗化来生成粗糙/细腻网格层级，有三种方法可以用来达成此目的。一，首先生成粗网格，然后再加密生成细网格，优点是粗-细网格关系非常好获取，缺点是细网格的分布要依赖于粗网格。二，使用非嵌套网格使得网格层级间的信息传递变得耗费较大（这个我确实没看明白，只好直译了）。这两种方法都不允许复杂区域有很好的分辨率。在第三种方法（推荐使用）中，由最精细的网格开始求解，然后通过把细网格单元做聚合处理来生成粗网格，如图所示，聚合的过程或者由单元的几何拓扑关系给出，或者是由与邻近单元的系数满足一定的判据来给出。以下的讨论都针对第三种方法。

粗网格是通过聚合算法把细网格单元汇聚起来生成的。对每个粗网格层级，该算法不断进行直到细网格层级的所有单元都与粗网格层级的单元相互关联为止。在该聚合过程中，细网格点被逐一访问，若某个单元被选中为种子（seed）单元，则其周围满足设定判据的邻居单元将会与该种子单元聚合形成一个粗单元。提前需要设定好融合成一个粗单元的细单元的最大数目。如果所选的种子单元没办法生成一个粗单元，那么把它加到其邻居单元中的最小粗单元中即可。

一个高效的聚合算法是直接聚合（DA）算法，由Mavriplis提出。在DA中，聚合由种子单元开始，将其与邻近的细单元，基于它们之间的几何关联强度做聚合。该过程仅需在求解开始前执行一次就好。

4.2 限制步和粗层次系数（The Restriction Step and Coarse Level Coefficients）

先在细网格上求解，执行几次迭代后，将误差转化或者说是限制到粗网格层级上，然后在粗网格上迭代几次，再往更高层级的网格上限制和求解，如此重复，直至达到最高或者说是最粗糙的网格层级。用(k)代表不同层级，则通过求解如下的系统方程可获得该层级的解
$${\mathbf{A}}^{(k)}{\mathbf{e}}^{(k)}={\mathbf{r}}^{(k)}$$
下一个粗糙层级是(k+1)，误差要限制到(k+1)层。令$G_I$代表聚合成粗网格层级(k+1)上的I单元的细网格层级(k)上的i单元的集合，那么，在(k+1)层级要求解的系统是
$${\mathbf{A}}^{(k+1)}{\mathbf{e}}^{(k+1)}={\mathbf{r}}^{(k+1)}$$
上式右端项的残差为
$${\mathbf{r}}^{(k+1)}={\mathbf{I}}_k^{k+1}{\mathbf{r}}^{(k)}$$
其中${\mathbf{I}}_k^{k+1}$为限制算子（即，插值矩阵），其反映的是细网格到粗网格的聚合过程中残差的传递或者说转化特性。在AMG中，限制算子直接用线性形式定义，即等于细网格残差的加和
$${\mathbf{r}}_I^{(k+1)}=\sum _{i\in G_I}{\mathbf{r}}_i^{(k)}$$
此外，粗网格单元的系数是由细网格单元的系数相加来构造的，再回顾下离散后所得的线性方程形式
$$a_C{\varphi }_C+\sum _{F\sim NB(C)}{a}_F{\varphi }_F=b_C$$
对于(k)层级，写成
$$a_i^{(k)}{\varphi }_i^{(k)}+\sum _{j\sim NB(i)}{a}_{ij}^{(k)}{\varphi }_j^{(k)}=b_i^{(k)}$$
其中$NB(i)$代表单元$i$的邻近单元。当然一开始上式是不会满足的，会有如下的残差存在
$$r_i^{(k)}=b_i^{(k)}-\left(a_i^{(k)}{\varphi }_i^{(k)}+\sum _{j\sim NB(i)}{a}_{ij}^{(k)}{\varphi }_j^{(k)}\right)$$
定义${\varphi }_I^{(k+1)}$代表粗网格单元$I$的解，且是要传递到细网格上去的，那么由粗网格给予细网格的修正可以写成
$${\varphi }_i^{\prime (k)}={\varphi }_I^{(k+1)}-{\varphi }_i^{(k)}$$
当然我们希望修正能使得在粗网格的单元$I$上的残差是0。将新残差用${\tilde{r}}_i^{(k)}$表示，其等于
$${\tilde{r}}_i^{(k)}=b_i^{(k)}-\left(a_i^{(k)}\left({\varphi }_i^{(k)}+{\varphi }_i^{\prime (k)}\right)+\sum _{j\sim NB(i)}{a}_{ij}^{(k)}\left({\varphi }_j^{(k)}+{\varphi }_j^{\prime (k)}\right)\right)$$
等效于
$$\begin{array}{rl}{\tilde{r}}_i^{(k)}& =\underset{r_i^{(k)}}{\underbrace{b_i^{(k)}-\left(a_i^{(k)}{\varphi }_i^{(k)}+\sum _{j\sim NB(i)}{a}_{ij}^{(k)}{\varphi }_j^{(k)}\right)}}-\left(a_i^{(k)}{\varphi }_i^{\prime (k)}+\sum _{j\sim NB(i)}{a}_{ij}^{(k)}{\varphi }_j^{\prime (k)}\right)\\ & =r_i^{(k)}-\left(a_i^{(k)}{\varphi }_i^{\prime (k)}+\sum _{j\sim NB(i)}{a}_{ij}^{(k)}{\varphi }_j^{\prime (k)}\right)\end{array}$$
强迫$I$的残差加和为0，即
$$\sum _{i\in G_I}{\tilde{\mathbf{r}}}_i^{(k)}=0$$
得到
$$0=\sum _{i\in G_I}{\tilde{\mathbf{r}}}_i^{(k)}-\left(\sum _{i\in G_I}{a}_i^{(k)}{\varphi }_i^{\prime (k)}+\sum _{i\in G_I}\sum _{j\sim NB(i)}{a}_{ij}^{(k)}{\varphi }_j^{\prime (k)}\right)$$
将上式重写成粗网格形式，即，粗网格的修正方程
$$a_I^{(k+1)}{\varphi }_I^{\prime (k+1)}+\sum _{J\sim NB(I)}{a}_{IJ}^{(k+1)}{\varphi }_J^{\prime (k+1)}=r_I^{(k+1)}$$
其中$a_I^{(k+1)}$、$a_{IJ}^{(k+1)}$、$r_I^{(k+1)}$是直接从细网格系数推导出来的，即
$$\begin{array}{rl}& 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 _{i\in G_I}\sum _{j\notin G_I,j\in NB(I)}{a}_{ij}^{(k)}\\ & r_I^{(k+1)}=\sum _{i\in G_I}{r}_i^{(k)}\end{array}$$

4.3 扩展步和细网格层级修正（The Prolongation Step and Fine Grid Level Corrections）

扩展算子用于把修正从粗网格转化到细网格层级上，有许多可用的算子，其中一个如图所示，为0阶扩展算子，其在细网格上产生同样的误差值，即，粗网格单元的误差将会被其在细网格层级上该单元的子单元所继承。

修正是从粗网格上的系统方程获取的，然后再插值到细网格层级上的
$${\mathbf{e}}^{(k)}={\mathbf{I}}_{k+1}^k{\mathbf{e}}^{(k+1)}$$
其中${\mathbf{I}}_{k+1}^k$为从粗网格到细网格的插值矩阵，最终，细网格上的解得以修正
$${\varphi }^{(k)}\leftarrow {\varphi }^{(k)}+{\mathbf{e}}^{(k)}$$
网格层级的层数取决于网格的多寡，即便是多层的网格系统，其流程也和图中所展示的差不太多。

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

遍历策略指的是在求解过程中粗网格的访问方式，也被称之为多重网格循环。通常在AMG方法中用的循环是V循环、W循环、F循环，如图

AMG循环中最简单的是V循环，其对每层次的网格访问一次，对于刚性特别大的系统，V循环可能无法加速求解过程。因此，需要在粗网格上做较多求解，W循环实际上是在粗网格层级上做了些小的V循环，这样一来，W循环由嵌套的粗细网格层级扫描构成，随着AMG层数的增加其复杂性也有所增加。F循环是W循环的变种，其可视为把W循环劈开一半，F循环的粗网格扫描没W循环多，但是比V循环要多，因此，其介于V和W循环之间。

感觉书中对多重网格算法的讲解不是很透彻，看完也没能理解该如何应用，如何求解具体问题，所以准备找些材料弄明白后再单独讲讲这部分内容。