线性方程组求解1.直接法Gauss消元法：function x=DelGauss(a,b)% Gauss消去法[n,m]=size(a);nb=length(b);det=1;%存储行列式值x=zeros(n,1);for k=1:n-1for i=k+1:nif a(k,k)==0returnendm=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1 %回代for j=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);endExample:>> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]';>> x=DelGauss(A,b)x =0.9739-0.00471.0010列主元Gauss消去法：function x=detGauss(a,b)% Gauss列主元消去法[n,m]=size(a);nb=length(b);det=1;%存储行列式值x=zeros(n,1);for k=1:n-1amax=0;% 选主元for i=k:nif abs(a(i,k))>amaxamax=abs(a(i,k));r=i;endendif amax<1e-10return;endif r>k %交换两行for j=k:nz=a(k,j);a(k,j)=a(r,j);a(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:n %进行消元m=a(i,k)/a(k,k);for j=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);for k=n:-1:1 %回代for j=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);endExample:>> x=detGauss(A,b)x =0.9739-0.00471.0010Gauss-Jordan消去法：function x=GaussJacobi(a,b)% Gauss-Jacobi消去法[n,m]=size(a);nb=length(b);x=zeros(n,1);for k=1:namax=0;% 选主元for i=k:nif abs(a(i,k))>amaxamax=abs(a(i,k));r=i;endendif amax<1e-10return;endif r>k %交换两行for j=k:nz=a(k,j);a(k,j)=a(r,j);a(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;end%进行消元b(k)=b(k)/a(k,k);for j=k+1:na(k,j)=a(k,j)/a(k,k);endfor i=1:nif i~=kfor j=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);endb(i)=b(i)-a(i,k)*b(k);endendendfor i=1:nx(i)=b(i);endExample:>> x=GaussJacobi(A,b) x =0.9739-0.00471.0010LU分解法：function [l,u]=lu(a)%LU分解n=length(a);l=eye(n);u=zeros(n);for i=1:nu(1,i)=a(1,i);endfor i=2:nl(i,1)=a(i,1)/u(1,1);endfor r=2:n%%%%for i=r:nuu=0;for k=1:r-1uu=uu+l(r,k)*u(k,i);endu(r,i)=a(r,i)-uu;end%%%%for i=r+1:nll=0;for k=1:r-1ll=ll+l(i,k)*u(k,r);endl(i,r)=(a(i,r)-ll)/u(r,r);end%%%%Endfunction x=lusolv(a,b)%LU分解求解线性方程组aX=b if length(a)~=length(b)error('Error in inputing!')return;endn=length(a);[l,u]=lu(a);y(1)=b(1);for i=2:nz=0;for k=1:i-1z=z+l(i,k)*y(k);endy(i)=b(i)-z;endx(n)=y(n)/u(n,n);for i=n-1:-1:1z=0;for k=i+1:nz=z+u(i,k)*x(k);endx(i)=(y(i)-z)/u(i,i);endExample:>> x=lusolv(A,b)x =0.9739 -0.0047 1.0010对称正定矩阵之Cholesky分解法：function L=Cholesky(A)%对对称正定矩阵A进行Cholesky分解n=length(A);L=zeros(n);for k=1:ndelta=A(k,k);for j=1:k-1delta=delta-L(k,j)^2;endif delta<1e-10return;endL(k,k)=sqrt(delta);for i=k+1:nL(i,k)=A(i,k);for j=1:k-1L(i,k)=L(i,k)-L(i,j)*L(k,j);endL(i,k)=L(i,k)/L(k,k);endendfunction x=Chol_Solve(A,b)%利用对称正定矩阵之Cholesky分解求解线性方程组Ax=b n=length(b);l=Cholesky(A);x=ones(1,n);y=ones(1,n);for i=1:nz=0;for k=1:i-1z=z+l(i,k)*y(k);endy(i)=(b(i)-z)/l(i,i);endfor i=n:-1:1for k=i+1:nz=z+l(k,i)*x(k);endx(i)=(y(i)-z)/l(i,i);endExample:>> a=[9 -36 30 ;-36 192 -180;30 -180 180]; >> b=[1 1 1]';>> x=Chol_Solve(a,b)x =1.8333 1.0833 0.7833对称正定矩阵之LDL’分解法：function [L,D]=LDL_Factor(A)%对称正定矩阵A进行LDL'分解n=length(A);L=eye(n);D=zeros(n);d=zeros(1,n);for k=1:nd(k)=A(k,k);for j=1:k-1d(k)=d(k)-L(k,j)*T(k,j);endif abs(d(k))<1e-10return;endfor i=k+1:nT(i,k)=A(i,k);for j=1:k-1T(i,k)=T(i,k)-T(i,j)*L(k,j);endL(i,k)=T(i,k)/d(k);endendD=diag(d);function x=LDL_Solve(A,b)%利用对对称正定矩阵A进行LDL'分解法求解线性方程组Ax=b n=length(b);[l,d]=LDL_Factor(A);y(1)=b(1);for i=2:nz=0;for k=1:i-1z=z+l(i,k)*y(k);endy(i)=b(i)-z;endx(n)=y(n)/d(n,n);for i=n-1:-1:1z=0;for k=i+1:nz=z+l(k,i)*x(k);endx(i)=y(i)/d(i,i)-z; endExample:>> x=LDL_Solve(a,b)x =1.8333 1.0833 0.78332.迭代法Richardson迭代法：function [x,n]=richason(A,b,x0,eps,M) %Richardson法求解线性方程组Ax=b %方程组系数矩阵:A%方程组之常数向量:b%迭代初始向量:X0%e解的精度控制:eps%迭代步数控制：M%返回值线性方程组的解：x%返回值迭代步数：nif(nargin == 3)eps = 1.0e-6;M = 200;elseif(nargin == 4)M = 200;endI =eye(size(A));x1=x0;x=(I-A)*x0+b;n=1;while(norm(x-x1)>eps)x1=x;x=(I-A)*x1+b;n = n + 1;if(n>=M)disp('Warning: 迭代次数太多，现在退出！');return;endendExample:>> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]';x0=[0 0 0]';>> [x,n]=richason(A,b,x0)x =0.9739-0.00471.0010n =5Jacobi迭代法：function [x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3eps= 1.0e-6;M = 200;elseif nargin<3errorreturnelseif nargin ==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; %迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=Jacobi(A,b,x0)x =0.97391.0010n =5Gauss-Seidel迭代法：function [x,n]=gauseidel(A,b,x0,eps,M) if nargin==3eps= 1.0e-6;M = 200;elseif nargin == 4M = 200;elseif nargin<3errorreturn;endD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵f=(D-L)\b;x=G*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=gauseidel(A,b,x0)x =0.9739-0.00471.0010n =4超松驰迭代法：function [x,n]=SOR(A,b,x0,w,eps,M) if nargin==4eps= 1.0e-6;M = 200;elseif nargin<4errorreturnelseif nargin ==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; %迭代次数while norm(x-x0)>=epsx0=x;x =B*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=SOR(A,b,x0,1)x =0.9739-0.00471.0010n =4对称逐次超松驰迭代法：function [x,n]=SSOR(A,b,x0,w,eps,M) if nargin==4eps= 1.0e-6;M = 200;elseif nargin<4errorreturnelseif nargin ==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; %迭代次数while norm(x-x0)>=epsx0=x;x12=B1*x0+f1;x =B2*x12+f2;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=SSOR(A,b,x0,1)x =0.9739-0.00471.0010n =3两步迭代法：function [x,n]=twostep(A,b,x0,eps,varargin)if nargin==3eps= 1.0e-6;M = 200;elseif nargin<3errorreturnelseif nargin ==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; %迭代次数while norm(x-x0)>=epsx0 =x;x12=B1*x0+f1;x =B2*x12+f2;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendExample:>> [x,n]=twostep(A,b,x0)x =0.9739-0.00471.0010n =3最速下降法：function [x,n]=fastdown(A,b,x0,eps) if(nargin == 3)eps = 1.0e-6;endr = b-A*x0;d = dot(r,r)/dot(A*r,r);x = x0+d*r;n=1;while(norm(x-x0)>eps)x0 = x;r = b-A*x0;d = dot(r,r)/dot(A*r,r);x = x0+d*r;n = n + 1;endExample:>> [x,n]=fastdown(A,b,x0)x =0.9739-0.00471.0010n =5共轭梯度法：function [x,n]=conjgrad(A,b,x0) if(nargin == 3)eps = 1.0e-6;endr1 = b-A*x0;p1 = r1;d = dot(r1,r1)/dot(p1,A*p1);x = x0+d*p1;r2 = r1-d*A*p1;f = dot(r2,r2)/dot(r1,r1);p2 = r2+f*p1;n = 1;for(i=1:(rank(A)-1))x0 = x;p1 = p2;r1 = r2;d = dot(r1,r1)/dot(p1,A*p1);x = x0+d*p1;r2 = r1-d*A*p1;f = dot(r2,r2)/dot(r1,r1);p2 = r2+f*p1;n = n + 1;endd = dot(r2,r2)/dot(p2,A*p2);x = x+d*p2;n = n + 1;Example:>> [x,n]=conjgrad(A,b,x0)x =0.9739-0.00471.0010n =4预处理的共轭梯度法：当AX=B为病态方程组时，共轭梯度法收敛很慢。