用matlab解线性方程组

1、仅供个人参考用matlab解线性方程组2008-04-1217:00一。高斯消去法1 .顺序高斯消去法直接编写命令文件a=d=n,n=size(a);c=n+1a(：,c)=d;fork=1:n-1a(k+1:n,k:c)=a(k+1:n,k:c)-(a(k+1:n,k)/a(k,k)*a(k,k:c);%消去endx=0000%回带x(n)=a(n,c)/a(n,n);forg=n-1:-1:1x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n)/a(g,g)end2 .列主高斯消去法*由于r,m=max(abs(a(k:n,k)”返回的行是“k:n,k”内的第几行，所以要通过修

2、正来把m改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。直接编写命令文件a=d=n,n=size(a);c=n+1a(:,c尸d;%(增广)fork=1:n-1r,m=max(abs(a(k:n,k);%选主m=m+k-1;%(修正操作行的值)if(a(m,k)=0)if(m=k)a(km,:)=a(mk,:);%换行enda(k+1:n,k:c)=a(k+1:n,k:c)-(a(k+1:n,k)/a(k,k)*a(k,k:c);%消去endendx=0000%回带x(n)=a(n,c)/a(n,n);forg=n-1:-1:1x(g)=(a(g,c)

3、-a(g,g+1:n)*x(g+1:n)/a(g,g)end不得用于商业用途仅供个人参考不得用于商业用途3 .分别用顺序高斯消去法和列主高斯消去法解方程组a*x=d,并比较结果a=0123;9112334;62.523.415.517.2;192.0112425.159.3d=1;1;1;1顺序高斯消去法：提示“Warning:Dividebyzero.”x=NaNNaNNaNNaN列主高斯消去法：x=-1.24602.87965.5206-4.3069由此可见列主高斯消去法可以解决顺序高斯消去法所不能解决的问题。4 .将上述矩阵中的“2”改为2.05，“34”改为“34.6”，“15.5”改

4、为“15.57”，“124”改为“124.7”再用列主高斯消去法计算，与上述结果比较。x=-0.80811.82263.5568-2.7047很明显虽然系数矩阵只有很小的变化但结果的变化却很大，所以系数矩阵是病态的。这是因为系数矩阵的条件数很大：cond(a)=2.0695e+003二。迭代法J迭代公式a=521;-142;2-310d=-12;20;3x=0;0;0;%stop=1.0e-4%L=-tril(a,-1)U=-triu(a,1)D=inv(diag(diag(a)X=D*(L+U)*x+D*d;%J初始向量迭代的精度迭代公式n=1;whilenorm(X-x,inf)=stop

5、x=X;X=D*(L+U)*x+D*d;n=n+1;endXnG-S迭代公式a=521;-142;2-310x=0;0;0;%d=-12;20;3stop=1.0e-4L=-tril(a,-1)U=-triu(a,1)D=(diag(diag(a)X=inv(D-L)*U*x+inv(D-L)*d;n=1;whilenorm(X-x,inf)=stop%时迭代中止否则继续初始向量%G-S迭代公式%时迭代中止否则继续x=X;X=inv(D-L)*U*x+inv(D-L)*d;n=n+1;endXnSO砒代公式a=521;-142;2-310x=0;0;0;%初始向量d=-12;20;3stop=

6、1.0e-4w=1.44;%松弛因子L=-tril(a,-1)U=-triu(a,1)D=(diag(diag(a)X=inv(D-w*L)*(1-w)*D+w*U)*x+w*inv(D-w*L)*d%SOR迭代公式n=1;whilenorm(X-x,inf)=stop%时迭代中止否则继续x=X;X=inv(D-w*L)*(1-w)*D+w*U)*x+w*inv(D-w*L)*d;n=n+1;endXn3.a*x=da=521;-142;2-310d=-12;20;3(1)考察用J、G-S及SORB此方程组的收敛性.(2) 分别用雅可比迭代法、高斯-赛德尔迭代法及逐次超松弛迭代法解此方程组。要

7、求时迭代中止，观察收敛速度。(1) J迭代矩阵的谱半径：max(abs(eig(D*(L+U)=0.506G-S迭代矩阵的谱半径：max(abs(eig(inv(D-L)*U)=0.2000SOR4代矩阵的谱半径：max(abs(eig(inv(D-w*L)*(1-w)*D+w*U)=0.9756所以用J、G-S及SOR匀收敛(且有)。(2)J：X=-4.00003.00002.0000n=18G-S：X=-4.00003.00002.0000n=8SOR()：X=-4.00003.00002.0000n=482可见高斯-赛德尔迭代法要比雅可比迭代公式的收敛速度快。同时，如果超松弛迭代法的选取

8、不当不但不会加速收敛还会对收敛速度造成很恶劣的结果。线性方程组求解1.直接法Gauss消元法：functionx=DelGauss(a,b)%Gauss消去法n,m=size(a);nb=length(b);det=1;%存储行列式值x=zeros(n,1);fork=1:n-1fori=k+1:nifa(k,k)=0returnendm=a(i,k)/a(k,k);forj=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);fork=n:-1:1%回代forj=k

9、+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);endExample: A=1.0170-0.00920.0095;-0.00920.99030.0136;0.00950.01360.9898; b=101; x=DelGauss(A,b)0.9739-0.0047列主元Gauss消去法：functionx=detGauss(a,b)%Gauss列主元消去法n,m=size(a);nb=length(b);det=1;%存储行列式值x=zeros(n,1);fork=1:n-1amax=0;%选主元fori=k:nifabs(a(i,k)amaxam

10、ax=abs(a(i,k);r=i;endendifamaxk%交换两行forj=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;endfori=k+1:n%进行消元m=a(i,k)/a(k,k);forj=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);fork=n:-1:1%回代forj=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k)

11、;endExample:x=detGauss(A,b)x=0.9739-0.00471.0010Gauss-Jordan消去法：functionx=GaussJacobi(a,b)%Gauss-Jacobi消去法n,m=size(a);nb=length(b);x=zeros(n,1);fork=1:namax=0;%选主元fori=k:nifabs(a(i,k)amaxamax=abs(a(i,k);r=i;endendifamaxk%交换两行forj=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%进行消元

12、b(k)=b(k)/a(k,k);forj=k+1:na(k,j)=a(k,j)/a(k,k);endfori=1:nifi=kforj=k+1:na(i,j)=a(i,j)-a(i,k)*a(k,j);endb(i)=b(i)-a(i,k)*b(k);endendendfori=1:nx(i)=b(i);endExample:x=GaussJacobi(A,b)x=0.9739-0.00471.0010LU分解法：functionl,u=lu(a)%LU分解n=length(a);l=eye(n);u=zeros(n);fori=1:nu(1,i)=a(1,i);endfori=2:nl(i

13、,1)=a(i,1)/u(1,1);endforr=2:n%fori=r:nuu=0;fork=1:r-1uu=uu+l(r,k)*u(k,i);endu(r,i)=a(r,i)-uu;end%fori=r+1:nll=0;fork=1:r-1ll=ll+l(i,k)*u(k,r);endl(i,r)=(a(i,r)-ll)/u(r,r);end%Endfunctionx=lusolv(a,b)%LU分解求解线性方程组aX=biflength(a)=length(b)error(Errorininputing!)return;endn=length(a);l,u=lu(a);y(1)=b(1)

14、;fori=2:nz=0;fork=1:i-1z=z+l(i,k)*y(k);endy(i)=b(i)-z;endx(n)=y(n)/u(n,n);fori=n-1:-1:1z=0;fork=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分解法：functionL=Cholesky(A)%对对称正定矩阵A进行Cholesky分解n=length(A);L=zeros(n);fork=1:ndelta=A(k,k);forj=1:k-

15、1delta=delta-L(k,j42;endifdelta a=9-3630;-36192-180;30-180180; b=111; x=Chol_Solve(a,b)x=1.83331.08330.7833对称正定矩阵之LDL分解法：functionL,D=LDL_Factor(A)%对称正定矩阵A进行LDL分解n=length(A);L=eye(n);D=zeros(n);d=zeros(1,n);T=zeros(n);fork=1:nd(k)=A(k,k);forj=1:k-1d(k)=d(k)-L(k,j)*T(k,j);endifabs(d(k)x=LDL_Solve(a,b)

16、x=1.83331.08330.78332.迭代法Richardson迭代法：functionx,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;

17、n=n+1;if(n=M)disp(Warning:迭代次数太多，现在退出！);return;endendExample:A=1.0170-0.00920.0095;-0.00920.99030.0136;0.00950.01360.9898;b=101;x0=000;x,n=richason(A,b,x0)不得用于商业用途仅供个人参考x=0.9739-0.00471.0010n=5Jacobi迭代法：functionx,n=jacobi(A,b,x0,eps,varargin)ifnargin=3eps=1.0e-6;M=200;elseifnargin=epsx0=x;x=B*x0+f;n

18、=n+1;if(n=M)disp(Warning:迭代次数太多，可能不收敛！);return;endendExample:x,n=Jacobi(A,b,x0)x=0.9739-0.00471.0010n=5Gauss-Seidel迭代法：functionx,n=gauseidel(A,b,x0,eps,M)ifnargin=3eps=1.0e-6;M=200;elseifnargin=4M=200;elseifnargin=epsx0=x;x=G*x0+f;n=n+1;if(n=M)disp（Warning:迭代次数太多，可能不收敛！）;return;endendExample:x,n=gau

19、seidel(A,b,x0)0.9739-0.00471.0010超松驰迭代法：functionx,n=SOR(A,b,x0,w,eps,M)ifnargin=4eps=1.0e-6;M=200;elseifnargin4errorreturnelseifnargin=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;wh

20、ilenorm(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对称逐次超松驰迭代法：functionx,n=SSOR(A,b,x0,w,eps,M)ifnargin=4eps=1.0e-6;M=200;elseifnargin4errorreturnelseifnargin=5仅供个人参考M=200;endif(w=2)error;return;end不得用于商业用途D=diag(diag(A)

21、;%求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;whilenorm(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

22、,x0,1)x=0.9739-0.00471.0010n=两步迭代法：functionx,n=twostep(A,b,x0,eps,varargin)ifnargin=3eps=1.0e-6;M=200;elseifnargin=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.0047n=3最速下降法：functionx,n=fastdown(A,b,x0,eps)if(nargin=3)eps

23、=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共轭梯度法：functionx,n=conjgrad(A,b,x0)if(nargin=3)eps=1.0e-6;r1=b-A*x0;p1=r1;d=dot(r1,r1)/dot(p1,A*p1);x=x0+d*p1;r2=r1-d*A*

24、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;end不得用于商业用途仅供个人参考d=dot(r2,r2)/dot(p2,A*p2);x=x+d*p2;n=n+1;Example:x,n=conjgrad(A,b,x0)0.9739-0.00471.0010预处理的共轭梯度法：当AX=B为病态方程组时，共轭梯度

25、法收敛很慢。预处理技术是在用共轭梯度法求解之前对系数矩阵做一些变换后再求解。Example:A=25-3001050-1400630;-3004800-1890026880-12600;1050-1890079380-11760056700;-140026880-117600179200-88200;630-1260056700-8820044100;b=53-10-2;x0=00000;M=pascal(5)%预处理矩阵x,flag,re,it=pcg(A,b,1.e-8,1000,M,M,x0)%flag=0表示在指定迭代次数之内按要求精度收敛%re表示相对误差%it表示迭代次数x=5.7

26、6672.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特殊解法解三对角线性方程组之追赶法：functionx=followup(A,b)n=rank(A);for(i=1:n)

27、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

28、)=(b(i,1)-c(i,1)*x(i+1,1)/d(i,1);endExample:A=2.51.00000;1 1.51.0000;01.00.51.000;001.00.51.00;0001.01.51.0;00001.02.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说明LT上三角阵UT下三角阵UHESS上三角HessenbergSYM实对称矩阵POSDEF正定矩阵RECT一般矩阵TRANSA指出求解的方程是Ax=b还是Ax=b,A为之共轲转置Example: A=4-12;-150;206;b=1-12; optt.SYM=true