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

1、线性方程组求解Matlab程序线性方程组求解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 %回代 fo

2、r j=k+1:n b(k)=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 k=1:n-1 amax

3、=0;% 选主元 for i=k:n if abs(a(i,k)amax amax=abs(a(i,k);r=i; end end if amaxk %交换两行 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*b(k); end det=det*a(k,k);enddet=det*a(n,n); for

4、 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;% 选主元 for i=k:n if abs(a(i,k)amax amax=abs(a(i,k);r=i; end end

5、 if amaxk %交换两行 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 end end for i=1:n x(i)=b(i); endExample: x=GaussJacobi(A,

6、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)*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

7、)=(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,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

8、(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-1 delta=delta-L(k,j)2; end if delta a=9 -36 30 ;-36 192 -180;30 -180; b=1 1 1; x=Chol_Solve(a,b) x = 1.8333 1.0833

9、0.7833 对称正定矩阵之LDL分解法：function L,D=LDL_Factor(A)%对称正定矩阵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) x=LDL_Solve(a,b) x = 1.8333 1.0833 0.78332.迭代法Richardson迭代法：function x,n=richason(A,b,x0,eps,M)%Richardson

10、法求解线性方程组 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; n = n + 1; if(n=M) disp(Warning: 迭代次数太多，现在退出！); return; endend E

11、xample: 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,varargin)if nargin=3 eps= 1.0e-6; M = 200;elseif nargin=eps x0=x; x=B*x0+f; n=n+1; if(n=M) disp(Warning: 迭代次数太多，可

12、能不收敛！); 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=eps x0=x; x=G*x0+f; n=n+1; if(n=M) disp(Warning: 迭代次数太多，可能不收敛！); return; endendExample: x,n=

13、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;elseif nargin4 error returnelseif nargin =5 M = 200;end if(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

14、=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(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 nargin4 error returnels

15、eif nargin =5 M = 200;end if(w=2) error; return;end D=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)=eps x0=x; x12=B1*x0+f1; x =B2*x12+

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

17、; return; endend Example: x,n=twostep(A,b,x0) x = 0.9739 -0.0047 1.0010 n = 3 最速下降法：function 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:

18、 x,n=fastdown(A,b,x0) x = 0.9739 -0.0047 1.0010 n = 5 共轭梯度法：function x,n=conjgrad(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,r

19、1)/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;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 105

20、0 -1400 630; -300 4800 -18900 26880 -12600; 1050 -18900 79380 -117600 56700; -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

21、1.1349 flag = 0 re = 5.7305e-012 it = 10 其他迭代法：函数 说明 x=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);fo

22、r(i=1:n) if(A(i,i)=0) disp(Error: 对角有元素为0！); 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(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.