二维泊松方程求解-SIP-最速下降法-共轭梯度

1. 直接解法：LU分解

在前面的内容中曾经提到，使用有限差分或有限体积法通过隐式离散得到$A\varphi =Q$的求解形式，其中$A$为系数矩阵。在一定条件下，$A$能够通过因式分解为$A=LU$，其中$L$为下三角矩阵，$U$为上三角矩阵。

这样的分解方式在高斯消元中十分有用，对$A\varphi =Q$的求解可分为以下两步
$$\begin{gathered} U\varphi =Y \\ LY=Q \end{gathered}$$

2. 迭代法：incomplete LU decomposition

如果存在一个与$A$近似的矩阵$M$，对$M$做LU分解，我们把这样的步骤称为$A$的不完全LU分解，ILU，即
$$M=LU=A+N$$
其中$N$为小量。

2.1. SIP

Stone提出了基于ILU且具有良好普适性的strongly implicit procedure (SIP)，接下来我们简要介绍一下该算法。

在ILU中，如果矩阵$A$中的$(i,j)$元素非零，那么$L$或$U$矩阵在该元素对应位置上的对角线非零。以下图所示的$A$为五对角矩阵为例说明，图中非零元素占据五条对角线，那么对应的$L$和$U$各自包含三条对角线。

令$U$主对角线上元素为1。可见，$L$与$U$的乘积不可能为五对角矩阵，而是七对角矩阵，因此不可能还原为$A$，而是$LU=M$，如下图所示

令两次相邻迭代的差值为$\delta$，则原代数方程组可写成
$$A{\varphi }^{n+1}-A{\varphi }^n=Q-A{\varphi }^n$$

$$A{\delta }^n=R^n$$
矩阵$R$为残差矩阵，为已知值。设想如果$A$能够完全LU分解，那么便可以通过直接解法解出${\delta }^n$，然后加到${\varphi }^n$求出${\varphi }^{n+1}$，整个过程只需一次迭代。可是根据前文的分析，五对角矩阵$A$无法完全LU分解，因此Stone的目标变成了尽可能的让$LU$的乘积$M$接近$A$。计算$M\delta$的乘积可得
$$\begin{gathered} M\delta =LU\delta =M_P{\delta }_P+M_S{\delta }_S+M_N{\delta }_N+M_E{\delta }_E+M_W{\delta }_W+ \\ {M}_{NW}{\delta }_{NW}+M_{SE}{\delta }_{SE} \end{gathered}$$

而$A\delta$的乘积只有五项，对比来看多出了最后两项。因此Stone认为应尽量地消去最后两项。对${\delta }_{NW}$做Taylor展开
$${\delta }_{NW}\approx {\delta }_P-\frac{\partial \delta }{\partial x}{|}_P\Delta x+\frac{\partial \delta }{\partial y}{|}_P\Delta y$$
$$\frac{\partial \delta }{\partial x}{|}_P\approx \frac{{\delta }_P-{\delta }_W}{\Delta x}$$
∂ δ ∂ y ∣ P ≈ δ N − δ P Δ y \frac{\partial \delta}{\partial y}\bigg|_P\approx \frac{\delta_N-\delta_P}{\Delta y} ∂y∂δ​