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

2020年06月10日 阅读数：1391

Chapter 10 Solving the System of Algebraic Equationshtml

1 引言

$\bold A \boldsymbol \phi = \bold b$

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

2 直接或Gauss消元方法

2.1 Gauss消元

$a_{11}\phi_1+a_{12}\phi_2=b_1 \\ ~\\ a_{21}\phi_1+a_{22}\phi_2=b_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=\frac {b_2-\displaystyle\frac{a_{21}}{a_{11}}b_1} {a_{22}-\displaystyle\frac{a_{21}}{a_{11}}a_{12}}$

$\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.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
}


2.6 LU分解

$\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 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$

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$

$l_{i1}u_{11}=a_{i1}\Rightarrow l_{i1}=\frac{a_{i1}}{u_{11}}~~~~i=2,3,...,N$

$l_{21}u_{1j}+u_{2j}=a_{2j} \Rightarrow u_{2j}=a_{2j} - l_{21}u_{1j}~~~~j=2,3,...,N$

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

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

$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 回代步

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

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

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

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


$\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 \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 \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_i = P_i \phi_{i+1} + Q_i,\quad\quad i=1,2,3,...,N$

$\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=-\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$

$P_1=-\frac{b_1}{a_1}\quad\quad Q_1=\frac{d_1}{a_2}$

$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）

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

$\phi_1=-\frac{b_1}{a_1}\phi_3-\frac{c_1}{a_1}\phi_2+\frac{f_1}{a_1}$

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

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

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

\begin{aligned} & \phi_N=R_N \\ & \phi_{N-1}=Q_{N-1}\phi_N+R_{N-1} \end{aligned}

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 迭代方法

$\bold A = \bold D + \bold L + \bold U$

$\bold A = \bold M - \bold N$

$(\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,...$

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=\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$

$\lim_{n\rightarrow\infin}\bold B^n=\lim_{n\rightarrow\infin} \underbrace{\bold B *\bold B *\bold B*... *\bold B}_{n~times}=\bold 0$

$\rho(\bold B) < 1$