數值計算：線性方程組的迭代解法 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