解线性方程组的列主元素高斯消去法和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=b LUX=b LY=b A=LU UX=Y；

四、实验内容：

解下列两个线性方程组

?3.016.031.99??x1??1???????4.16?1.23??x2???1? （1） ?1.27?0.987?4.819.34??x??1????3????7?10???32.099999（2） ?5?1??21?1??x1??8??????62??x2??5.900001??? 5?1??x3??5??????????02??x4??1?0a、用你熟悉的算法语言编写程序用列主元高斯消去法和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.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];b=[1;1;1];[RA,RB,n,X]=liezhuY(A,b) D = -0.0305 RA =3 RB =3 n = 3 X =

1.5926 -0.6319 -0.4936 （２）． >>A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];b=[8;5.900001;5;1];[RA,RB,n,X]=liezhuY(A,b) D =-762.0001 RA = 4 RB =4 n =4 X =

0.0000 -1.0000 1.0000 1.0000

ＬＵ分解法程序及结果 （１）

>> A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];h1=zhjLU(A) D = -0.0305 RA =3 U =

3.0100 6.0300 1.9900 0 4.1600 -2.0696

0 0 5.3107 L =

1.0000 0 0 0.4219 1.0000 0 0.3279 -1.6316 1.0000 h1 =

3.0100 4.8635 -0.0305 （２）

>> A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];h1=zhjLU(A) D =-762.0001 RA = 4 U =

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.7333 L =

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

b、将方程组（1）中系数3.01改为3.00，0.987改为0.990，用列主元高斯消去法求解变换后的方程组，输出列主元行交换次序，解向量x及detA，并与（1）中结果比较。 线性方程组数据改变后程序及结果 >>A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 9.34];b=[1;1;1];A(1,1)=3;A(1,3)=0.990;[RA,RB,n,X]=liezhuY(A,b) D = 9.8547 RA = 3 RB =3 n =3 X =

-4.0264 1.9193

1.5210

>> X=[1.5926 -0.6319 -0.4936]';X1=[-4.0264 1.9193 1.5210]';wucha=X1-X wucha =

-5.6190

-4.81

2.5512 2.0146

c、将方程组（2）中的2.099999改为2.1，5.900001改为5.9，用列主元高斯消去法求解变换后的方程组，输出解向量x及detA，并与（1）中的结果比较。 线性方程组数据改变后的程序及结果

>>A=[10 -7 0 1;-3 2.099999 6 2;5 -1 5 2];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.0000 RA = 4 RB =4 n = 4 X =

-0.0000 -1.0000 1.0000 1.0000

>> X=[0;-1;1;1]';X1=[0;-1;1;1]';wucha=X1-X wucha =

0 0 0 0

-1;2

1

0

d、用MATLAB的内部函数inv求出系数矩阵的逆矩阵，再输入命令x=inv(A)*b,即可求出上述各个方程组的解，并与列主元高斯消去法和LU分解法求出的解进行比较，体会选主元的方法具有良好的数值稳定性。用MATLAB的内部函数det求出系数行列式的值，并与（1）、（2）、（3）中输出的系数行列式的值进行比较。 线性方程组求解的程序及结果 （1）

>> A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];B=inv(A) B =

-1.0783 2.1571 0.5138 0.4281 -0.8560 -0.2039 0.3344 -0.6688 -0.1592

>> b=[1;1;1];X=inv(A)*b X =

1.5926 -0.6319 -0.4936

>> A=[3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];A(1,1)=3;A(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 >> b=[1;1;1];X=inv(A)*b X =

-4.0264 1.9193 1.5210 （2）

>> A=[10 -7 0 1;-3 2.099999 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.0108 0.1063 0.0724 -0.0755 0.1024 0.1575 -0.1890 0.1969

>> b=[8;5.900001;5;1];X=inv(A)*b X =

0 -1.0000 1.0000 1.0000

>> A=[10 -7 0 1;-3 2.099999 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.1024 0.1575 -0.1890 0.1969 >> b=[8;5.900001;5;1];b(2,1)=5.9;X=inv(A)*b X =

-0.0000 -1.0000 1.0000

1.0000

六、实验结果分析

通过本次的解线性方程组数值分析实验，不仅让我熟悉了Matlab软件的功能和它强大的计算能力，而且大大提高了我解决数值分析实际问题的能力，进一步掌握了列主元高斯消去法解方程组的思想和有效性。通过对以上方程的计算和比较，方程组（1）具有病态性。当系数矩阵有微小的变化时，所得到的解与原方程组的解有很大的相对误差。方程组（2）中当系数矩阵A和b有较小变化时，所得到的解没有相对误差。所以方程组（2）是良性的。采用Matlab内部函数inv求解逆矩阵，然后通过X=inv（a)*b也可以求解方程组的解，但

没有高斯列主元消去法稳定性高。Det函数求解方程组系数矩阵的行列式所得到的结果和三角法所得到的结果相同，有方便快捷明了的特点。

七、附录

列主元消去法

function [RA,RB,n,X]=liezhuY(A,b)

B=[A b];n=length(b);RA=rank(A); D=det(A)

RB=rank(B);zhica=RB-RA; if zhica>0, return

warning off MATLAB:return_outside_of_loop end if RA==RB if RA==n

X=zeros(n,1);C=zeros(1,1+n); for p=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; 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); end end

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); end else end end

LU分解法

function h1=zhjLU(A) [n n]=size(A);RA=rank(A); D=det(A) if RA~=n return end if RA==n for p=1:n

h(p)=det(A(1:p,1:p)); end

h1=h(1:n); for i=1:n

if h(1,i)==0 return end end

if h(1,i)~=0 for j=1:n

U(1,j)=A(1,j); end for k=2: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 h1;RA,U,L end end

本文来源：https://www.bwwdw.com/article/ju2g.html