1. 问题描述

本算例来自B站Up主“Red-Green鲤鱼”的系列教程。本文主要介绍计算代数方程组的三种点迭代方法。

1.1. 泊松方程

含有二阶偏导数的偏微分方程：
$\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}=f(x,y)$

当 $f=0$ 时，上述方程被称为拉普拉斯方程。许多物理过程都可以用泊松方程来描述，如热传导方程
$\frac{\partial \phi}{\partial t}=\frac{\partial^2 \phi}{\partial x^2}+\frac{\partial^2 \phi}{\partial y^2}$

在求解不可压缩流动的NS方程时，通常将已知压力场代入动量方程来预估速度场，然后将预估的速度场代入连续性方程中，由于预估的速度场中包含压力梯度项，因此代入连续性方程后会得到 $\nabla \cdot(\nabla P)$ ，即压力泊松方程

1.2. 算例

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

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

比较数值解与解析解在定义域 $x\in [-1,1], y\in [-1,1]$ 上的差别。边界条件为Dirichlet，边界上的值由解析解求出。

2. 区域离散和方程离散

$x$ 方向设置 $N$ 个结点，编号从 $1-N$ ， $y$ 方向上设置 $M$ 个结点，编号从 $1-M$ ，结点间距分别为 $\Delta x$ 和 $\Delta y$ ，将 $(i,j)$ 号结点记为 $P$ ，该结点上下左右四个结点分别记为 $N,S,W,E$

采用有限差分法计算泊松方程，二阶偏导数项采用二阶中心差分离散，
$\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. 代数方程组求解

上一篇文章介绍了代数方程组求解方法中的直接解法TDMA，本文介绍另一大类求解方法：迭代法。迭代法的思想也可以概况为“预测-校正”，给出初始值，通过不断迭代逐步改进，直到达到一定精度要求为止。
该方法需要首先构造迭代方式；其次是所构造的迭代序列是否收敛，如果收敛则要进一步提高收敛速度。

迭代法可分为点迭代、块迭代、交替方向迭代法以及强隐迭代法。在点迭代法中，每一步计算只能改进求解区域中一个结点的值，且该值是由一个显函数形式由其余各点的已知值来确定，因而点迭代法又称为显式迭代法。

下面将讨论点迭代法的三种实施方式。

3.1. 雅可比迭代

$\phi_P^{n+1}=\frac{1}{4}(\phi_W^n+\phi_E^n+\phi_N^n+\phi_S^n-h^2f_P)$
上式中上标 $n$ 为当前预测值， $n+1$ 为代入迭代方程后的校正值

3.2. 高斯-赛德尔迭代

$\phi_P^{n+1}=\frac{1}{4}(\phi_W^{n+1}+\phi_E^n+\phi_N^n+\phi_S^{n+1}-h^2f_P)$

在逐点计算过程中， $W$ 和 $S$ 点的值在本次迭代过程中已知，因此将已知值代入迭代方程中

3.3. SOR迭代

Successive Over Relaxation，逐次超松弛。SOR迭代法收敛的充要条件是松弛因子 $0<\beta <2$ ，当 $\beta>1$ 时能够起到加速收敛的效果。
$\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$
当 $\beta=1$ ，SOR迭代退化为Gauss-Seidel. 在《Computational Methods for Fluid Dynamics》第四版5.3.3节给出了一些关于SOR的讨论。

上述方程变换形式可得

$\phi^{n+1}=\phi^n+\beta(\phi_{GS}^{n+1}-\phi^n)$
其中 $\phi_{GS}$ 为Gauss-Seidel法求出的值。进一步移项， $\beta$ 用其他几项表示，可以得出 $\beta>1$ 能够加速迭代，即缩小初始值变化到最终值所需时间。

3.4. 迭代收敛标准

本文使用如下标准来定义收敛标准
$Residual=\frac{1}{N}\sum_{i=1}^{N}\frac{|\phi_i^{n+1}-\phi_i^{n}|}{|\phi_i^{n+1}+\varepsilon|}<10^{-5}$
此处 $N$ 为全局结点个数，上式表示结点平均相对偏差小于 $10^{-5}$ 时，认为达到收敛。 $\varepsilon$ 为极小值，防止分母为0

除此之外，《数值传热学》还介绍了其他收敛标准。

3.5. 迭代法收敛的分析

《数值传热学》第二版7.4节写道，对于如下形式的方程：
$a_PT_P=\sum a_{nb}T_{nb} +b$

Jacobi与Gauss-Seidel迭代法收敛的一个充分条件是：系数矩阵不可约且按行或按列弱对角占优。其中“弱对角占优”需满足：
$\frac{\sum |a_{nb}|}{|a_P|}\leq 1 \quad or \quad |a_P| \ge \sum |a_{nb}|$
对各行成立，且其中至少对一行不等号成立。在本文的算例中，系数矩阵的第一行和最后一行对应为不等号，其他各行均是等号成立，即 $|a_P| = \sum |a_{nb}|$