解线性方程组的直接方法和迭代法.doc
解线性方程组的直接方法和迭代法用LU分解及列主元高斯消去法解线性方程组 .输出中系数分解法的矩阵及，解向量及；列主元法的行交换次序，解向量及；比较两种方法所得的结果.2线性方程组的及为 则解.用Matlab内部函数求及的所有特征值和若令 ，求解，输出向量和，从理论结果和实际计算两方面分析线性方程组解的相对误差及的相对误差的关系。 3给出线性方程组，其中系数矩阵为希尔伯特矩阵： ，假设，若取，分别用雅可比迭代法及SOR迭代法求解，比较计算结果。结果1. >>A=10 -7 0 1;-3 2. 6 2;5 -1 5 -1;2 1 0 2>> L,U,y,x=zhjLU(A,b)结果：L = 1.0e+006 * 0.0000 0 0 0 -0.0000 0.0000 0 0 0.0000 -2.5000 0.0000 0 0.0000 -2.4000 0.0000 0.0000U = 1.0e+007 * 0.0000 -0.0000 0 0.0000 0 -0.0000 0.0000 0.0000 0 0 1.5000 0.5750 0 -0.0000 0 0.0000y = 1.0e+007 * 0.0000 0.0000 2.0750 0.0000x = 0.0000 -1.0000 1.0000 1.0000列主元消去法：>> b=8 5.9 5 1>> b=b'>> RA,RB,n,X=liezhu(A,b)结果：RA =4RB =4n =4X = 0.0000 -1.0000 1.0000 1.00002. >> det(A)ans =1.0000>> eig(A)ans = 0.0102 0.8431 3.8581 30.2887用高斯列主元消去法求向量>> A=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98>> b=32;23;33;31结果为：X = -9.5863 18.3741 -3.22583.5240>> X-1ans = -10.5863 17.3741 -4.2258 2.5240求得：>> norm(ans)ans =20.9322求得：>> ans/norm(X)ans = 0.9842求得：>> norm(A)/norm(a)ans =1.00203.生成6阶希尔伯特矩阵以及对应的b>> Hilbert(6)H = 1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909b = 2.4500 1.5929 1.2179 0.9956 0.84560.7365利用雅克比迭代计算：>> jacdd(H,b,2,1,4,3,6,7,2,0.,10000)误差为时，超松弛迭代法结果为松弛因子迭代次数近似解1.0620(1.0005,1.0045,0.9626,1.0441,1.0285,0.9583)1.25588(0.9997,1.0134,0.9362,1.0706,1.0230,0.9555)1.5539(0.9991,1.0211,0.9082,1.1146,0.9899,0.9656)>> x1,n1=SOR(H,b,0;0;0;0;0;0,1.5)>> x2,n2=SOR(H,b,0;0;0;0;0;0,1.25)>> x3,n3=SOR(H,b,0;0;0;0;0;0,1)误差为时，超松弛迭代法结果为松弛因子迭代次数近似解1.0426(0.9970,1.0417,0.8967,1.0155,1.0654,1.0505,0.9991,0.9309)1.251157(0.9971,1.0353,0.9174,1.0167,1.0408,1.0378,1.0022,0.9508)1.51701(0.9980,1.0232,0.9484,1.0045,1.0260,1.0324,1.0037,0.9623)>> x1,n1=SOR(H,b,0;0;0;0;0;0;0;0,1)>> x2,n2=SOR(H,b,0;0;0;0;0;0;0;0,1.25)>> x3,n3=SOR(H,b,0;0;0;0;0;0;0;0,1.5)源程序：1. LU分解function L,U,y,x=zhjLU(A,b)n n =size(A); RA=rank(A); if RA=ndisp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'), RA,hl=det(A);returnendif RA=n for p=1:n h(p)=det(A(1:p, 1:p); end hl=h(1:n); for i=1:n if h(1,i)=0 disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'), hl;RA return end end if h(1,i)=0 disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.') for j=1:n U(1,j)=A(1,j); end for k=1:n for i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if i>j L(1,1)=1;L(2,1)=A(2,1)/U(1,1); L(i,1)=A(i,1)/U(1,1); L(i,k)=(A(i,k)- L(i,1:k-1)*U(1:k-1,k)/U(k,k); else U(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j); end end end end hl; endend%求解x，yy=zeros(n,1);x=zeros(n,1);y(1)=b(1);for i=2:n y(i)=b(i)-L(i,1:i-1)*y(1:i-1); endx(n)=y(n)/U(n,n);for i=n-1:-1:1 x(i)=(y(i)-U(i,i+1:n)*x(i+1:n)/U(i,i);end列主元高斯消去function RA,RB,n,X=liezhu(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意：因为RA=RB，所以此方程组无解.')returnendif RA=RB if RA=ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1Y,j=max(abs(B(p:n,p); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse disp('请注意：因为RA=RB<n，所以此方程组有无穷多解.')endend3.生成希尔伯特矩阵：function Hilbert(n)H=ones(n,n);for i=1:n for j=1:n H(i,j)=1/(i+j-1); endendH x1=ones(n);x=x1(:,1);b=H*x用雅可比迭代解线性方程组的MATLAB主程序function X=jacdd(A,b,X0,P,wucha,max1)n m=size(A); for j=1:ma(j)=sum(abs(A(:,j)-2*(abs(A(j,j);endfor i=1:nif a(i)>=0disp('请注意：系数矩阵A不是严格对角占优的，此雅可比迭代不一定收敛')returnendendif a(i)<0disp('请注意：系数矩阵A是严格对角占优的，此方程组有唯一解，且雅可比迭代收敛 ')endfor k=1:max1k for j=1:mX(j)=(b(j)-A(j,1:j-1,j+1:m)*X0(1: j-1,j+1:m)/A(j,j);endX,djwcX=norm(X'-X0,P); xdwcX=djwcX/(norm(X',P)+eps); X0=X'X1=Ab;if (djwcX<wucha)&(xdwcX<wucha) disp('请注意：雅可比迭代收敛,此方程组的精确解jX和近似解X如下：')returnendendif (djwcX>wucha)&(xdwcX>wucha)disp('请注意：雅可比迭代次数已经超过最大迭代次数max1 ')enda,X=X;jX=X1',超松弛法迭代求解程序function x,n=SOR(A,b,x0,w,eps,M)if nargin=4 eps= 1.0e-4; M = 2000;elseif nargin<4 error returnelseif nargin =5 M = 2000;end if(w<=0 | 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)>=eps x0=x; x =B*x0+f; n=n+1; if(n>=M) disp('Warning: 迭代次数太多，可能不收敛！'); return; endend
