线性方程组求解Matlab程序

1、线性方程组求解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-1 for i=k+1:n if a(k,k)=0 return end m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k);enddet=det*a(n,n); for k=n:-1:1 %回代 for j=k+1:n b(k)=

2、b(k)-a(k,j)*x(j); end x(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.0047 1.0010列主元Gauss消去法：function x=detGauss(a,b)% Gauss列主元消去法n,m=size(a);nb=length(b);det=1;%存储行列式值x=zeros(n,1);for

3、 k=1:n-1 amax=0;% 选主元 for i=k:n if abs(a(i,k)>amax amax=abs(a(i,k);r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n z=a(k,j);a(k,j)=a(r,j);a(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det; end for i=k+1:n %进行消元 m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*

4、b(k); end det=det*a(k,k);enddet=det*a(n,n); for k=n:-1:1 %回代 for j=k+1:n b(k)=b(k)-a(k,j)*x(j); end x(k)=b(k)/a(k,k);endExample:>> x=detGauss(A,b) x = 0.9739 -0.00471.0010 Gauss-Jordan消去法：function x=GaussJacobi(a,b)% Gauss-Jacobi消去法n,m=size(a);nb=length(b);x=zeros(n,1);for k=1:n amax=0;% 选主元 f

5、or i=k:n if abs(a(i,k)>amax amax=abs(a(i,k);r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n z=a(k,j);a(k,j)=a(r,j);a(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z; end %进行消元 b(k)=b(k)/a(k,k); for j=k+1:n a(k,j)=a(k,j)/a(k,k); end for i=1:n if i=k for j=k+1:n a(i,j)=a(i,j)-a(i,k)*a(k,j);

6、 end b(i)=b(i)-a(i,k)*b(k); end end end for i=1:n x(i)=b(i); endExample:>> x=GaussJacobi(A,b) x = 0.9739 -0.0047 1.0010LU分解法：function l,u=lu(a)%LU分解n=length(a);l=eye(n);u=zeros(n);for i=1:n u(1,i)=a(1,i);endfor i=2:n l(i,1)=a(i,1)/u(1,1);end for r=2:n % for i=r:n uu=0; for k=1:r-1 uu=uu+l(r,k)

7、*u(k,i); end u(r,i)=a(r,i)-uu; end % for i=r+1:n ll=0; for k=1:r-1 ll=ll+l(i,k)*u(k,r); end l(i,r)=(a(i,r)-ll)/u(r,r); end %Endfunction x=lusolv(a,b)%LU分解求解线性方程组aX=bif length(a)=length(b) error('Error in inputing!') return;endn=length(a);l,u=lu(a);y(1)=b(1);for i=2:n z=0; for k=1:i-1 z=z+l(i

8、,k)*y(k); end y(i)=b(i)-z;end x(n)=y(n)/u(n,n);for i=n-1:-1:1 z=0; for k=i+1:n z=z+u(i,k)*x(k); end x(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:n delta=A(k,k); for j=1:k-

9、1 delta=delta-L(k,j)2; end if delta<1e-10 return; end L(k,k)=sqrt(delta); for i=k+1:n L(i,k)=A(i,k); for j=1:k-1 L(i,k)=L(i,k)-L(i,j)*L(k,j); end L(i,k)=L(i,k)/L(k,k); endendfunction x=Chol_Solve(A,b)%利用对称正定矩阵之Cholesky分解求解线性方程组Ax=bn=length(b);l=Cholesky(A);x=ones(1,n);y=ones(1,n);for i=1:n z=0; f

10、or k=1:i-1 z=z+l(i,k)*y(k); end y(i)=(b(i)-z)/l(i,i);endfor i=n:-1:1 z=0; for k=i+1:n z=z+l(k,i)*x(k); end x(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)%对

11、称正定矩阵A进行LDL'分解n=length(A);L=eye(n);D=zeros(n);d=zeros(1,n);T=zeros(n);for k=1:n d(k)=A(k,k); for j=1:k-1 d(k)=d(k)-L(k,j)*T(k,j); end if abs(d(k)<1e-10 return; end for i=k+1:n T(i,k)=A(i,k); for j=1:k-1 T(i,k)=T(i,k)-T(i,j)*L(k,j); end L(i,k)=T(i,k)/d(k);endendD=diag(d); function x=LDL_Solve(

12、A,b)%利用对对称正定矩阵A进行LDL'分解法求解线性方程组Ax=bn=length(b);l,d=LDL_Factor(A); y(1)=b(1);for i=2:n z=0; for k=1:i-1 z=z+l(i,k)*y(k); end y(i)=b(i)-z;end x(n)=y(n)/d(n,n);for i=n-1:-1:1 z=0; for k=i+1:n z=z+l(k,i)*x(k); end x(i)=y(i)/d(i,i)-z;endExample: >> x=LDL_Solve(a,b) x = 1.8333 1.0833 0.78332.迭代法

13、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;end I =eye(size(A);x1=x0;x=(I-A)*x0+b;n=1; while(norm(x-x1)>eps) x1=x; x=(I-A)*x1+b;

14、n = n + 1; if(n>=M) disp('Warning: 迭代次数太多，现在退出！'); return; endend Example:>> 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.0047 1.0010 n = 5 Jacobi迭代法：function x,n=jacobi(A,b,x0,eps,var

15、argin)if nargin=3 eps= 1.0e-6; M = 200;elseif nargin<3 error returnelseif nargin =5 M = varargin1;end D=diag(diag(A); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=D(L+U);f=Db;x=B*x0+f;n=1; %迭代次数 while norm(x-x0)>=eps x0=x; x=B*x0+f; n=n+1; if(n>=M) disp('Warning: 迭代次数太多，可能不收敛！

16、'); return; endend Example:>> x,n=Jacobi(A,b,x0) x = 0.9739 -0.0047 1.0010 n = 5 Gauss-Seidel迭代法：function x,n=gauseidel(A,b,x0,eps,M)if nargin=3 eps= 1.0e-6; M = 200;elseif nargin = 4 M = 200;elseif nargin<3 error return;end D=diag(diag(A); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %

17、求A的上三角阵G=(D-L)U;f=(D-L)b;x=G*x0+f;n=1; %迭代次数while norm(x-x0)>=eps x0=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.0047 1.0010 n = 4 超松驰迭代法：function x,n=SOR(A,b,x0,w,eps,M)if nargin=4 eps= 1.0e-6; M = 200

18、;elseif nargin<4 error returnelseif nargin =5 M = 200;end if(w<=0 | w>=2) error; return;end D=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)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=M) disp(&

19、#39;Warning: 迭代次数太多，可能不收敛！'); return; endend Example:>> x,n=SOR(A,b,x0,1) x = 0.9739 -0.0047 1.0010 n = 4 对称逐次超松驰迭代法：function x,n=SSOR(A,b,x0,w,eps,M)if nargin=4 eps= 1.0e-6; M = 200;elseif nargin<4 error returnelseif nargin =5 M = 200;end if(w<=0 | w>=2) error; return;end D=diag(

20、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)>=eps x0=x; x12=B1*x0+f1; x =B2*x12+f2; n=n+1; if(n>=M) disp('Warning: 迭代次数太多，可能不收敛！

21、'); return; endend Example:>> x,n=SSOR(A,b,x0,1) x = 0.9739 -0.0047 1.0010 n = 3 两步迭代法：function x,n=twostep(A,b,x0,eps,varargin)if nargin=3 eps= 1.0e-6; M = 200;elseif nargin<3 error returnelseif nargin =5 M = varargin1;end D=diag(diag(A); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求

22、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)>=eps x0 =x; x12=B1*x0+f1; x =B2*x12+f2; n=n+1; if(n>=M) disp('Warning: 迭代次数太多，可能不收敛！'); return; endend Example:>> x,n=twostep(A,b,x0) x = 0.9739 -0.0047 1.0010 n = 3 最速下降法：functio

23、n x,n=fastdown(A,b,x0,eps)if(nargin = 3) eps = 1.0e-6;end r = 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;end Example:>> x,n=fastdown(A,b,x0) x = 0.9739 -0.0047 1.0010 n = 5 共轭梯度法：function x,n=conjgr

24、ad(A,b,x0)if(nargin = 3) eps = 1.0e-6;end r1 = 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 = r

25、2+f*p1; n = n + 1;end d = 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.0047 1.0010 n = 4 预处理的共轭梯度法： 当AX=B为病态方程组时，共轭梯度法收敛很慢。预处理技术是在用共轭梯度法求解之前对系数矩阵做一些变换后再求解。Example:A=25 -300 1050 -1400 630; -300 4800 -18900 26880 -12600; 1050 -18900 79380 -117600 5

26、6700; -1400 26880 -117600 179200 -88200; 630 -12600 56700 -88200 44100;b=5 3 -1 0 -2'x0=0 0 0 0 0'M=pascal(5)%预处理矩阵x,flag,re,it=pcg(A,b,1.e-8,1000,M,M,x0)%flag=0表示在指定迭代次数之内按要求精度收敛%re表示相对误差%it表示迭代次数>> x = 5.7667 2.9167 1.9310 1.4333 1.1349 flag = 0 re = 5.7305e-012 it = 10 其他迭代法：函数 说明 x

27、=symmlq(A,b) 线性方程组的LQ解法 x=bicg(A,b) 线性方程组的双共轭梯度法 x=bicgstab(A,b) 线性方程组的稳定双共轭梯度法 x=lsqr(A,b) 线性方程组的共轭梯度LSQR解法 x=gmres(A,b) 线性方程组的广义最小残差解法 x=minres(A,b) 线性方程组的最小残差解法 x=qmr(A,b) 线性方程组的准最小残差解法 3特殊解法解三对角线性方程组之追赶法：function x=followup(A,b)n = rank(A);for(i=1:n) if(A(i,i)=0) disp('Error: 对角有元素为0！');

28、 return; endend; d = ones(n,1);a = ones(n-1,1);c = ones(n-1); for(i=1:n-1) a(i,1)=A(i+1,i); c(i,1)=A(i,i+1); d(i,1)=A(i,i);endd(n,1) = A(n,n); for(i=2:n) d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1)*c(i-1,1); b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1)*b(i-1,1);endx(n,1) = b(n,1)/d(n,1); for(i=(n-1):-1:1) x(i,1) = (b(

29、i,1)-c(i,1)*x(i+1,1)/d(i,1);end Example:>> A=2.5 1.0 0 0 0 0; 1 1.5 1.0 0 0 0; 0 1.0 0.5 1.0 0 0; 0 0 1.0 0.5 1.0 0; 0 0 0 1.0 1.5 1.0; 0 0 0 0 1.0 2.5;>> b=ones(6,1);>> x=followup(A,b) x = 0.4615 -0.1538 0.7692 0.7692 -0.1538 0.4615 快速求解法：通用求解线性方程组的函数：x=linsolve(A,b,options)其意义为快速求解方程组Ax=b，其中