数值计算：线性方程组的迭代解法 01 静态迭代法

对于线性方程组的迭代求解方法可以分为两类，静态迭代方法与非静态迭代方法，两者区别在于，前者构造简单，迭代步长与方向恒定，但是收敛条件限制较大，收敛速度较慢。而非静态方法构造格式更复杂，收敛速度更快。本文主要记录静态迭代方法

静态迭代法

考虑以下线性方程组

\[\boldsymbol Ax=\boldsymbol b \]

对于工程计算中产生的大型稀疏矩阵A（非奇异），迭代法是求解此类方程组的最佳方法。根据此方程构造其一阶定常迭代格式迭

\[\begin{cases} x^{(0)}\qquad \text {初始值}\\ \boldsymbol x^{k+1}=\boldsymbol Bx^{k}+d \end{cases} \]

迭代方法的收敛条件可以简单认作：\(\boldsymbol B\)的谱半径小于\(1\)，即\(\rho(B)<1\)

收敛速度可以简单认作：\(R(\boldsymbol B)=-\text{ln}\rho(\boldsymbol B)\)

---

对于方程\(\boldsymbol Ax=\boldsymbol b\)，一般选择将\(\boldsymbol A\)进行分裂，\(\boldsymbol {A=M-N}\)，可以选择合适的\(\boldsymbol M\)使得\(\boldsymbol {M}x=d\)易于求解，

\[\boldsymbol Ax=b \Leftrightarrow x=\boldsymbol M^{-1}\boldsymbol Nx+\boldsymbol M^{-1}b \]

Jacobi迭代法

将系数矩阵\(\boldsymbol A\)分解

\[\boldsymbol {A=D-L-U} \]

\(\boldsymbol {D,L,U}\)分别为对角矩阵、下三角矩阵、上三角矩阵，则可以构造得到\(Jacobi\)迭代格式

\[\begin{cases} x^{(0)}\qquad \text {初始值}\\ \boldsymbol x^{(k+1)}=\boldsymbol {D^{-1}(L+U)}x^{(k)}+\boldsymbol{D^{-1}b} \end{cases} \]

function [x,iter,Residual]=Jacobi(A,b,x0,epsilon,iter_max)
%系数矩阵A,b,初始值x0,误差限制epsilon，最大迭代步数iter_max
iter=0;
x=x0;
%D=diag(diag(A));    %对角线
invD=diag(diag(A).^(-1));
U=triu(A,1);        %上三角
L=tril(A,-1);       %下三角
B=-invD*(L+U);
d=invD*b;
if(max(abs(eig(B)))>=1)
    error('can not convergent!')
end
Residual=1e20;
while sqrt(Residual) >= epsilon && iter < iter_max
    iter=iter+1;
    x_new=B*x+d;
    Residual=norm(x_new-x,2);
    x=x_new;
end
end

Gauss-Seidel迭代法

与\(Jacobi\)迭代法不同的是，迭代过程中，更新计算第\(i+1\)个分量时，使用了已经更新的第\(i\)个变量。迭代形式为

\[\begin{cases} x^{(0)}\qquad \text {初始值}\\ \boldsymbol x^{(k+1)}=\boldsymbol {{(D-L)}^{-1}U}x^{(k)}+\boldsymbol{(D-L)^{-1}b} \end{cases} \]

SOR迭代

逐次超松弛迭代法（Successive Over Relaxation method）

选取分裂的下三角矩阵\(M\)中带有松弛因子

\[\boldsymbol M =\frac{1}{\omega}(\boldsymbol D-\omega \boldsymbol L),\quad \omega>0 \]

带有松弛因子的迭代格式为

\[\begin{cases} x^{(0)}\qquad \text {初始值}\\ \boldsymbol x^{(k+1)}=\boldsymbol {{(D-\omega L)}^{-1}((1-\omega)D+\omega U)}x^{(k)}+ \boldsymbol{\omega(D-\omega L)^{-1}b} \end{cases} \]

显然\(\omega=1\)时，\(SOR\)方法为\(Gauss-Seidel\)迭代法

\(\omega>1\)，称为超松弛迭代；\(\omega<1\)，称为亚松弛迭代；

---

function [x,iter,Residual]=SOR(A,b,x0,omega,epsilon,iter_max)
%系数矩阵A,b,初始值x0,误差限制epsilon,松弛因子omega,最大迭代步数iter_max
iter=0;
x=x0;
D=diag(diag(A));    %对角线
U=triu(A,1);        %上三角
L=tril(A,-1);       %下三角
B=inv(D-omega*L)*((1-omega)*D+omega*U);
d=omega*inv(D-omega*L)*b;
Residual=1e20;
while sqrt(Residual) >= epsilon && iter < iter_max
    iter=iter+1;
    x_new=B*x+d;
    Residual=norm(x_new-x,2);
    x=x_new;
end
end