数值分析第二章实验.docx

数值分析第二章实验.docx

《数值分析第二章实验.docx》由会员分享，可在线阅读，更多相关《数值分析第二章实验.docx（11页珍藏版）》请在冰豆网上搜索。

数值分析第二章实验

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

一、实验目的：

通过数值实验，从中体会解线性方程组选主元的必要性和LU分解法的优点，以及方程组系数矩阵和右端向量的微小变化对解向量的影响。

二、实验内容：

解下列两个线性方程组

三、实验要求：

（1）用你熟悉的算法语言编写程序用列主元高斯消去法和LU分解求解上述两个方程组，输出Ax=b中矩阵A及向量b,A=LU分解的L及U，detA及解向量x.

（2）将方程组

中系数3.01改为3.00，0.987改为0.990，用列主元高斯消去法求解变换后的方程组，输出列主元行交换次序，解向量x及detA，并与

中结果比较。

（3）将方程组

中的2.099999改为2.1，5.900001改为5.9，用列主元高斯消去法求解变换后的方程组，输出解向量x及detA，并与

中的结果比较。

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

用MATLAB的内部函数det求出系数行列式的值，并与

（1）、

（2）、（3）中输出的系数行列式的值进行比较。

四、实验过程：

（1）列主元高斯消去法的主程序为

function[RA,RB,n,X]=liezhu（A,b）

B=[Ab];n=length（b）;RA=rank（A）;

RB=rank（B）;zhica=RB-RA;

ifzhica>0,

disp（'请注意：

因为RA~=RB，所以此方程组无解.'）

return

end

ifRA==RB

ifRA==n

disp（'请注意：

因为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:

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）;

forq=n-1:

-1:

1

X（q）=（b（q）-sum（A（q,q+1:

n）*X（q+1:

n）））/A（q,q）;

end

else

disp（'请注意：

因为RA=RB

end

end

解方程组

（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]=liezhu（A,b）

运行后输出结果为

请注意：

因为RA=RB=n，所以此方程组有唯一解.

RA=3RB=3n=3

X=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，所以此方程组有唯一解.

RA=4RB=4n=4

X=0.0000

-1.0000

1.0000

1.0000

LU分解法及MATLAB主程序为

functionhl=zhjLU（A）

[nn]=size（A）;RA=rank（A）;

ifRA~=n

disp（'请注意：

因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下:

'）,RA,hl=det（A）;

return

end

ifRA==n

forp=1:

n

h（p）=det（A（1:

p,1:

p））;

end

hl=h（1:

n）;

fori=1:

n

ifh（1,i）==0

disp（'请注意：

因为A的r阶主子式等于零，所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：

'）,hl;RA

return

end

end

ifh（1,i）~=0

disp（'请注意：

因为A的各阶主子式都不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下：

'）

forj=1:

n

U（1,j）=A（1,j）;

end

fork=2:

n

fori=2:

n

forj=2:

n

L（1,1）=1;L（i,i）=1;

ifi>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;RA,U,L

end

end

解方程组

（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依次如下：

RA=3

U=3.01006.03001.9990

04.1600-2.0734

005.3016

L=1.000000

0.42191.00000

0.3279-1.63161.0000

h1=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依次如下：

RA=4

U=10.0000-7.000001.0000

02.10006.00002.3000

00-2.1429-4.2381

0-0.0000012.7333

L=1.0000000

-0.30001.000000

0.50001.19051.0000-0.0000

0.20001.14293.20001.0000

h1=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=3

X=-4.0264

1.9193

1.5210

hi=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=4

X=0.0000

-1.0000

1.0000

1.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=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.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.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.0264

1.9193

1.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.2690

0.01080.10630.0724-0.0755

0.10240.1575-0.18900.1969

在MATLAB工作窗口输入

b=[8;5.900001;5;1];x=inv（A）*b

运行后输出结果为

x=0

-1.0000

1.0000

1.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.2690

0.01080.10630.0724-0.0755

0.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.0000

1.0000

1.0000

用det函数求方程组

（1）系数行列式

在MATLAB工作窗口输入程序为

A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];det（A）

运行输出结果为

ans=-0.1225

在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;det（A）

运行输出结果为

ans=9.8547

用det函数求方程组

（2）的系数行列式

在MATLAB工作窗口输入程序为

A=[10-701;-32.09999962;5-15-1;2102];det（A）

运行结果输出为

ans=-762.0001

在MATLAB工作窗口输入程序为

A=[10-701;-32.09999962;5-15-1;2102];A（2,2）=2.1;det（A）

运行结果输出为

ans=-762.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函数求方程组系数矩阵的行列式时所得结果和高斯消去法和三角法所得结果相同，具有方便快捷的优点。