线性方程组求解matlab程序.docx

线性方程组求解 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); det=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); end Example: 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 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 amax1e-10 return; end if rk %交换两行 for j=k:n z=a(k,j);a(k,j)=a(r,j);a(r,j)=z; z=b(k);b(k)=b(r);b(r)=z;det=-det; 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*b(k); end det=det*a(k,k); end det=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); end Example: x=detGauss(A,b) 0.9739 -0.0047 1.0010 Gauss-Jorda n 消去法： function x=GaussJacobi(a,b) % Gauss-Jacobi 消去法 [n,m]=size(a); nb=length(b); x=zeros(n,1); for k=1:n amax=0;% 选主元 for i=k:n if abs(a(i,k))amax amax=abs(a(i,k));r=i; end end if amax1e-10 return; if rk %交换两行 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); end b(i)=b(i)-a(i,k)*b(k); end for i=1:n x(i)=b(i); end Example: x=GaussJacobi(A,b) x = 0.9739 -0.0047 1.0010 LU 分解法： 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); end for 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)*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 %%%% function x=lusolv(a,b) %LU 分解求解线性方程组 aX=b if length(a)~=length(b) error(Error in inputi