对于线性方程组的迭代求解方法可以分为两类,静态迭代方法与非静态迭代方法,两者区别在于,前者构造简单,迭代步长与方向恒定,但是收敛条件限制较大,收敛速度较慢。而非静态方法构造格式更复杂,收敛速度更快。本文主要记录静态迭代方法
考虑以下线性方程组
\[\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 \]将系数矩阵\(\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
与\(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} \]逐次超松弛迭代法(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