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

1、解线性方程组的列主元素高斯消去法和LU分解法一、实验目的：通过数值实验，从中体会解线性方程组选主元的必要性和LU分解法的优点，以及方程组系数矩阵和右端向量的微小变化对解向量的影响。二、实验内容：解下列两个线性方程组（1）（2）三、实验要求：（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、2.改为2.1，5.改为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=A b; n=length(b); RA=rank(A)

3、; RB=rank(B);zhica=RB-RA;D=det(A)if zhica0,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)

4、;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 dispendend解方程组（1）在MATLAB工作窗口输入A=3.01 6.03 1.999;1.27 4.16 -1.23;0.987 -4.81 9.34;b=1;1;1;RA,RB,n,X=liezhuY(A,b)运行后输出结果为请注意：因为RA=RB=n，所以此方程组有唯一解.D=-0.1225RA =3 RB =3 n =3X = 397.8654

5、 -157.6242 -123.1120解方程组（2）在MATLAB工作窗口输入 A=10 -7 0 1;-3 2. 6 2;5 -1 5 -1;2 1 0 2;b=8;5.;5;1;RA,RB,n,X=liezhu(A,b)运行后输出结果为请注意：因为RA=RB=n，所以此方程组有唯一解.D=-762.0000RA =4 RB =4 n =4X =0.0000 -1.0000 1.0000 1.0000LU分解法及MATLAB主程序为function hl=zhjLU(A)n n =size(A); RA=rank(A);D=det（A）if RA=ndisp(请注意：因为A的n阶行列式hl

6、等于零，所以A不能进行LU分解.A的秩RA如下:), RA,hl=det(A);returnendif RA=n for p=1:nh(p)=det(A(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp(请注意：因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：), hl;RAreturnendend if h(1,i)=0 disp(请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：)for j=1:nU(1,j)=A(1,j);endfor k=2:n

7、for i=2:n for j=2:n L(1,1)=1;L(i,i)=1; if ijL(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.01 6.03 1.999;1.27 4.16 -1.23;0.987 -4.81 9.34;h1=zhjLU(A)运行输出结果为请注意

8、：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：D=9.8547RA =3 U =3.0100 6.0300 1.9990 0 4.1600 -2.0734 0 0 5.3016L =1.0000 0 0 0.4219 1.0000 0 0.3279 -1.6316 1.0000h1 =3.0100 4.8635 -0.1225解方程组（2）在MATLAB工作窗口输入A=10 -7 0 1;-3 2. 6 2;5 -1 5 -1;2 1 0 2;h1=zhjLU(A)运行后输出结果为请注意：因为A的各阶主子式都不等于零，所以A能进行LU分解.A的

9、秩RA和各阶顺序主子式值hl依次如下：D=-762.0000RA =4U =10.0000 -7.0000 0 1.0000 0 2.1000 6.0000 2.3000 0 0 -2.1429 -4.2381 0 -0.0000 0 12.7333L =1.0000 0 0 0 -0.3000 1.0000 0 0 0.5000 1.1905 1.0000 -0.0000 0.2000 1.1429 3.2000 1.0000 h1 =10.0000 -0.0000 -150.0001 -762.0001（2）在MATLAB工作窗口输入A=3.01 6.03 1.999;1.27 4.16

10、-1.23;0.987 -4.81 9.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 =3 RB =3 n =3 X = -4.0264 1.9193 1.5210hi = 3.0000 4.8219 9.8547在MATLAB工作窗口输入x=397.8654;-157.6242;-123.1120;x1=-4.0264;1.9193;1.5210;wucha=x1-x运行后输出结果为wucha =-401.8918 159.5435 124.6330（3）在MATLAB工作窗

11、口输入A=10 -7 0 1;-3 2. 6 2;5 -1 5 -1;2 1 0 2;A(2,2)=2.1;b(2,1)=5.9;b=8;5.;5;1;RA,RB,n,X=liezhu(A,b)运行后输出结果为请注意：因为RA=RB=n，所以此方程组有唯一解.RA =4 RB =4 n =4X =0.0000 -1.0000 1.00001.0000 h1 =10.0000 -0.0000 -150.0000 -762.0000在MATLAB工作窗口输入x=0;-1;1;1;x1=0;-1;1;1;wucha=x1-x运行后输出结果为wucha = 0 0 0 0（4）解方程组（1）在MATL

12、AB工作窗口输入A=3.01 6.03 1.999;1.27 4.16 -1.23;0.987 -4.81 9.34;B=inv(A)运行后结果为B =-268.9293 538.3418 128.4529 106.7599 -213.4281 -50.9561 83.3992 -166.8022 -39.7090在MATLAB工作窗口输入b=1;1;1;x=inv(A)*b运行后结果为x =397.8654 -157.6242 -123.1120在MATLAB工作窗口输入A=3.01 6.03 1.999;1.27 4.16 -1.23;0.987 -4.81 9.34;A(1,1)=3;A

13、(1,3)=0.990;B=inv(A)运行输出结果为B = 3.3424 -6.1983 -1.1705 -1.3269 2.7442 0.5020 -1.0365 2.0682 0.4893在MATLAB工作窗口输入b=1;1;1;x=inv(A)*b运行后输出结果为x =-4.0264 1.9193 1.5210解方程组（2）在MATLAB工作窗口输入A=10 -7 0 1;-3 2. 6 2;5 -1 5 -1;2 1 0 2;B=inv(A)运行后结果为B =-0.0223 -0.0984 0.1181 0.1686 -0.1601 -0.1181 0.1417 0.2690 0.0

14、108 0.1063 0.0724 -0.0755 0.1024 0.1575 -0.1890 0.1969在MATLAB工作窗口输入b=8;5.;5;1;x=inv(A)*b运行后输出结果为x = 0 -1.0000 1.0000 1.0000在MATLAB工作窗口输入A=10 -7 0 1;-3 2. 6 2;5 -1 5 -1;2 1 0 2;A(2,2)=2.1;B=inv(A)运行后输出结果为B =-0.0223 -0.0984 0.1181 0.1686 -0.1601 -0.1181 0.1417 0.2690 0.0108 0.1063 0.0724 -0.0755 0.102

15、4 0.1575 -0.1890 0.1969在MATLAB工作窗口输入b=8;5.;5;1;b(2,1)=5.9;x=inv(A)*b运行后输出结果为x =-0.0000 -1.0000 1.00001.0000五、实验结果分析： 实验的数学原理很容易理解，也容易上手。把运算的结果带入原方程组，可以发现符合的还是比较好。这说明列主元消去法计算这类方程的有效性。当A可逆时，能够将计算进行到底，列主元法就能确保算法的稳定，而且计算量不大。直接三角消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法。通过以上的计算比较，方程组（1）具有严重的病态性。当系数矩阵有微小的变