5 线性方程组的迭代解法
5.1 引言
直接法是最简单也最容易想到的方法，除直接法之外，还有基于迭代的方式求解的方法，基本思想是将原始的方程组\(Ax=b\)，转换为：
\[x^{(k+1)}=Bx^{(k)}+g
\]
的形式进行迭代，如果向量序列\(\{x^{(k)}\}\)能够收敛到某个向量，那么最终的收敛向量就是所要求的最终解。
由于实际应用中迭代次数总是有限的，因此通过迭代法得到的解只是理论值的一个近似值。为了保证近似值足够接近理论值，我们通过设定一个误差上限，使得最终结果能够满足精度要求。
5.2 雅可比（Jacobi）迭代
假设有方程组\(Ax=b\)：
\[\left\{\begin{aligned}
a_{11}x_1+a_{12}x_2&+\cdots+a_{1n}x_n=b_1\\
a_{21}x_1+a_{22}x_2&+\cdots+a_{2n}x_n=b_2\\
&\cdots\\
a_{n1}x_1+a_{n2}x_2&+\cdots+a_{nn}x_n=b_n
\end{aligned}\right.
\]
并假设\(A\)是非奇异矩阵，\(a_{ii}\ne0(1\le i\le n)\)，那么
\[\left\{\begin{aligned}
x_1&=-\frac{1}{a_{11}}(a_{12}x_2+\cdots+a_{1n}x_n-b_1)\\
x_2&=-\frac{1}{a_{22}}(a_{21}x_1+a_{23}x_3+\cdots+a_{2n}x_n-b_2)\\
&\qquad\qquad\vdots\\
x_n&=-\frac{1}{a_{nn}}(a_{n1}x_1+\cdots+a_{n(n-1)}x_{n-1}-b_n)\\
\end{aligned}\right.
\]
将上述式子改写为迭代格式，就得到了Jacobi迭代格式：
\[\left\{\begin{aligned}
x_1^{(k+1)}&=-\frac{1}{a_{11}}(a_{12}x_2^{(k)}+\cdots+a_{1n}x_n^{(k)}-b_1)\\
x_2^{(k+1)}&=-\frac{1}{a_{22}}(a_{21}x_1^{(k)}+a_{23}x_3^{(k)}+\cdots+a_{2n}x_n^{(k)}-b_2)\\
&\qquad\qquad\vdots\\
x_n^{(k+1)}&=-\frac{1}{a_{nn}}(a_{n1}x_1^{(k)}+\cdots+a_{n(n-1)}x_{n-1}^{(k)}-b_n)\\
\end{aligned}\right.\tag{1}
\]
也可以写成
\[x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum_{j=1\\j\ne i}^na_{ij}x_j^{(k)}),
i=1,2,\cdots,n
\]
为了利用计算机并行计算的能力，将上述迭代格式表示成矩阵的形式，令\(A=D-L-U\)，其中
\[D=\left[\begin{matrix}
a_{11}&&0\\
&\ddots&\\
0&&a_{nn}
\end{matrix}\right],
L=\left[\begin{matrix}
0&&&0\\
-a_{21}&0&&\\
\vdots&\ddots&\ddots&\\
-a_{n1}&\cdots&-a_{n(n-1)}&0
\end{matrix}\right],
U=\left[\begin{matrix}
0&-a_{12}&\cdots&-a_{1n}\\
&0&\ddots&\vdots\\
&&\ddots&-a_{(n-1)n}\\
0&&&0
\end{matrix}\right]
\]
即\(D\)包含主对角线上的元素，其余位置全为0。\(L\)包含主对角线下方区域所有元素（取相反数），其余位置全为0。\(U\)包含主对角线上方区域所有元素（取相反数），其余位置全为0。
因此，方程组\(Ax=b\)可以改写为
\[\begin{aligned}
Ax=b&\Leftrightarrow (D-L-U)x=b\\
&\Leftrightarrow Dx=(L+U)x+b\\
&\Leftrightarrow x=D^{-1}(L+U)x+D^{-1}b
\end{aligned}
\]
令\(B_J=D^{-1}(L+U),g_J=D^{-1}b\)，那么Jacobi迭代的矩阵形式为
\[x^{(k+1)}=B_Jx^{(k)}+g_J
\]
Jacobi矩阵迭代格式也是可以直接根据迭代格式\((1)\)写出来的，请读者好好思考其中的对应关系。
5.3 高斯-赛德尔（Gauss-Seidel）迭代
为了加快收敛速度，同时节省计算机内存，可以对Jacobi迭代进行如下改进：每计算出一个分量的近似值，就将其应用到下一个分量的计算过程中，即迭代格式变为
\[\left\{\begin{aligned}
x_1^{(k+1)}&=-\frac{1}{a_{11}}(a_{12}x_2^{(k)}+\cdots+a_{1n}x_n^{(k)}-b_1)\\
x_2^{(k+1)}&=-\frac{1}{a_{22}}(a_{21}x_1^{(k+1)}+a_{23}x_3^{(k)}+\cdots+a_{2n}x_n^{(k)}-b_2)\\
&\qquad\qquad\vdots\\
x_n^{(k+1)}&=-\frac{1}{a_{nn}}(a_{n1}x_1^{(k+1)}+\cdots+a_{n(n-1)}x_{n-1}^{(k+1)}-b_n)\\
\end{aligned}\right.
\]
也可以将上述格式归结为
\[x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i+1}^{n}a_{ij}x_j^{(k)}),i=1,2,\cdots,n
\]
我们把Jacobi迭代的矩阵格式写成
\[Dx^{(k+1)}=Lx^{(k)}+Ux^{(k)}+b
\]
那么高斯-赛德尔迭代就是将其变为
\[Dx^{(k+1)}=Lx^{(k+1)}+Ux^{(k)}+b
\]
因此高斯-赛德尔迭代的矩阵格式为
\[x^{(k+1)}=B_{S}x^{(k)}+g_S
\]
其中\(B_S=(D-L)^{-1}U,g_S=(D-L)^{-1}b\)。
5.4 SOR（Successive Over Relaxation）迭代
SOR迭代是对高斯-赛德尔迭代的进一步改进，其基本思想是先利用当前值通过高斯-赛德尔迭代计算一个估计值，然后再将这个估计值和当前值做加权平均，以获得一个更好的结果。其迭代格式为
\[x_i^{(k+1)}=\omega{\hat{x}_i^{(k+1)}}+(1-\omega)x_i^{(k)}
\]
其中\(\hat{x}_i^{(k+1)}\)是通过高斯-赛德尔迭代计算出来的结果，因此
\[x_i^{(k+1)}=(1-\omega)x_i^{(k)}+\frac{\omega}{a_{ii}}(b_i-\sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum_{j=i+1}^{n}a_{ij}x_j^{(k)}),i=1,2,\cdots,n
\]
其中\(\omega\)称为松弛因子。当\(\omega>1\)时，称为超松弛；\(\omega<1\)时，称为低松弛；当\(\omega=1\)时，就是高斯-赛德尔迭代。
根据上述迭代格式，我们可以得到
\[x^{(k+1)}=(1-\omega)x^{(k)}+\omega{D^{-1}b}+\omega{D^{-1}Lx^{(k+1)}}+\omega{D^{-1}Ux^{(k)}}
\]
整理之后可以得到
\[x^{(k+1)}=B_\omega x^{(k)}+g_\omega
\]
其中\(B_\omega=(D-\omega{L})^{-1}[(1-\omega)D+\omega{U}],g_\omega=\omega{(D-\omega L)^{-1}}b\)。
5.5 收敛性判别
上面三种迭代格式都可以写成\(x^{(k+1)}=Bx^{(k)}+g\)的形式，因此，在这一部分将探讨在什么情况下，上述三种迭代格式能够收敛。并给出误差估计式。
5.5.1 收敛的充要条件
在给出收敛的充要条件前，这里先引入几个定义和定理。
定义1 设\(\lambda_i(1\le i\le n)\)是矩阵\(A\)特征值，则称\(\rho(A)=\max\{|\lambda_i|\}\)为矩阵\(A\)的谱半径。
定义2 设有矩阵序列\(\{B_k=[b_{ij}^{(k)}]_{n\times n}\}\)，以及矩阵\(B=[b_{ij}]_{n\times n}\)。如果
\[\lim_{k\rightarrow\infty}b_{ij}^{(k)}=b_{ij}\quad(i,j=1,2,\cdots,n)
\]
则称矩阵序列\(\{B_k\}\)收敛于矩阵\(B\)，记为\(\lim_{k\rightarrow\infty}{B_k}=B\)。
定理1 \(\lim_{k\rightarrow\infty}{B_k}=B\)的充要条件是\(\lim_{k\rightarrow\infty}||{B_k-B}||=0\)。
定理2 \(\lim_{k\rightarrow\infty}{B_k}=O\)的充要条件是\(\rho(B)<1\)。
上述定理的证明涉及到Jordan标准型理论，具体证明可以参考矩阵分析等书籍。
定理3 对于任意\(A\in\mathbb{R}^{n\times n}\)，都有\(\rho(A)\le||A||\)。
证明：设\(\lambda_i\)和\(x_i\)分别是矩阵\(A\)对应的一个特征值和特征向量，则有\(Ax_i=\lambda_i{x_i}\)，等式两边取范数，可以得到
\[||\lambda_i||\cdot||x_i||=||\lambda_ix_i||=||Ax_i||\le||A||\cdot||x_i||
\]
因为\(||x_i||\ne0\)，所以\(||A||\ge||\lambda_i||\)。对于任意的特征值该结论都是成立的，因此
\[\rho(A)=\max_{1\le i\le n}{|\lambda_i|}\le||A||
\]
定理4 迭代格式\(x^{(k+1)}=Bx^{(k)}+g\)对任意的初始向量\(x^{(0)}\in\mathbb{R}^{n}\)产生的向量序列\(\{x^{(k)}\}\)都收敛的充要条件是\(\rho(B)<1\)。
必要性证明：
设\(\lim{x^{(k)}}=x^*\)，则有\(x^*=Bx^*+g\)，那么
\[x^{(k+1)}-x^*=(Bx^{(k)}+g)-(Bx^*+g)=B(x^{(k)}-x^*)=\cdots=B^{k+1}(x^{(0)}-x^*)\tag{4-1}
\]
因此
\[\lim_{k\rightarrow\infty}(x^{(k+1)}-x^*)=
(x^{(0)}-x^*)\lim_{k\rightarrow\infty}B^{k+1}=\mathbf0
\]
因为\(x^{(0)}\)可以是任意的初始向量，为了使得上述极限式恒成立，就必须要求\(\lim_{k\rightarrow\infty}{B_k}=O\)。
因此，根据定理2可知，\(\rho(B)<1\)。
充分性证明：
如果\(\rho(B)<1\)，那么1不是\(B\)的特征值，进而\(I-B\)是可逆矩阵（因为不存在为0的特征值），于是\((I-B)x=g\)会有唯一解\(x^*\)，即\(x^*=Bx^*+g\)成立。基于这一点，我们可以推出\((4-1)\)式是成立的。
根据定理2可知，\(\lim_{k\rightarrow\infty}{B^{k+1}}=O\)成立。于是对于任意的初始向量\(x^{(0)}\)，根据\((4-1)\)式都有
\[\lim_{k\rightarrow\infty}(x^{(k+1)}-x^*)=
(x^{(0)}-x^*)\lim_{k\rightarrow\infty}B^{k+1}=\mathbf0
\]
因此\(\lim{x^{(k)}}=x^*\)成立。
综上可知，定理4得证。
因此，判断雅可比迭代、高斯-赛德尔迭代、SOR迭代是否收敛，等价于判断\(\rho(B_J),\rho(B_S),\rho(B_\omega)\)的值是否小于1。
5.5.2 收敛的充分条件
定理4虽然给出了收敛的充要条件，但是谱半径不容易计算。因此，在这一部分给出几个判断收敛性的充分条件。
定理5 若矩阵\(B\)的某个范数满足\(||B||<1\)，则迭代格式\(x^{(k+1)}=Bx^{(k)}+g\)对任意的初始向量\(x^{(0)}\in\mathbb{R}^{n}\)产生的向量序列\(\{x^{(k)}\}\)都收敛。并且有
\[||x^{(k)}-x^*||\le\frac{||B||}{1-||B||}||x^{(k)}-x^{(k-1)}||\\
||x^{(k)}-x^*||\le\frac{||B||^k}{1-||B||}||x^{(1)}-x^{(0)}||
\]
定理5指出了控制迭代误差和次数的方法。
证明：
根据定理3和定理4直接就可以知道当\(||B||<1\)时，迭代格式收敛。
\[\begin{aligned}
||x^{(k)}-x^*||&=||(Bx^{(k-1)}+g)-(Bx^*+g)||\\
&=||B(x^{(k-1)}-x^*)||\\
&=||B(x^{(k-1)}-x^{(k)}+x^{(k)}-x^*)||\\
&\le\lVert{B}\rVert\lVert{x^{(k)}-x^{(k-1)}}\rVert
+\lVert{B}\rVert\lVert{x^{(k)}-x^*}\rVert
\end{aligned}
\]
所以
\[(1-||B||)\lVert{x^{(k)}-x^*}\rVert\le\lVert{B}\rVert\lVert{x^{(k)}-x^{(k-1)}}\rVert
\]
因此有
\[||x^{(k)}-x^*||\le\frac{||B||}{1-||B||}||x^{(k)}-x^{(k-1)}||
\]
又因为
\[\begin{aligned}
||x^{(k)}-x^{(k-1)}||&=||(Bx^{(k-1)}+g)-(Bx^{(k-2)}+g)||\\
&=||B(x^{(k-1)}-x^{(k-2)})||\\
&\le\lVert{B}\rVert\lVert{x^{(k-1)}-x^{(k-2)}}\rVert\\
&\le\cdots\\
&\le\lVert{B}\rVert^{k-1}\lVert{x^{(1)}-x^{(0)}}\rVert
\end{aligned}
\]
于是有
\[||x^{(k)}-x^*||\le\frac{||B||}{1-||B||}||x^{(k)}-x^{(k-1)}||
\le\frac{||B||^k}{1-||B||}||x^{(1)}-x^{(0)}||
\]
定理6 SOR迭代收敛的必要条件是\(0\lt\omega\lt2\)。
证明：如果SOR迭代收敛，那么\(\rho(B_\omega)<1\)，并且有如下结论
\[|\det{B_\omega}|=|\lambda_1\lambda_2\cdots\lambda_n|\le\rho(B_\omega)^n<1
\]
其中\(\lambda_1,\cdots,\lambda_n\)是\(B_\omega\)的\(n\)个特征值。
\[\begin{aligned}
\det{B_\omega}&=\det\{{(D-\omega{L})^{-1}}[(1-\omega)D+\omega U]\}\\
&=\frac{\det[{(1-\omega)D+\omega U}]}{\det{(D-\omega L)}}
\end{aligned}
\]
考虑到\(D-\omega L\)和\((1-\omega)D+\omega U\)实际上是对角矩阵，所以
\[|\det{B_\omega}|=\frac{|(1-\omega)a_{11}a_{22}\cdots a_{nn}|}
{|a_{11}a_{22}\cdots a_{nn}|}=|1-\omega|<1
\]
于是\(0\lt\omega\lt2\)。
下面在介绍两个判断收敛性的充分非必要条件。
定义3 如果矩阵\(A=[a_{ij}]_{n\times n}\)满足以下条件
\[|a_{ii}|>\sum_{j=1\\j\ne i}^n{|a_{ij}|}
\]
那么就称\(A\)是严格对角占优的。
定理7 如果矩阵\(A\)是严格对角占优的，那么\(Ax=b\)方程的雅可比迭代和高斯-赛德尔迭代对于任意初始向量都是收敛的。
定理8 如果矩阵\(A\)是对称正定的，那么\(Ax=b\)方程的高斯-赛德尔迭代对于任意初始向量都是收敛的。
上述几个充分条件相对来说都还是比较好用的，并且对于任意初始向量收敛，这个性质是非常好的。这意味着即使在某一个计算过程中\(x^{(k)}\)发生了溢出，这也不是太大的问题。因为，这相当于我们之前所计算的结果全部作废，重新从一个新的初始向量开始迭代而已。
