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

解线性方程组的列主元素高斯消去法和LU分解法一、实验目的：通过数值实验，从中体会解线性方程组选主元的必要性和LU分解法的优点，以及方程组系数矩阵和右端向量的微小变化对解向量的影响。二、实验内容：解下列两个线性方程组（1）11134.981.4987.023.116.427.199.103.601.3321xxx（2）15900001.582012151526099999.23107104321xxxx三、实验要求：（1）用你熟悉的算法语言编写程序用列主元高斯消去法和LU分解求解上述两个方程组，输出Ax=b中矩阵A及向量b,A=LU分解的L及U，detA及解向量x.（2）将方程组（1）中系数3.01改为3.00，0.987改为0.990，用列主元高斯消去法求解变换后的方程组，输出列主元行交换次序，解向量x及detA，并与（1）中结果比较。（3）将方程组（2）中的2.099999改为2.1，5.900001改为5.9，用列主元高斯消去法求解变换后的方程组，输出解向量x及detA，并与（1）中的结果比较。（4）用MATLAB的内部函数inv求出系数矩阵的逆矩阵，再输入命令x=inv(A)*b,即可求出上述各个方程组的解，并与列主元高斯消去法和LU分解法求出的解进行比较，体会选主元的方法具有良好的数值稳定性。用MATLAB的内部函数det求出系数行列式的值，并与（1）、（2）、（3）中输出的系数行列式的值进行比较。四、实验过程：（1）列主元高斯消去法的主程序为function[RA,RB,n,X]=liezhuY(A,b)B=[Ab];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;D=det(A)ifzhica0,disp('请注意：因为RA~=RB，所以此方程组无解.')returnendifRA==RBifRA==ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.')X=zeros(n,1);C=zeros(1,n+1);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:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RBn，所以此方程组有无穷多解.')endend解方程组（1）在MATLAB工作窗口输入A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];b=[1;1;1];[RA,RB,n,X]=liezhuY(A,b)运行后输出结果为请注意：因为RA=RB=n，所以此方程组有唯一解.D=-0.1225RA=3RB=3n=3X=397.8654-157.6242-123.1120解方程组（2）在MATLAB工作窗口输入A=[10-701;-32.09999962;5-15-1;2102];b=[8;5.900001;5;1];[RA,RB,n,X]=liezhu(A,b)运行后输出结果为请注意：因为RA=RB=n，所以此方程组有唯一解.D=-762.0000RA=4RB=4n=4X=0.0000-1.00001.00001.0000LU分解法及MATLAB主程序为functionhl=zhjLU(A)[nn]=size(A);RA=rank(A);D=det（A）ifRA~=ndisp('请注意：因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:'),RA,hl=det(A);returnendifRA==nforp=1:nh(p)=det(A(1:p,1:p));endhl=h(1:n);fori=1:nifh(1,i)==0disp('请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：'),hl;RAreturnendendifh(1,i)~=0disp('请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：')forj=1:nU(1,j)=A(1,j);endfork=2:nfori=2:nforj=2:nL(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);endendendendhl;RA,U,Lendend解方程组（1）在MATLAB工作窗口输入A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];h1=zhjLU(A)运行输出结果为请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：D=9.8547RA=3U=3.01006.03001.999004.1600-2.0734005.3016L=1.0000000.42191.000000.3279-1.63161.0000h1=3.01004.8635-0.1225解方程组（2）在MATLAB工作窗口输入A=[10-701;-32.09999962;5-15-1;2102];h1=zhjLU(A)运行后输出结果为请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：D=-762.0000RA=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.0001（2）在MATLAB工作窗口输入A=[3.016.031.999;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]=liezhu(A,b)请注意：因为RA=RB=n，所以此方程组有唯一解.RA=3RB=3n=3X=-4.02641.91931.5210hi=3.00004.82199.8547在MATLAB工作窗口输入x=[397.8654;-157.6242;-123.1120]';x1=[-4.0264;1.9193;1.5210]';wucha=x1-x运行后输出结果为wucha=-401.8918159.5435124.6330（3）在MATLAB工作窗口输入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]=liezhu(A,b)运行后输出结果为请注意：因为RA=RB=n，所以此方程组有唯一解.RA=4RB=4n=4X=0.0000-1.00001.00001.0000h1=10.0000-0.0000-150.0000-762.0000在MATLAB工作窗口输入x=[0;-1;1;1]';x1=[0;-1;1;1]';wucha=x1-x运行后输出结果为wucha=0000（4）解方程组（1）在MATLAB工作窗口输入A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];B=inv(A)运行后结果为B=-268.9293538.3418128.4529106.7599-213.4281-50.956183.3992-166.8022-39.7090在MATLAB工作窗口输入b=[1;1;1];x=inv(A)*b运行后结果为x=397.8654-157.6242-123.1120在MATLAB工作窗口输入A=[3.016.031.999;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.4893在MATLAB工作窗口输入b=[1;1;1];x=inv(A)*b运行后输出结果为x=-4.02641.91931.5210解方程组（2）在MATLAB工作窗口输入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.1969在MATLAB工作窗口输入b=[8;5.900001;5;1];x=inv(A)*b运行后输出结果为x=0-1.00001.00001.0000在MATLAB工作窗口输入A=[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.1969在MATLAB工作窗口输入b=[8;5.900001;5;1];b(2,1)=5.9;x=inv(A)*b运行后输出结果为x=-0.0000-1.00001.00001.0000五、实验结果分析：实验的数学原理很容易理解，也容易上手。把运算的结果带入原方程组，可以发现符合的还是比较好。这说明列主元消去法计算这类方程的有效性。当A可逆时，能够将计算进行到底，列主元法就能确保算法的稳定，而且计算量不大。直接三角消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法。通过以上的计算比较，方程组（1）具有严重的病态性。当系数矩阵有微小的变化时，wucha=-401.8918159.5435124.6330，所得的解与原方程组的解有很大的相对误差。方程组（2）中当系数矩阵A和b有微小变化时，wucha=0000，所得的解与方程组的解没有相对误差。所以方程组（2）是良性的。用MATLAB内部函数inv通过求逆矩阵，然后通过x=inv（A）*b也可以求出方程组的解，但是没有列主元高斯消去法具有良好的稳定性。det函数求方程组系数矩阵的行列式时所得结果和高斯消去法和三角法所得结果相同，具有方便快捷的优点。