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

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

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

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

1 引言

线性求解器求解的是以下形式的代数方程组
$\bold A \boldsymbol \phi = \bold b$
其中 $\bold A$是单元的系数矩阵 $a_{ij}$， $\boldsymbol \phi$是未知变量 $\phi$的向量，而 $\bold b$则是源项 $b_i$的矢量。使用矩阵编号，该方程的展开形式为
$\begin{bmatrix} a_{11} & a_{12} & ... & a_{1~N-1} & a_{1~N} \\ a_{21} & a_{22} & ... & a_{2~N-1} & a_{2~N} \\ ... & ... & ... & ... & ... \\ a_{N1} & a_{N2} & ... & a_{N~N-1} & a_{N~N} \end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \\ ... \\ \phi_{N-1} \\ \phi_N \end{bmatrix} = \begin{bmatrix} b_1 \\ b_2 \\ ... \\ b_{N-1} \\ b_N \end{bmatrix}$
通常来讲呢，矩阵 $\bold A$中的每一行表明计算域1个单元上定义的方程，非0系数是与该单元紧邻的邻居单元的，系数 $a_{ij}$反映的是存储在 $i$单元控制体形心的 $\phi_i$与其邻近单元形心的 $\phi_j$的关系紧密程度。因为一个单元只和少数几个单元相邻，且其邻接关系是与离散区域上单元的剖分编码方式相关的，因此大部分的系数其实是0，这就使得矩阵 $\bold A$老是个稀疏矩阵（即，非零元素的个数不多，在总体矩阵中占比很小，绝大多数的元素都是0）。更进一步，若是用的是结构化的网格，那么矩阵 $\bold A$将会是带状矩阵，即，其仅含几条非零的对角线元素。所以，求解该类系统时，可针对其特性选用对应的高效方法。

如上所述，求解代数方程组系统的方法能够大致上分为两大类，即直接和迭代方法。在直接方法中，矩阵 $\bold A$求逆，解 $\bold\phi$经过 $\bold\phi=\bold A^{-1}\bold b$来求出。当矩阵 $\bold A$很是大的时候，应用直接线性解法来求解CFD问题是不切实际的，由于这些CFD问题一般包含非线性系统方程， 即它们的系数是与解相关联的，那么使用迭代方法是切合实际的。

另外一方面，在迭代代数解法中，求解算法将重复屡次直至达到预期的收敛水平，而不须要在每次迭代的过程当中都得到彻底收敛的解。

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

2 直接或Gauss消元方法

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

2.1 Gauss消元

假设有以下2变量 $\phi_1$和 $\phi_2$的线性方程组
$a_{11}\phi_1+a_{12}\phi_2=b_1 \\ ~\\ a_{21}\phi_1+a_{22}\phi_2=b_2$
为了消去 $\phi_1$，将第1个式子乘上 $a_{21}/a_{11}$，并减去第2个式子，可得
$\left(a_{22}-\frac{a_{21}}{a_{11}}a_{12}\right)\phi_2=b_2-\frac{a_{21}}{a_{11}}b_1$
从而能够直接求得 $\phi_2$
$\phi_2=\frac {b_2-\displaystyle\frac{a_{21}}{a_{11}}b_1} {a_{22}-\displaystyle\frac{a_{21}}{a_{11}}a_{12}}$
再把这个求得的 $\phi_2$代入到方程最初的第1个式子里，能够求得 $\phi_1$的值
$\phi_1=\frac{b_1}{a_{11}}-\frac{a_{12}}{a_{11}}\frac {b_2-\displaystyle\frac{a_{21}}{a_{11}}b_1} {a_{22}-\displaystyle\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分解执行以后，对于右端项 $\bold b$不一样的线性系统就想求解多少次就求解多少次，而不须要再作消元处理了（至关于把第1次的消元处理作了LU分解），然而在Gauss消元方法中，消元是始终须要进行的。

基于前面消元处理过程，是将本来的矩阵 $\bold A$转化为上三角矩阵，即
$\begin{bmatrix} u_{11} & u_{12} & u_{13} & ... & u_{1~N-1} & u_{1~N} \\ 0 & u_{22} & u_{23} &... & u_{2~N-1} & u_{2~N} \\ 0 & 0 & u_{33} & ... & u_{3~N-1} & u_{3~N} \\ ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & ... & 0 & u_{N~N} \end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \\ \phi_3 \\ ... \\ \phi_N \end{bmatrix} = \begin{bmatrix} c_1 \\ c_2 \\ c_3 \\ ... \\ c_N \end{bmatrix}$
使用缩写形式，即
$\bold U \boldsymbol \phi - \bold c = \bold 0$
令 $\bold L$为单位下三角矩阵（对角线上的元素是1），即
$\bold L = \begin{bmatrix} 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_{N~N-1} & 1 \end{bmatrix}$
将方程 $\bold U \boldsymbol \phi - \bold c = \bold 0$左乘以 $\bold L$，则变回了最初的方程
$\bold L(\bold U \boldsymbol \phi - \bold c) = \bold L \bold U \boldsymbol \phi - \bold L \bold c=\bold A \boldsymbol \phi - \bold b$
这样，变获得
$\bold L \bold U = \bold A$
和
$\bold L \bold c = \bold b$
即，矩阵 $\bold A$写成了下三角和上三角矩阵的乘积形式，即 $LU$分解。

2.7 分解步骤

将 $\bold L \bold U = \bold A$写成展开形式
$\begin{bmatrix} 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_{N~N-1} & 1 \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} & ... & u_{1~N-1} & u_{1~N} \\ 0 & u_{22} & u_{23} &... & u_{2~N-1} & u_{2~N} \\ 0 & 0 & u_{33} & ... & u_{3~N-1} & u_{3~N} \\ ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & ... & 0 & u_{N~N} \end{bmatrix} = \begin{bmatrix} a_{11} & a_{12} & a_{13} & ... & a_{1~N-1} & a_{1~N} \\ a_{21} & a_{22} & a_{23} &... & a_{2~N-1} & a_{2~N} \\ a_{31} & a_{32} & a_{33} &... & a_{3~N-1} & a_{3~N} \\ ... & ... & ... & ... & ... & ... \\ a_{N1} & a_{N2} & a_{N3} & ... & a_{N~N-1} & a_{N~N} \end{bmatrix}$
$\bold L$的第1行和 $\bold U$的每一列（第1到N列）相乘，获得
$u_{1j}=a_{1j}~~~~j=1,2,3,...,N$
接下来，让 $\bold L$的第2到第N行与 $\bold U$的第1列相乘，获得
$l_{i1}u_{11}=a_{i1}\Rightarrow l_{i1}=\frac{a_{i1}}{u_{11}}~~~~i=2,3,...,N$
重复该过程，让 $\bold L$的第2行与 $\bold U$的第2到第N列相乘，获得
$l_{21}u_{1j}+u_{2j}=a_{2j} \Rightarrow u_{2j}=a_{2j} - l_{21}u_{1j}~~~~j=2,3,...,N$
而后呢，再把 $\bold L$的第3到N行与 $\bold U$的第2列相乘，获得
$l_{i1}u_{12}+l_{i2}u_{22}=a{i2} \Rightarrow l_{i2}=\frac{a{i2}-l_{i1}u_{12}}{u_{22}}~~~~i=3,4,...,N$
如此不断重复，直至解出全部 $\bold L$和 $\bold U$的元素值。

通常来讲， $\bold L$的第i行与 $\bold U$的第i到第N列相乘，可得
$u_{ij}=a_{ij}-\sum_{k=1}^{i-1}l_{ik}u_{kj}~~~~j=i,i+1,...,N$
而 $\bold L$的第i+1到第N行与 $\bold U$的第i列相乘，则可得
$l_{ki}=\frac{a_{ki}-\sum_{j=1}^{i-1}l_{kj}u_{ji}}{u_{ii}}~~~~k=i+1,i+2,...,N$
对于 $\bold L$的第N行，其系数与 $\bold 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 回代步

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

第1步是求解 $\bold L \bold c = \bold b$，经过前向回代，能够直接得出
$\begin{aligned} & c_1 = b_1 \\ & c_i = b_i - \sum_{j=1}^{i-1}l_{ij}c_j,\quad\quad i=2,3,...,N \end{aligned}$

第2步是求解 $\bold U \boldsymbol \phi = \bold c$，采用反向回代，能够直接获得
$\begin{aligned} & \phi_N = \frac{C_N}{u_{NN}} \\ & \phi_i = \frac{c_i - \displaystyle \sum_{j=i+1}^{N}u_{ij}\phi_j}{u_{ii}},\quad\quad i=N-1,N-2,...,3,2,1 \end{aligned}$

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

2.10 LU分解和Gauss消元

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

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

若是把 $\bold L$和 $\bold U$放在原来的矩阵 $\bold A$中，那么直接这么作就行了，至关于在Gauss消元的算法中，把系数存放到了下三角位置上，注意，不须要对 $\bold b$作处理，由于LU分解就是为了求解不一样的右端项 $\bold 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
	}
}

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

---

例1 用LU分解方法求解以下线性代数方程组
$\begin{bmatrix} 3 & -1 & 0 & 0 \\ -2 & 6 & -1 & 0 \\ 0 & -2 & 6 & -1 \\ 0 & 0 & -2 & 7 \\ \end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \\ \phi_3 \\ \phi_4 \end{bmatrix} =\begin{bmatrix} 3 \\ 4 \\ 5 \\ -3 \end{bmatrix}$

解
求解过程从略，答案直接给出
$\bold L=\begin{bmatrix} 1 & 0 & 0 & 0 \\ -\displaystyle\frac{2}{3} & 1 & 0 & 0 \\ 0 & -\displaystyle\frac{3}{8} & 1 & 0 \\ 0 & 0 & -\displaystyle\frac{16}{45} & 1 \\ \end{bmatrix}, \quad\quad \bold U =\begin{bmatrix} 3 & -1 & 0 & 0 \\ 0 & \displaystyle\frac{16}{3} & -1 & 0 \\ 0 & 0 & \displaystyle\frac{45}{8} & -1 \\ 0 & 0 & 0 & \displaystyle\frac{299}{45} \\ \end{bmatrix}$
$\bold c=\begin{bmatrix} 3 \\ 6 \\ 29/4\\ -19/45 \end{bmatrix},\quad\quad \boldsymbol \phi=\begin{bmatrix} 435/299 \\ 408/299 \\ 382/299 \\ -19/299 \end{bmatrix}$

---

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

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

2.13 三对角矩阵算法（TDMA）

一维问题导出的方程每每具备三对角形式，即
$a_i \phi_i + b_i \phi_{i+1}+c_i \phi_{i-1} = d_i, \quad\quad i=1,2,3,...,N, \quad\quad c_1=b_N=0$
对于 $i=1$，可直接推得 $\phi_1$和 $\phi_2$的关系
$i=1 \Rightarrow a_1\phi_1=-b_1\phi_2 + d_1 \Rightarrow \phi_1=-\frac{b_1}{a_1}\phi_2+\frac{d_1}{a_2}$
一样对于 $i=2$，可导出 $\phi_2$和 $\phi_3$的关系
$i=2 \Rightarrow a_2\phi_2 = -b_2\phi_3 - c_2\phi_1 + d_2 \Rightarrow \phi_2=-\frac{a_1b_2}{a_1a_2-c_2b_1}\phi_3 + \frac{d_2a_1-c_2d_1}{a_1a_2-c_2b_1}$
一样的操做可用于 $\phi_3$到 $\phi_N$，假设 $\phi_i$可用 $\phi_{i+1}$表示为
$\phi_i = P_i \phi_{i+1} + Q_i,\quad\quad i=1,2,3,...,N$
假设对 $i-1$上式已知，那么对 $i$的推导以下
$\left. \begin{matrix} \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 \end{matrix} \right\} \Rightarrow \phi_i=-\frac{b_i}{a_i+c_iP_{i-1}}\phi_{i+1}+\frac{d_i-c_iQ_{i-1}}{a_i+c_iP_{i-1}}$
从而推得 $P_i$和 $Q_i$为
$P_i=-\frac{b_i}{a_i+c_iP_{i-1}} \quad\quad Q_i=\frac{d_i-c_iQ_{i-1}}{a_i+c_iP_{i-1}} \quad\quad i=2,3,4,...,N$
对于 $i=1$，可直接算得 $P_1$和 $Q_1$
$P_1=-\frac{b_1}{a_1}\quad\quad Q_1=\frac{d_1}{a_2}$
而对于 $i=N$，因为 $b_N=0$，有
$b_N=0\Rightarrow P_N=0 \Rightarrow \phi_N=Q_N$
TDMA算法总结以下

1. 计算 $P_1$和 $Q_1$，用
   $P_1=-\frac{b_1}{a_1}\quad\quad Q_1=\frac{d_1}{a_2}$
2. 计算 $P_i$和 $Q_i$， $i=2,3,...,N$，用
   $P_i=-\frac{b_i}{a_i+c_iP_{i-1}} \quad\quad Q_i=\frac{d_i-c_iQ_{i-1}}{a_i+c_iP_{i-1}}$
3. 求出 $\phi_N$，用
   $\phi_N=Q_N$
4. 反向依次求出 $\phi_i$， $i=N-1,N-2,...,3,2,1$，用
   $\phi_i = P_i \phi_{i+1} + Q_i,\quad\quad i=N-1,N-2,...,3,2,1$

2.14 五对角矩阵算法（PDMA）

仍是1维问题，只不过离散格式用的精度更高，让i节点与i+1，i+2，i-1，i-2都有关系，这样子出来的离散矩阵即是5对角形式了，即
$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\quad i=1,2,3,...,N$
对于边界的几个节点，有
$d_1=e_1=e_2=0 \\ b_{N-1}=b_{N}=c_{N}=0$
对于 $i=1$，有
$\phi_1=-\frac{b_1}{a_1}\phi_3-\frac{c_1}{a_1}\phi_2+\frac{f_1}{a_1}$
对于 $i=2$，有
$\phi_2=-\frac{a_1b_2}{a_1a_2-d_2c_1}\phi_4-\frac{a_1c_2-b_1d_2}{a_1a_2-d_2c_1}\phi_3+\frac{a_1f_2-d_2f_1}{a_1a_2-d_2c_1}$
该过程可重复下去，获取第 $i$个变量 $\phi_i$的通有形式
$\phi_i=P_i\phi_{i+2}+Q_i\phi_{i+1}+R_i\quad\quad i=1,2,3,...,N$
在 $\phi_{i-1}$和 $\phi_{i-2}$已知的状况下，可求得 $\phi_i$的形式为
$\begin{aligned} \phi_i=&-\frac{b_i}{a_i+e_iP_{i-2}+(d_i+e_iQ_{i-2})Q_{i-1}}\phi_{i+2} \\ &-\frac{c_i+(d_i+e_iQ_{i-2})P_{i-1}}{a_i+e_iP_{i-2}+(d_i+e_iQ_{i-2})Q_{i-1}}\phi_{i+1} \\ &-\frac{f_i-e_iR_{i-2}-(d_i+e_iQ_{i-2})R_{i-1}}{a_i+e_iP_{i-2}+(d_i+e_iQ_{i-2})Q_{i-1}} \end{aligned}$
比较可得
$P_i=-\frac{b_i}{a_i+e_iP_{i-2}+(d_i+e_iQ_{i-2})Q_{i-1}} \\ Q_i= -\frac{c_i+(d_i+e_iQ_{i-2})P_{i-1}}{a_i+e_iP_{i-2}+(d_i+e_iQ_{i-2})Q_{i-1}} \\ R_i= -\frac{f_i-e_iR_{i-2}-(d_i+e_iQ_{i-2})R_{i-1}}{a_i+e_iP_{i-2}+(d_i+e_iQ_{i-2})Q_{i-1}}$
前两个节点，即 $i=1,2$的值为
$\begin{aligned} & 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_2Q_1}, \quad Q_2 = -\frac{c_2+d_2P_1}{a_2+d_2Q_1}, \quad R_2 = \frac{f_2-d_2R_1}{a_2+d_2Q_1} \end{aligned}$
因为 $b_{N-1}=b_{N}=c_{N}=0$，因此 $P_{N-1}=P_N=Q_N=0$，那么 $\phi_N$和 $\phi_{N-1}$的方程为
$\begin{aligned} & \phi_N=R_N \\ & \phi_{N-1}=Q_{N-1}\phi_N+R_{N-1} \end{aligned}$
将PDMA算法的流程汇总以下

1. 计算 $P_1,~Q_1,~R_1,~P_2,~Q_2,~R_2$，使用
   $\begin{aligned} & 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_2Q_1}, \quad Q_2 = -\frac{c_2+d_2P_1}{a_2+d_2Q_1}, \quad R_2 = \frac{f_2-d_2R_1}{a_2+d_2Q_1} \end{aligned}$
2. 对于 $i=3,4,5,...,N$计算 $P_i, ~Q_i,~ R_i$，用
   $P_i=-\frac{b_i}{a_i+e_iP_{i-2}+(d_i+e_iQ_{i-2})Q_{i-1}} \\ Q_i= -\frac{c_i+(d_i+e_iQ_{i-2})P_{i-1}}{a_i+e_iP_{i-2}+(d_i+e_iQ_{i-2})Q_{i-1}} \\ R_i= -\frac{f_i-e_iR_{i-2}-(d_i+e_iQ_{i-2})R_{i-1}}{a_i+e_iP_{i-2}+(d_i+e_iQ_{i-2})Q_{i-1}}$
3. 计算 $\phi_N$和 $\phi_{N-1}$，用
   $\begin{aligned} & \phi_N=R_N \\ & \phi_{N-1}=Q_{N-1}\phi_N+R_{N-1} \end{aligned}$
4. 对于 $i=N-2,N-3,...,3,2,1$，计算 $\phi_i$，用
   $\phi_i=P_i\phi_{i+2}+Q_i\phi_{i+1}+R_i$

3 迭代方法

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

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

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

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

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

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

在详细描述各类迭代方法以前，先讲下迭代方法须要知足什么样的一些特性，才能让其获得收敛的解。

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

B. 从 $\boldsymbol \phi^{0} \neq \boldsymbol \phi$开始，该方法应保证 $\boldsymbol \phi^{n}$随着迭代次数 $n$的增长最终收敛到 $\boldsymbol \phi$。由于 $\boldsymbol \phi^{n}$能够展开成 $\boldsymbol \phi^{0}$的以下形式
$\boldsymbol \phi^{n} = \bold B^n \boldsymbol \phi^{0} + \sum_{i=0}^{n-1}\bold B^i \bold C \bold b$
那么，若要收敛成立，则须要让 $\bold B$知足（即让 $\boldsymbol \phi^{0}$的影响彻底消失掉）
$\lim_{n\rightarrow\infin}\bold B^n=\lim_{n\rightarrow\infin} \underbrace{\bold B *\bold B *\bold B*... *\bold B}_{n~times}=\bold 0$
这意味着 $\bold B$的谱半径要小于1，即
$\rho(\bold B) < 1$
该条件保证了迭代方法是有自我纠正能力的，即，其对于精确解的任何不利的错误扰动是robust（健壮的，鲁棒的）的。

关于该条件的深层理解，能够经过定义精确值与迭代值之间的错误矢量