清华大学数值分析实验报告.doc

清华大学数值分析实验报告.doc

《清华大学数值分析实验报告.doc》由会员分享，可在线阅读，更多相关《清华大学数值分析实验报告.doc（38页珍藏版）》请在冰豆网上搜索。

数值分析实验报告

一、实验3.1

题目：

考虑线性方程组，，,编制一个能自动选取主元,又能手动选取主元的求解线性代数方程组的Gauss消去过程。

（1）取矩阵，,则方程有解。

取计算矩阵的条件数。

分别用顺序Gauss消元、列主元Gauss消元和完全选主元Gauss消元方法求解，结果如何？

（2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元,结果又如何？

分析实验的结果。

（3）取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用.

（4）选取其他你感兴趣的问题或者随机生成的矩阵,计算其条件数，重复上述实验，观察记录并分析实验的结果。

1.算法介绍

首先,分析各种算法消去过程的计算公式，

顺序高斯消去法:

第k步消去中，设增广矩阵中的元素（若等于零则可以判定系数矩阵为奇异矩阵，停止计算），则对k行以下各行计算,分别用乘以增广矩阵的第行并加到第行，则可将增广矩阵中第列中以下的元素消为零；重复此方法，从第1步进行到第n-1步,则可以得到最终的增广矩阵，即；

列主元高斯消去法：

第k步消去中，在增广矩阵中的子方阵中,选取使得，当时，对中第行与第行交换,然后按照和顺序消去法相同的步骤进行。

重复此方法，从第1步进行第n—1步，就可以得到最终的增广矩阵，即；

完全主元高斯消去法：

第k步消去中，在增广矩阵中对应的子方阵中，选取使得,若或，则对中第行与第行、第列与第列交换，然后按照和顺序消去法相同的步骤进行即可。

重复此方法，从第1步进行到第n—1步,就可以得到最终的增广矩阵,即；

接下来，分析回代过程求解的公式,容易看出，对上述任一种消元法,均有以下计算公式:

2。

实验程序的设计

一、输入实验要求及初始条件；

二、计算系数矩阵A的条件数及方程组的理论解；

三、对各不同方法编程计算,并输出最终计算结果.

3.计算结果及分析

（1）

先计算系数矩阵的条件数，结果如下，

可知系数矩阵的条件数较大,故此问题属于病态问题,b或A的扰动都可能引起解的较大误差；

采用顺序高斯消去法，计算结果为:

最终解为x=（1.000000000000000，1。

000000000000000，1.000000000000000，1。

000000000000001,0.999999999999998，1。

000000000000004,0。

999999999999993,1.000000000000012，0。

999999999999979，1。

000000000000028）T

使用无穷范数衡量误差,得到=2.842170943040401e—14，可以发现，采用顺序高斯消元法求得的解与精确解之间误差较小。

通过进一步观察，可以发现,按照顺序高斯消去法计算时，其选取的主元值和矩阵中其他元素大小相近，因此顺序高斯消去法方式并没有对结果造成特别大的影响。

若采用列主元高斯消元法，则结果为：

最终解为x=（1。

000000000000000,1。

000000000000000，1。

000000000000000，1.000000000000000,1。

000000000000000，1。

000000000000000,1。

000000000000000,1。

000000000000000，1.000000000000000,1.000000000000000）T

同样使用无穷范数衡量误差，有=0;

若使用完全主元高斯消元法,则结果为

最终解x=（1.000000000000000,1。

000000000000000，1。

000000000000000，1。

000000000000000，1.000000000000000,1。

000000000000000，1。

000000000000000,1。

000000000000000，1。

000000000000000，1。

000000000000000）T

同样使用无穷范数衡量误差，有=0；

（2）

若每步都选取模最小或尽可能小的元素为主元，则计算结果为

最终解x=（1。

0000000000000001.0000000000000001。

0000000000000001。

0000000000000010。

9999999999999981.0000000000000040。

9999999999999931。

0000000000000120.9999999999999791.000000000000028）T

使用无穷范数衡量误差,有为2。

842170943040401e-14；而完全主元消去法的误差为=0。

从

（1）和

（2）的实验结果可以发现，列主元消去法和完全主元消去法都得到了精确解，而顺序高斯消去法和以模尽量小的元素为主元的消去法没有得到精确解.在后两种消去法中，由于程序计算时的舍入误差，对最终结果产生了一定的影响，但由于方程组的维度较低，并且元素之间相差不大,所以误差仍比较小。

为进一步分析，计算上述4种方法每步选取的主元数值，并列表进行比较，结果如下：

第n次消元

顺序

列主元

完全主元

模最小

1

6.000000000000000

8

8

6.000000000000000

2

4.666666666666667

8

8

4.666666666666667

3

4.285714285714286

8

8

4。

285714285714286

4

4。

133333333333333

8

8

4。

133333333333333

5

4。

064516129032258

8

8

4.064516129032258

6

4。

031746031746032

8

8

4.031746031746032

7

4。

015748031496063

8

8

4.015748031496063

8

4。

007843137254902

8

8

4。

007843137254902

9

4。

003913894324853

8

8

4。

003913894324853

10

4。

001955034213099

0.015617370605469

0.015617370605469

4。

001955034213099

从上表可以发现，对这个方程组而言，顺序高斯消去选取的主元恰好事模尽量小的元素，而由于列主元和完全主元选取的元素为8，与4在数量级上差别小,所以计算过程中的累积误差也较小,最终4种方法的输出结果均较为精确.

在这里，具体解释一下顺序法与模最小法的计算结果完全一致的原因。

该矩阵在消元过程中，每次选取主元的一列只有两个非零元素,对角线上的元素为4左右，而其正下方的元素为8，该列其余位置的元素均为0。

在这样的情况下，默认的主元也就是该列最小的主元，因此两种方法所得到的计算结果是一致的。

理论上说，完全高斯消去法的误差最小，其次是列主元高斯消去法，而选取模最小的元素作为主元时的误差最大，但是由于方程组的特殊性（元素相差不大并且维度不高），这个理论现象在这里并没有充分体现出来。

（3）

时，重复上述实验过程，各种方法的计算结果如下所示,在这里，仍采用无穷范数衡量绝对误差。

顺序高斯消去法

列主元高斯消去

完全主元高斯消去

选取模最小或尽可能小元素作为主元消去

X

1。

0000000000000001.0000000000000001.0000000000000001.0000000000000010。

9999999999999981.0000000000000040。

9999999999999931.0000000000000140.9999999999999721.0000000000000570.999999999999886

1。

0000000000002270.9999999999995471.0000000000009020。

9999999999982091。

0000000000035240.9999999999931791。

0000000000127320.9999999999781731。

000000000029102

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1.0000000000000001。

0000000000000001。

0000000000000001.0000000000000010。

9999999999999981.0000000000000040。

9999999999999931。

0000000000000140.9999999999999721。

0000000000000570.999999999999886

1。

0000000000002270.9999999999995471。

0000000000009020.9999999999982091.0000000000035240。

9999999999931791.0000000000127320.9999999999781731.000000000029102

2.910205409989430e—11

0

0

2.910205409989430e-11

可以看出,此时列主元和完全主元的计算结果仍为精确值,而顺序高斯消去和模尽可能小方法仍然产生了一定的误差，并且两者的误差一致。

与n=10时候的误差比相比,n=20时的误差增长了大约1000倍，这是由于计算过程中舍入误差的不断累积所致。

所以,如果进一步增加矩阵的维数，应该可以看出更明显的现象。

（4）

不同矩阵维度下的误差如下,在这里,为方便起见，选取2—条件数对不同维度的系数矩阵进行比较。

维度

条件数

顺序消去

列主元

完全主元

模尽量小

1.7e+3

2.84e-14

0

0

2。

84e—14

1。

8e+6

2。

91e-11

0

0

2。

91e—11

5。

7e+7

9.31e—10

0

0

9。

31e—10

1.8e+9

2。

98e-08

0

0

2.98e-08

1。

9e+12

3.05e—05

0

0

3.05e—05

3。

8e+16

3.28e+04

3。

88e-12

3.88e—12

3.28e+04

8.5e+16

3.52e+13

4。

2e—3

4。

2e-3

3。

52e+13

从上表可以看出，随着维度的增加,不同方法对计算误差的影响逐渐体现,并且增长较快，这是由于舍入误差逐步累计而造成的。

不过，方法二与方法三在维度小于40的情况下都得到了精确解，这两种方法的累计误差远比方法一和方法四慢；同样地，出于与前面相同的原因，方法一与方法四的计算结果保持一致，方法二与方法三的计算结果保持一致.

4.结论

本文矩阵中的元素差别不大，模最大和模最小的元素并没有数量级上的差异，因此,不同的主元选取方式对计算结果的影响在维度较低的情况下并不明显，四种方法都足够精确.

对比四种方法，可以发现采用列主元高斯消去或者完全主元高斯消去法，可以尽量抑制误差，算法最为精确.不过，对于低阶的矩阵来说，四种方法求解出来的结果误差均较小。

另外，由于完全选主元方法在选主元的过程中计算量较大，而且可以发现列主元法已经可以达到很高的精确程度，因而在实际计算中可以选用列主元法进行计算。

附录:

程序代码

clear

clc；

formatlong;

%方法选择

n=input（’矩阵A阶数:

n='）；

disp（'选取求解方式’）;

disp（’1顺序Gauss消元法,2列主元Gauss消元法，3完全选主元Gauss消元法,4模最小或近可能小的元素作为主元'）；

a=input（'求解方式序号：

'）;

%赋值A和b

A=zeros（n,n）;

b=zeros（n,1）；

fori=1：

n

A（i，i）=6;

ifi〉1

A（i,i-1）=8；

end

ifi〈n

A（i，i+1）=1;

end

end

fori=1：

n

forj=1：

n

b（i）=b（i）+A（i，j）；

end

end

disp（'给定系数矩阵为：

’）;

A

disp（'右端向量为：

’）;

b

％求条件数及理论解

disp（’线性方程组的精确解：

’）;

X=（A\b）’

fprintf（’矩阵A的1-条件数:

％f\n',cond（A,1））；

fprintf（’矩阵A的2—条件数:

%f\n’，cond（A））;

fprintf（'矩阵A的无穷-条件数：

％f\n’，cond（A，inf））;

％顺序Gauss消元法

ifa==1

A1=A；b1=b；

fork=1:

n

ifA1（k，k）==0

disp（’主元为零，顺序Gauss消元法无法进行’）；

break

end

fprintf（’第%d次消元所选取的主元：

％g\n',k，A1（k，k））

%disp（’此次消元后系数矩阵为：

’）；

％A1

forp=k+1：

n

l=A1（p，k）/A1（k，k）;

A1（p，k：

n）=A1（p，k:

n）—l*A1（k,k：

n）;

b1（p）=b1（p）—l*b1（k）；

end

end

x1（n）=b1（n）/A1（n，n）;

fork=n-1：

—1:

1

forw=k+1：

n

b1（k）=b1（k）-A1（k，w）*x1（w）;

end

x1（k）=b1（k）/A1（k,k）;

end

disp（'顺序Gauss消元法解为：

’）;

disp（x1）;

disp（'所求解与精确解之差的无穷—范数为'）;

norm（x1—X,inf）

end

％列主元Gauss消元法

ifa==2

A2=A;b2=b；

fork=1：

n

[max_i,max_j]=find（A2（：

，k）==max（abs（A2（k:

n，k））））;

ifmax_i~=k

A2_change=A2（k，：

）；

A2（k，：

）=A2（max_i,:

）；

A2（max_i，:

）=A2_change;

b2_change=b2（k）;

b2（k）=b2（max_i）；

b2（max_i）=b2_change;

end

ifA2（k，k）==0

disp（'主元为零,列主元Gauss消元法无法进行'）;

break

end

fprintf（'第％d次消元所选取的主元：

%g\n’,k,A2（k,k））

％disp（’此次消元后系数矩阵为：

'）;

%A2

forp=k+1：

n

l=A2（p，k）/A2（k，k）；

A2（p，k：

n）=A2（p，k:

n）—l＊A2（k,k：

n）；

b2（p）=b2（p）—l*b2（k）;

end

end

x2（n）=b2（n）/A2（n,n）；

fork=n—1:

-1:

1

forw=k+1：

n

b2（k）=b2（k）—A2（k，w）*x2（w）；

end

x2（k）=b2（k）/A2（k，k）；

end

disp（'列主元Gauss消元法解为:

’）；

disp（x2）；

disp（'所求解与精确解之差的无穷—范数为’）；

norm（x2-X,inf）

end

%完全选主元Gauss消元法

ifa==3

A3=A;b3=b;

fork=1：

n

VV=eye（n）；

［max_i,max_j］=find（A3（k:

n,k：

n）==max（max（abs（A3（k:

n,k：

n）））））；

ifnumel（max_i）==0

[max_i，max_j］=find（A3（k:

n，k：

n）==-max（max（abs（A3（k：

n，k：

n）））））；

end

W=eye（n）；

W（max_i

（1）+k—1，max_i

（1）+k-1）=0；

W（k,k）=0;

W（max_i

（1）+k—1,k）=1;

W（k,max_i

（1）+k—1）=1;

V=eye（n）；

V（k,k）=0;

V（max_j

（1）+k—1，max_j

（1）+k-1）=0；

V（k，max_j

（1）+k—1）=1;

V（max_j

（1）+k-1,k）=1;

A3=W＊A3＊V;

b3=W＊b3;

VV=VV*V;

ifA3（k,k）==0

disp（'主元为零，完全选主元Gauss消元法无法进行'）；

break

end

fprintf（'第％d次消元所选取的主元:

%g\n'，k,A3（k,k））

％disp（’此次消元后系数矩阵为：

'）；

％A3

forp=k+1:

n

l=A3（p，k）/A3（k，k）;

A3（p，k:

n）=A3（p,k：

n）-l＊A3（k，k:

n）；

b3（p）=b3（p）—l＊b3（k）;

end

end

x3（n）=b3（n）/A3（n,n）；

fork=n-1:

—1:

1

forw=k+1：

n

b3（k）=b3（k）-A3（k,w）*x3（w）；

end

x3（k）=b3（k）/A3（k，k）;

end

disp（'完全选主元Gauss消元法解为：

’）；

disp（x3）;

disp（'所求解与精确解之差的无穷—范数为'）；

norm（x3—X，inf）

end

%模最小或近可能小的元素作为主元

ifa==4

A4=A;b4=b；

fork=1：

n

AA=A4；

AA（AA==0）=NaN；

[min_i，j]=find（AA（k:

n,k）==min（abs（AA（k:

n,k））））;

ifnumel（min_i）==0

[min_i，j]=find（AA（k:

n，k）==-min（abs（AA（k：

n，k:

n））））;

end

W=eye（n）；

W（min_i

（1）+k—1,min_i

（1）+k-1）=0；

W（k，k）=0；

W（min_i

（1）+k-1，k）=1;

W（k,min_i

（1）+k-1）=1;

A4=W*A4;

b4=W*b4;

ifA4（k,k）==0

disp（'主元为零，模最小Gauss消元法无法进行'）;

break

end

fprintf（’第%d次消元所选取的主元：

%g\n'，k,A4（k,k））

％A4

forp=k+1：

n

l=A4（p，k）/A4（k，k）;

A4（p,k：

n）=A4（p,k:

n）—l*A4（k,k:

n）；

b4（p）=b4（p）—l＊b4（k）；

end

end

x4（n）=b4（n）/A4（n,n）;

fork=n—1:

-1:

1

forw=k+1:

n

b4（k）=b4（k）-A4（k,w）*x4（w）；

end

x4（k）=b4（k）/A4（k，k）;

end

disp（'模最小Gauss消元法解为:

'）;

disp（x4）；

disp（'所求解与精确解之差的无穷-范数为’）;

norm（x4-X,inf）

end

二、实验3.3

题目：

考虑方程组的解,其中系数矩阵H为Hilbert矩阵：

这是一个著名的病态问题.通过首先给定解（例如取为各个分量均为1）再计算出右端的办法给出确定的问题。

（1）选择问题的维数为6,分别用Gauss消去法（即LU分解）、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？

将计算结果与问题的解比较,结论如何。

（2）逐步增大问题的维数，仍用上述的方法来解它们，计算的结果如何？

计算的结果说明的什么?

（3）讨论病态问题求解的算法。

1.算法设计

对任意线性方程组，分析各种方法的计算公式如下，

（1）Gauss消去法：

首先对系数矩阵进行LU分解，有,则原方程转化为，令，则原方程可以分为两步回代求解：

具体方法这里不再赘述.

（2）J迭代法:

首先分解，再构造迭代矩阵,其中，进行迭代计算，直到误差满足要求。

（3）GS迭代法：

首先分解，再构造迭代矩阵，其中，进行迭代计算，直到误差满足要求。

（4）SOR迭代法：

首先分解，再构造迭代矩阵，其中，进行迭代计算，直到误差满足要求.

2。

实验过程

一、根据维度n确定矩阵H的