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（健壯的，魯棒的）的。

關於該條件的深層理解，能夠經過定義精確值與迭代值之間的錯誤矢量