计算方法--解线性方程组的迭代法

雅可比迭代法（Jacobi）

对于方程组$\mathbf{A}\mathbf{x}=\mathbf{b}$

将线性方程组的矩阵A分解为
$$A=L+D+U$$
其中
$$\mathbf{L}=\left(\begin{array}{ccccc}0& 0& \cdots & 0& 0\\ a_{21}& 0& \cdots & 0& 0\\ a_{31}& a_{32}& \cdots & 0& 0\\ ⋮& ⋮& & ⋮& ⋮\\ a_{n1}& a_{n2}& \cdots & a_{n,n-1}& 0\end{array}\right)$$
$$\mathbf{D}=\left(\begin{array}{ccccc}{a}_{11}& 0& 0& \cdots & 0\\ 0& a_{22}& 0& \cdots & 0\\ 0& 0& a_{33}& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & a_{nn}\end{array}\right)$$
$$\mathbf{U}=\left(\begin{array}{ccccc}0& a_{12}& a_{13}& \cdots & a_{1n}\\ 0& 0& a_{23}& \cdots & a_{2n}\\ ⋮& ⋮& ⋮& & ⋮\\ 0& 0& 0& \cdots & a_{n-1,n}\\ 0& 0& 0& \cdots & 0\end{array}\right)$$

迭代公式的矩阵形式

$${\mathbf{x}}^{\mathbf{k}+1}={\mathbf{B}}_{\mathbf{J}}{\mathbf{X}}^{\mathbf{k}}+{\mathbf{f}}_{\mathbf{J}}$$
根据如上形式的分解，得：
$$\begin{gathered} {\mathbf{B}}_{\mathbf{J}}=\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A} \\ {\mathbf{f}}_{\mathbf{J}}={\mathbf{D}}^{-1}\mathbf{b} \end{gathered}$$
即
$${\mathbf{x}}^{\mathbf{k}+1}=(\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}){\mathbf{X}}^{\mathbf{k}}+{\mathbf{D}}^{-1}\mathbf{b}$$

编程计算公式：

$$x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum _{j=1,j\ne i}^n{a}_{ij}{x}_j^{(k)}\right),i=1,2,\cdots ,n$$

迭代思路：

将第$i$个方程中的$x_i$提出来，作为$x_i$的迭代式。例如：
$$a_3{x}_3=b-(a_1{x}_1+a_2{x}_2+a_4{x}_4+...)$$

高斯-赛德尔迭代法（Gauss-Seidel）

从直观上看，在收敛的前提下，新的分量$x_1^{k+1},x_2^{k+1},x_3^{k+1},...,x_n^{k+1}$会更好，更精确。根据此思路，有：
$${\mathbf{X}}^{(k+1)}=-(\mathbf{D}+\mathbf{L}{)}^{-1}\mathbf{U}{\mathbf{X}}^{(k)}+(\mathbf{D}+\mathbf{L}{)}^{-1}\mathbf{b}$$
其分量形式为：
$$x_i^{(k+1)}=\frac{1}{a_{ii}}\left(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=i+1}^n{a}_{ij}{x}_j^{(k)}\right),i=1,2,\cdots ,n$$

---

无论是雅可比迭代法还是高斯-赛德尔迭代法，只有在迭代法收敛的条件下才能进行迭代运算。下面介绍判定迭代法收敛的方法。

迭代法的收敛性

迭代法收敛性基本定义

对 $\mathbf{x}=\mathbf{B}\mathbf{x}+\mathbf{f}$ 迭代法收敛的充分必要条件为 $$谱半径\rho (\mathbf{B})<1$$
与初始向量的选取和右端向量 无关。

注：相关定理指出，当雅可比矩阵${\mathbf{B}}_{\mathbf{J}}=\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}$为非负矩阵，雅可比迭代法和高斯-赛德尔迭代法同时收敛或发散。

收敛速度

$$R(B)=-ln(\rho (B))$$
为迭代法的收敛速度。

---

一般，当n较大时，迭代矩阵B的特征值计算较为复杂，迭代法收敛性基本定义条件较难验证，因此可利用定理结论$$\rho (b)\le ||B||$$作为$\rho (B)$上界的一种估计。
于是得如下定理

迭代法充分条件1

若$$||B||<1$$则迭代公式$$x=Bx+f$$收敛
（注：是充分条件，不是必要条件）

在实际计算中可将$||X^{(k)}-X^{(k-1)}|{|}_{\infty }<\epsilon$作为终止条件。

迭代法充分条件2

若$Ax=b$的系数矩阵$A$为 严格对角占优矩阵 或 不可约对角占优矩阵，则雅可比迭代法和高斯-赛德尔迭代法收敛。

迭代法其他收敛条件

• 若$Ax=b$的系数矩阵$A$对称正定，则高斯-赛德尔迭代法收敛。
• 若$A$对称正定,$2D-A$也对称正定（D为A的对角元组成的对角阵），则雅可比迭代法收敛。
• 若$A$对称正定,$2D-A$非正定，则雅可比迭代法 不收敛。

---

逐次超松弛迭代是一种加速迭代的技术，下面以Jacobi迭代法来阐述这种收敛技术。

JOR迭代法

将
$${\mathbf{x}}^{\mathbf{k}+1}=(\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}){\mathbf{X}}^{\mathbf{k}}+{\mathbf{D}}^{-1}\mathbf{b}$$
加入松弛因子 $\omega$ ，得：

$${\mathbf{x}}^{\mathbf{k}+1}=(\mathbf{I}-\omega {\mathbf{D}}^{-1}\mathbf{A}){\mathbf{X}}^{\mathbf{k}}+\omega {\mathbf{D}}^{-1}\mathbf{b}$$
其分量形式为：
$$x_i^{(k+1)}=x_i^{(k)}+\frac{\omega }{a_{ii}}\left(b_i-\sum _{j=1}^n{a}_{ij}{x}_j^{(k)}\right),i=1,2,\cdots ,n$$

当 Jacobi 迭代方法

1. 迭代矩阵 ${\mathbf{B}}_J$ 的特征值为实数
2. ${\mathbf{B}}_J$ 的谱半径 $\rho \left({\mathbf{B}}_J\right)<1$

JOR 迭代方法有最佳松弛因子为
$${\omega }_{\mathrm{o}\mathrm{p}\mathrm{t}}=\frac{2}{2-{\lambda }_{\max }^{B_J}-{\lambda }_{\min }^{B_J}}$$

JOB收敛准则

若Jacobi迭代法收敛，则松弛因子 $\omega \in [0,1)$ 的JOB迭代法也收敛。

SOR迭代法

以下, 我们将 Gauss-Seidel 迭代法改造为一类 逐次超松弛 迭代格式.
将Gauss-Seidel 迭代法
$${\mathbf{X}}^{(k+1)}=-(\mathbf{D}+\mathbf{L}{)}^{-1}\mathbf{U}{\mathbf{X}}^{(k)}+(\mathbf{D}+\mathbf{L}{)}^{-1}\mathbf{b}$$
转化为如下形式
$$\mathbf{D}{\mathbf{X}}^{(k+1)}=\mathbf{D}{\mathbf{X}}^{(k)}+\left[\mathbf{b}-\mathbf{L}{\mathbf{X}}^{(k+1)}-(\mathbf{D}+\mathbf{U}){\mathbf{X}}^{(k)}\right]$$

今在上式右端引人松弛因子 $\omega >0$ , 化简后得 :
$${\mathbf{X}}^{(k+1)}=(\mathbf{D}+\omega \mathbf{L}{)}^{-1}\left\{[(1-\omega )\mathbf{D}-\omega \mathbf{U}]{\mathbf{X}}^{(k)}+\omega \mathbf{b}\right\}$$

其分量形式为：
$$x_i^{(k+1)}=x_i^{(k)}+\frac{\omega }{a_{ii}}\left(b_i-\sum _{j=1}^{i-1}{a}_{ij}{x}_j^{(k+1)}-\sum _{j=i}^n{a}_{ij}{x}_j^{(k)}\right)，i=1,2,\cdots ,n$$
SOR 迭代方法有最佳松弛因子为
$${\omega }_{\mathrm{o}\mathrm{p}\mathrm{t}}=\frac{2}{1+\sqrt{1-{\rho }^2\left[-D^{-1}(L+U)\right]}}$$