解线性方程组的列主元素高斯消去法和LU分解法

...页脚.数值试验报告分析一、实验名称：解线性程组的列主元素高斯消去法和LU分解法二、实验目的及要求：通过数值实验，从中体会解线性程组选主元的必要性和LU分解法的优点，以及程组系数矩阵和右端向量的微小变化对解向量的影响。三、算法描述：本次试验采用的是高斯列主元消去法和LU分解法求解线性程组的解。其中，高斯消去法的基本思想是避免接近于零的数作分母；能进行到底的条件:当A可逆时，列主元Gauss(高斯）消去法一定能进行到底。优点:具有很好的数值稳定性；具有与顺序Gauss消去法相同的计算量。列主元Gauss(高斯）消去法的精度显著高于顺序Gauss(高斯）消去法。注意:省去换列的步骤，每次仅选一列中最大的元。矩阵的三角分解法是A=LU,L是下三角阵，U是上三角阵,Doolittle分解:L是单位下三角阵，U是上三角阵；Crout分解:L是下三角阵，U是单位上三角阵。矩阵三角分解的条件是矩阵A有唯一的Doolittle分解的充要条件是A的前n-1顺序主子式非零；矩阵A有唯一的Crout分解的充要条件是A的前n-1顺序主子式非零。三角分解的实现是通过(1)Doolittle分解的实现；(2)Doolittle分解的缺点：条件苛刻，且不具有数值稳定性。（3）用Doolittle分解求解程组:AX=bLUX=bLY=bA=LUUX=Y；四、实验容：解下列两个线性程组（1）11134.981.4987.023.116.427.199.103.601.3321xxx（2）15900001.582012151526099999.23107104321xxxxa、用你熟悉的算法语言编写程序用列主元高斯消去法和LU分解求解上述两个程组，输出...页脚.Ax=b中矩阵A及向量b,A=LU分解的L及U，detA及解向量x.b、将程组（1）中系数3.01改为3.00，0.987改为0.990，用列主元高斯消去法求解变换后的程组，输出列主元行交换次序，解向量x及detA，并与（1）中结果比较。c、将程组（2）中的2.099999改为2.1，5.900001改为5.9，用列主元高斯消去法求解变换后的程组，输出解向量x及detA，并与（1）中的结果比较。d、用MATLAB的部函数inv求出系数矩阵的逆矩阵，再输入命令x=inv(A)*b,即可求出上述各个程组的解，并与列主元高斯消去法和LU分解法求出的解进行比较，体会选主元的法具有良好的数值稳定性。用MATLAB的部函数det求出系数行列式的值，并与（1）、（2）、（3）中输出的系数行列式的值进行比较。五、实验结果：a、用你熟悉的算法语言编写程序用列主元高斯消去法和LU分解求解上述两个程组，输出Ax=b中矩阵A及向量b,A=LU分解的L及U，detA及解向量x.Gauss(高斯）列主元法程序及结果（１）.A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34];b=[1;1;1];[RA,RB,n,X]=liezhuY(A,b)D=-0.0305RA=3RB=3n=3X=1.5926-0.6319-0.4936（２）．A=[10-701;-32.09999962;5-15-1;2102];b=[8;5.900001;5;1];[RA,RB,n,X]=liezhuY(A,b)D=-762.0001RA=4RB=4n=4X=0.0000-1.00001.00001.0000ＬＵ分解法程序及结果（１）A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34];h1=zhjLU(A)...页脚.D=-0.0305RA=3U=3.01006.03001.990004.1600-2.0696005.3107L=1.0000000.42191.000000.3279-1.63161.0000h1=3.01004.8635-0.0305（２）A=[10-701;-32.09999962;5-15-1;2102];h1=zhjLU(A)D=-762.0001RA=4U=10.0000-7.000001.000002.10006.00002.300000-2.1429-4.23810-0.0000012.7333L=1.0000000-0.30001.0000000.50001.19051.0000-0.00000.20001.14293.20001.0000h1=10.0000-0.0000-150.0001-762.0001b、将程组（1）中系数3.01改为3.00，0.987改为0.990，用列主元高斯消去法求解变换后的程组，输出列主元行交换次序，解向量x及detA，并与（1）中结果比较。线性程组数据改变后程序及结果A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34];b=[1;1;1];A(1,1)=3;A(1,3)=0.990;[RA,RB,n,X]=liezhuY(A,b)D=9.8547RA=3RB=3n=3X=-4.0264...页脚.1.91931.5210X=[1.5926-0.6319-0.4936]';X1=[-4.02641.91931.5210]';wucha=X1-Xwucha=-5.61902.55122.0146c、将程组（2）中的2.099999改为2.1，5.900001改为5.9，用列主元高斯消去法求解变换后的程组，输出解向量x及detA，并与（1）中的结果比较。线性程组数据改变后的程序及结果A=[10-701;-32.09999962;5-15-1;2102];A(2,2)=2.1;b(2,1)=5.9;b=[8;5.900001;5;1];[RA,RB,n,X]=liezhuY(A,b)D=-762.0000RA=4RB=4n=4X=-0.0000-1.00001.00001.0000X=[0;-1;1;1]';X1=[0;-1;1;1]';wucha=X1-Xwucha=0000d、用MATLAB的部函数inv求出系数矩阵的逆矩阵，再输入命令x=inv(A)*b,即可求出上述各个程组的解，并与列主元高斯消去法和LU分解法求出的解进行比较，体会选主元的法具有良好的数值稳定性。用MATLAB的部函数det求出系数行列式的值，并与（1）、（2）、（3）中输出的系数行列式的值进行比较。线性程组求解的程序及结果（1）A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34];B=inv(A)B=-1.07832.15710.51380.4281-0.8560-0.20390.3344-0.6688-0.1592b=[1;1;1];X=inv(A)*bX=1.5926...页脚.-0.6319-0.4936A=[3.016.031.99;1.274.16-1.23;0.987-4.819.34];A(1,1)=3;A(1,3)=0.990;B=inv(A)B=3.3424-6.1983-1.1705-1.32692.74420.5020-1.03652.06820.4893b=[1;1;1];X=inv(A)*bX=-4.02641.91931.5210（2）A=[10-701;-32.09999962;5-15-1;2102];B=inv(A)B=-0.0223-0.09840.11810.1686-0.1601-0.11810.14170.26900.01080.10630.0724-0.07550.10240.1575-0.18900.1969b=[8;5.900001;5;1];X=inv(A)*bX=0-1.00001.00001.0000A=[10-701;-32.09999962;5-15-1;2102];A(2,2)=2.1;B=inv(A)B=-0.0223-0.09840.11810.1686-0.1601-0.11810.14170.26900.01080.10630.0724-0.07550.10240.1575-0.18900.1969b=[8;5.900001;5;1];b(2,1)=5.9;X=inv(A)*bX=-0.0000-1.00001.00001.0000六、实验结果分析通过本次的解线性程组数值分析实验，不仅让我熟悉了Matlab软件的功能和它强大的计算能力，而且大大提高了我解决数值分析实际问题的能力，进一步掌握了列主元高斯消去...页脚.法解程组的思想和有效性。通过对以上程的计算和比较，程组（1）具有病态性。当系数矩阵有微小的变化时，所得到的解与原程组的解有很大的相对误差。程组（2）中当系数矩阵A和b有较小变化时，所得到的解没有相对误差。所以程组（2）是良性的。采用Matlab部函数inv求解逆矩阵，然后通过X=inv（a)*b也可以求解程组的解，但没有高斯列主元消去法稳定性高。Det函数求解程组系数矩阵的行列式所得到的结果和三角法所得到的结果相同，有便快捷明了的特点。七、附录列主元消去法function[RA,RB,n,X]=liezhuY(A,b)B=[Ab];n=length(b);RA=rank(A);D=det(A)RB=rank(B);zhica=RB-RA;ifzhica0,returnwarningoffMATLAB:return_outside_of_loopendifRA==RBifRA==nX=zeros(n,1);C=zeros(1,1+n);forp=1:n-1[Y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;fork=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);forq=n-1:-1:1...页脚.X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelseendendLU分解法functionh1=zhjLU(A)[nn]=size(A);RA=rank(A);D=det(A)ifRA~=nreturnendifRA==nforp=1:nh(p)=det(A(1:p,1:p));endh1=h(1:n);fori=1:nifh(1,i)==0returnendendifh(1,i)~=0forj=1:nU(1,j)=A(1,j);endfork=2:nfori=2:nforj=2:n...页脚.L(1,1)=1;L(i,i)=1;ifijL(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);elseU(k,j)=A(k,j)-L(k,1:k-1)*U(1:k-1,j);endendendendh1;RA,U,Lendend