线性方程求解.doc

线性方程组求解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 amaxk%交换两行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 amaxk%交换两行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=bif 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.00471.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 a=9 -36 30 ;-36 192 -180;30 -180 180; b=1 1 1; x=Chol_Solve(a,b)x =1.83331.08330.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:nd(k)=A(k,k);for j=1:k-1d(k)=d(k)-L(k,j)*T(k,j);endif abs(d(k) x=LDL_Solve(a,b)x =1.83331.08330.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=epsx0=x;x=B*x0+f;n=n+1;if(n=M)disp(Warning:迭代次数太多，可能不收敛！);return;endendExample: x,n=Jacobi(A,b,x0)x =0.9739-0.00471.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=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 nargin4errorreturnelseif nargin =5M= 200;endif(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 nargin4errorreturnelseif nargin =5M= 200;endif(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=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为病态方程组时，共轭梯度法收敛很慢。预处理技术是在用共轭梯度法求解之前对系数矩阵做一些变换后再求解。Example:A=25 -300 1050 -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.76672.91671.93101.43331.1349flag =0re =5.7305e-012it =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);for(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);endExample: 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.15380.76920.7692-0.15380.4615快速求解法：通用求解线性方程组的函数：x=linsolve(A,b,options)其意义为快速求解方程组Ax=b，其中A之结构由决定，内容如下表：options说