1. 問(wèn)題描述

本算例來(lái)自B站Up主“Red-Green鯉魚(yú)”的系列教程。本文主要介紹計(jì)算代數(shù)方程組的三種點(diǎn)迭代方法。

1.1. 泊松方程

含有二階偏導(dǎo)數(shù)的偏微分方程：
$\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}=f(x,y)$

當(dāng) $f=0$ 時(shí)，上述方程被稱為拉普拉斯方程。許多物理過(guò)程都可以用泊松方程來(lái)描述，如熱傳導(dǎo)方程
$\frac{\partial \phi}{\partial t}=\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}$

在求解不可壓縮流動(dòng)的NS方程時(shí)，通常將已知壓力場(chǎng)代入動(dòng)量方程來(lái)預(yù)估速度場(chǎng)，然后將預(yù)估的速度場(chǎng)代入連續(xù)性方程中，由于預(yù)估的速度場(chǎng)中包含壓力梯度項(xiàng)，因此代入連續(xù)性方程后會(huì)得到 $\nabla \cdot(\nabla P)$ ，即壓力泊松方程

1.2. 算例

令上述泊松方程的源項(xiàng)為如下形式：
$f(x,y)=-4\cdot sin(x-y)\cdot e^{(x-y)}$

其解析解為
$\phi(x,y)=cos(x-y)\cdot e^{(x-y)}$

比較數(shù)值解與解析解在定義域 $x\in [-1,1], y\in [-1,1]$ 上的差別。邊界條件為Dirichlet，邊界上的值由解析解求出。

2. 區(qū)域離散和方程離散

$x$ 方向設(shè)置 $N$ 個(gè)結(jié)點(diǎn)，編號(hào)從 $1-N$ ， $y$ 方向上設(shè)置 $M$ 個(gè)結(jié)點(diǎn)，編號(hào)從 $1-M$ ，結(jié)點(diǎn)間距分別為 $\Delta x$ 和 $\Delta y$ ，將 $(i,j)$ 號(hào)結(jié)點(diǎn)記為 $P$ ，該結(jié)點(diǎn)上下左右四個(gè)結(jié)點(diǎn)分別記為 $N,S,W,E$

采用有限差分法計(jì)算泊松方程，二階偏導(dǎo)數(shù)項(xiàng)采用二階中心差分離散，
$\frac{\partial^2 \phi}{\partial x^2}\bigg|_P=\frac{\phi_W+\phi_E-2\phi_P}{\Delta x^2}$

$\frac{\partial^2 \phi}{\partial y^2}\bigg|_P=\frac{\phi_N+\phi_S-2\phi_P}{\Delta y^2}$

$\frac{\phi_W+\phi_E-2\phi_P}{\Delta x^2}+\frac{\phi_N+\phi_S-2\phi_P}{\Delta y^2}=f_P$

$\phi_P=\frac{(\phi_W+\phi_E)\Delta y^2+(\phi_N+\phi_S)\Delta x^2-f_P \Delta x^2 \Delta y^2}{2(\Delta x^2+\Delta y^2)}$

令 $\Delta x=\Delta y=h$ ，則上式化為
$\phi_P=\frac{1}{4}(\phi_W+\phi_E+\phi_N+\phi_S-h^2f_P)$

2.1. 邊界條件

左邊界
$\phi_{1,j}=cos(-1-y_j)\cdot e^{(-1-y_j)}$

右邊界
$\phi_{N,j}=cos(1-y_j)\cdot e^{(1-y_j)}$

下邊界
$\phi_{i,1}=cos(x_i+1)\cdot e^{(x_i+1)}$

上邊界
$\phi_{i,N}=cos(x_i-1)\cdot e^{(x_i-1)}$

3. 代數(shù)方程組求解

上一篇文章介紹了代數(shù)方程組求解方法中的直接解法TDMA，本文介紹另一大類求解方法：迭代法。迭代法的思想也可以概況為“預(yù)測(cè)-校正”，給出初始值，通過(guò)不斷迭代逐步改進(jìn)，直到達(dá)到一定精度要求為止。
該方法需要首先構(gòu)造迭代方式；其次是所構(gòu)造的迭代序列是否收斂，如果收斂則要進(jìn)一步提高收斂速度。

迭代法可分為點(diǎn)迭代、塊迭代、交替方向迭代法以及強(qiáng)隱迭代法。在點(diǎn)迭代法中，每一步計(jì)算只能改進(jìn)求解區(qū)域中一個(gè)結(jié)點(diǎn)的值，且該值是由一個(gè)顯函數(shù)形式由其余各點(diǎn)的已知值來(lái)確定，因而點(diǎn)迭代法又稱為顯式迭代法。

下面將討論點(diǎn)迭代法的三種實(shí)施方式。

3.1. 雅可比迭代

$\phi_P^{n+1}=\frac{1}{4}(\phi_W^n+\phi_E^n+\phi_N^n+\phi_S^n-h^2f_P)$
上式中上標(biāo) $n$ 為當(dāng)前預(yù)測(cè)值， $n+1$ 為代入迭代方程后的校正值

3.2. 高斯-賽德?tīng)柕?/h2>
$\phi_P^{n+1}=\frac{1}{4}(\phi_W^{n+1}+\phi_E^n+\phi_N^n+\phi_S^{n+1}-h^2f_P)$

在逐點(diǎn)計(jì)算過(guò)程中， $W$ 和 $S$ 點(diǎn)的值在本次迭代過(guò)程中已知，因此將已知值代入迭代方程中

3.3. SOR迭代

Successive Over Relaxation，逐次超松弛。SOR迭代法收斂的充要條件是松弛因子 $0<\beta <2$ ，當(dāng) $\beta>1$ 時(shí)能夠起到加速收斂的效果。
$\phi_P^{n+1}=\frac{\beta}{4}(\phi_W^{n+1}+\phi_E^n+\phi_N^n+\phi_S^{n+1}-h^2f_P)+(1-\beta)\phi_P^n$
當(dāng) $\beta=1$ ，SOR迭代退化為Gauss-Seidel. 在《Computational Methods for Fluid Dynamics》第四版5.3.3節(jié)給出了一些關(guān)于SOR的討論。

上述方程變換形式可得

$\phi^{n+1}=\phi^n+\beta(\phi_{GS}^{n+1}-\phi^n)$
其中 $\phi_{GS}$ 為Gauss-Seidel法求出的值。進(jìn)一步移項(xiàng)， $\beta$ 用其他幾項(xiàng)表示，可以得出 $\beta>1$ 能夠加速迭代，即縮小初始值變化到最終值所需時(shí)間。

3.4. 迭代收斂標(biāo)準(zhǔn)

本文使用如下標(biāo)準(zhǔn)來(lái)定義收斂標(biāo)準(zhǔn)
$Residual=\frac{1}{N}\sum_{i=1}^{N}\frac{|\phi_i^{n+1}-\phi_i^{n}|}{|\phi_i^{n+1}+\varepsilon|}<10^{-5}$
此處 $N$ 為全局結(jié)點(diǎn)個(gè)數(shù)，上式表示結(jié)點(diǎn)平均相對(duì)偏差小于 $10^{-5}$ 時(shí)，認(rèn)為達(dá)到收斂。 $\varepsilon$ 為極小值，防止分母為0

除此之外，《數(shù)值傳熱學(xué)》還介紹了其他收斂標(biāo)準(zhǔn)。

3.5. 迭代法收斂的分析

《數(shù)值傳熱學(xué)》第二版7.4節(jié)寫(xiě)道，對(duì)于如下形式的方程：
$a_PT_P=\sum a_{nb}T_{nb} +b$

Jacobi與Gauss-Seidel迭代法收斂的一個(gè)充分條件是：系數(shù)矩陣不可約且按行或按列弱對(duì)角占優(yōu)。其中“弱對(duì)角占優(yōu)”需滿足：
$\frac{\sum |a_{nb}|}{|a_P|}\leq 1 \quad or \quad |a_P| \ge \sum |a_{nb}|$
對(duì)各行成立，且其中至少對(duì)一行不等號(hào)成立。在本文的算例中，系數(shù)矩陣的第一行和最后一行對(duì)應(yīng)為不等號(hào)，其他各行均是等號(hào)成立，即 $|a_P| = \sum |a_{nb}|$