MATLAB代码--解线性方程组的迭代法

解线性方程组的迭代法1.rs里查森迭代法求线性方程组Ax=b的解function[x,n]=rs(A,b,x0,eps,M)if(nargin==3)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==4)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%迭代过程while(toleps)x=(I-A)*x0+b;n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0);x0=x;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend2.crs里查森参数迭代法求线性方程组Ax=b的解function[x,n]=crs(A,b,x0,w,eps,M)if(nargin==4)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==5)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%迭代过程while(toleps)x=(I-w*A)*x0+w*b;n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0);x0=x;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend3.grs里查森迭代法求线性方程组Ax=b的解function[x,n]=grs(A,b,x0,W,eps,M)if(nargin==4)eps=1.0e-6;%eps表示迭代精度M=10000;%M表示迭代步数的限制值elseif(nargin==5)M=10000;endI=eye(size(A));n=0;x=x0;tol=1;%前后两次迭代结果误差%迭代过程while(toleps)x=(I-W*A)*x0+W*b;%迭代公式n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0);x0=x;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend4.jacobi雅可比迭代法求线性方程组Ax=b的解function[x,n]=jacobi(A,b,x0,eps,varargin)ifnargin==3eps=1.0e-6;M=200;elseifnargin3errorreturnelseifnargin==5M=varargin{1};endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;n=1;%迭代次数whilenorm(x-x0)=epsx0=x;x=B*x0+f;n=n+1;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend5.gauseidel高斯-赛德尔迭代法求线性方程组Ax=b的解function[x,n]=gauseidel(A,b,x0,eps,M)ifnargin==3eps=1.0e-6;M=200;elseifnargin==4M=200;elseifnargin3errorreturn;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1;%迭代次数whilenorm(x-x0)=epsx0=x;x=G*x0+f;n=n+1;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend6.SOR超松弛迭代法求线性方程组Ax=b的解function[x,n]=SOR(A,b,x0,w,eps,M)ifnargin==4eps=1.0e-6;M=200;elseifnargin4errorreturnelseifnargin==5M=200;endif(w=0||w=2)error;return;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=B*x0+f;n=1;%迭代次数whilenorm(x-x0)=epsx0=x;x=B*x0+f;n=n+1;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend7.SSOR对称逐次超松弛迭代法求线性方程组Ax=b的解function[x,n]=SSOR(A,b,x0,w,eps,M)ifnargin==4eps=1.0e-6;M=200;elseifnargin4errorreturnelseifnargin==5M=200;endif(w=0||w=2)error;return;endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B1=inv(D-L*w)*((1-w)*D+w*U);B2=inv(D-U*w)*((1-w)*D+w*L);f1=w*inv((D-L*w))*b;f2=w*inv((D-U*w))*b;x12=B1*x0+f1;x=B2*x12+f2;n=1;%迭代次数whilenorm(x-x0)=epsx0=x;x12=B1*x0+f1;x=B2*x12+f2;n=n+1;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend8.JOR雅可比超松弛迭代法求线性方程组Ax=b的解function[x,n]=JOR(A,b,x0,w,eps,M)ifnargin==4eps=1.0e-6;M=10000;elseifnargin==5M=10000;endif(w=0||w=2)%收敛条件要求error;return;endD=diag(diag(A));%求A的对角矩阵B=w*inv(D);%迭代过程x=x0;n=0;%迭代次数tol=1;%迭代过程whiletol=epsx=x0-B*(A*x0-b);n=n+1;tol=norm(x-x0);x0=x;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend9.twostep两步迭代法求线性方程组Ax=b的解function[x,n]=twostep(A,b,x0,eps,varargin)ifnargin==3eps=1.0e-6;M=200;elseifnargin3errorreturnelseifnargin==5M=varargin{1};endD=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角阵U=-triu(A,1);%求A的上三角阵B1=(D-L)\U;B2=(D-U)\L;f1=(D-L)\b;f2=(D-U)\b;x12=B1*x0+f1;x=B2*x12+f2;n=1;%迭代次数whilenorm(x-x0)=epsx0=x;x12=B1*x0+f1;x=B2*x12+f2;n=n+1;if(n=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend10.fastdown最速下降法求线性方程组Ax=b的解function[x,n]=fastdown(A,b,x0,eps)if(nargin==3)eps=1.0e-6;endx=x0;n=0;tol=1;while(toleps)%以下过程可参考算法流程r=b-A*x0;d=dot(r,r)/dot(A*r,r);x=x0+d*r;tol=norm(x-x0);x0=x;n=n+1;end11.conjgrad共轭梯度法求线性方程组Ax=b的解function[x,n]=conjgrad(A,b,x0)r1=b-A*x0;p=r1;n=0;fori=1:rank(A)%以下过程可参考算法流程if(dot(p,A*p)1.0e-50)%循环结束条件break;endalpha=dot(r1,r1)/dot(p,A*p);x=x0+alpha*p;r2=r1-alpha*A*p;if(r21.0e-50)%循环结束条件break;endbelta=dot(r2,r2)/dot(r1,r1);p=r2+belta*p;n=n+1;end12.preconjgrad预处理共轭梯度法求线性方程组Ax=b的解function[x,n]=preconjgrad(A,b,x0,M,eps)ifnargin==4eps=1.0e-6;endr1=b-A*x0;iM=inv(M);z1=iM*r1;p=z1;n=0;tol=1;whiletol=epsalpha=dot(r1,z1)/dot(p,A*p);x=x0+alpha*p;r2=r1-alpha*A*p;z2=iM*r2;belta=dot(r2,z2)/dot(r1,z1);p=z2+belta*p;n=n+1;tol=norm(x-x0);x0=x;%更新迭代值r1=r2;z1=z2;end13.BJ块雅克比迭代法求线性方程组Ax=b的解function[x,N]=BJ(A,b,x0,d,eps,M)ifnargin==4eps=1.0e-6;M=10000;elseifnargin4errorreturnelseifnargin==5M=10000;%参数的默认值endNS=size(A);n=NS(1,1);if(sum(d)~=n)disp('分块错误！');return;endbnum=length(d);bs=ones(bnum,1);fori=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);fori=1:bnumendb=bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):endb,bs(i,1):endb);%求A的对角分块矩阵endfori=1:bnumendb=bs(i,1)+d(i,1)-1;invDB(bs(i,1):endb,bs(i,1):endb)=inv(DB(bs(i,1):endb,bs(i,1):endb));%求A的对角分块矩阵的逆矩阵endN=0;tol=1;whiletol=epsx=invDB*(DB-A)*x0+invDB*b;%由于LB+DB=DB-AN=N+1;%迭代步数tol=norm(x-x0);%前后两步迭代结果的误差x0=x;if(N=M)disp('Warning:迭代次数太多，可能不收敛！');return;endend14.BGS块高斯-赛德尔迭代法求线性方程组Ax=b的解function[x,N]=BGS(A,b,x0,d,eps,M)ifnargin==4eps=1.0e-6;M=10000;elseifnargin4errorreturnelseifnargin==5M=10000;endNS=size(A);n=NS(1,1);bnum=length(d);bs=ones(bnum,1);fori=1:(bnum-1)bs(i+1,1)=sum(d(1:i))+1;%获得对角线上每个分块矩阵元素索引的起始值endDB=zeros(n,n);fori=1:bnumendb=bs(i,1)+d(i,1)-1;DB(bs(i,1):endb,bs(i,1):endb)=A(bs(i,1):end