数值分析实验报告-清华大学--线性代数方程组的数值解法.docx

资源描述

数值分析实验报告-清华大学--线性代数方程组的数值解法.docx

《数值分析实验报告-清华大学--线性代数方程组的数值解法.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告-清华大学--线性代数方程组的数值解法.docx（15页珍藏版）》请在冰点文库上搜索。

数值分析实验报告-清华大学--线性代数方程组的数值解法.docx

--本页仅作为文档封面，使用时请直接删除即可--

--内页可以根据需求调整合适字体及大小--

数值分析实验报告-清华大学--线性代数方程组的数值解法（总15页）

线性代数方程组的数值解法

实验1.主元的选取与算法的稳定性

问题提出：

Gauss消去法是我们在线性代数中已经熟悉的。

但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？

Gauss消去法从理论算法到数值算法，其关键是主元的选择。

主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验内容：

考虑线性方程组

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

实验要求：

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

取n=10计算矩阵的条件数。

让程序自动选取主元，结果如何？

（2）现选择程序中手动选取主元的功能。

每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。

若每步消去过程总选取按模最大的元素作为主元，结果又如何？

分析实验的结果。

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

（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。

重复上述实验，观察记录并分析实验结果。

程序清单

n=input（'矩阵A的阶数:

n='）;

A=6*diag（ones（1,n））+diag（ones（1,n-1）,1）+8*diag（ones（1,n-1）,-1）;

b=A*ones（n,1）;

p=input（'计算条件数使用p-范数,p='）;

cond_A=cond（A,p）

[m,n]=size（A）;

Ab=[Ab];

r=input（'选主元方式（0:

自动;1:

手动），r='）;

fori=1:

n-1

switchr

case（0）

[aii,ip]=max（abs（Ab（i:

n,i）））;

ip=ip+i-1;

case

（1）

ip=input（['第',num2str（i）,'步消元，请输入第',num2str（i）,'列所选元素所处的行数：

']）;

end;

Ab（[iip],:

）=Ab（[ipi],:

）;

aii=Ab（i,i）;

fork=i+1:

Ab（k,i:

n+1）=Ab（k,i:

n+1）-（Ab（k,i）/aii）*Ab（i,i:

n+1）;

end;

ifr==1

end

end;

x=zeros（n,1）;

x（n）=Ab（n,n+1）/Ab（n,n）;

fori=n-1:

-1:

x（i）=（Ab（i,n+1）-Ab（i,i+1:

n）*x（i+1:

n））/Ab（i,i）;

end

运行结果

（1）n=10，矩阵的条件数及自动选主元

Cond（A,1）=×103

Cond（A,2）=×103

Cond（A,inf）=×103

程序自动选择主元（列主元）

a.输入数据

矩阵A的阶数:

n=10

计算条件数使用p-范数,p=1

选主元方式（0:

自动;1:

手动），r=0

b.计算结果

x=[1，1，1，1，1，1，1，1，1，1]T

（2）n=10，手动选主元

a.每步消去过程总选取按模最小或按模尽可能小的元素作为主元

矩阵A的阶数:

n=10

计算条件数使用p-范数,p=1

选主元方式（0:

自动;1:

手动），r=1

第1步消元，请输入第1列所选元素所处的行数：

第2步消元，请输入第2列所选元素所处的行数：

…（实际选择时，第k步选择主元处于第k行）

最终计算得

x=[,,,,,,,,,]T

b.每步消去过程总选取按模最大的元素作为主元

矩阵A的阶数:

n=10

计算条件数使用p-范数,p=1

选主元方式（0:

自动;1:

手动），r=1

第1步消元，请输入第1列所选元素所处的行数：

第2步消元，请输入第2列所选元素所处的行数：

…（实际选择时，第k步选择主元处于第k+1行）

最终计算得

x=[1，1，1，1，1，1，1，1，1，1]T

（3）n=20，手动选主元

a.每步消去过程总选取按模最小或按模尽可能小的元素作为主元

矩阵A的阶数:

n=20

计算条件数使用p-范数,p=1

选主元方式（0:

自动;1:

手动），r=1

第1步消元，请输入第1列所选元素所处的行数：

第2步消元，请输入第2列所选元素所处的行数：

…（实际选择时，第k步选择主元处于第k行）

最终计算得

x=[，，，，，，，，，，，，，，，，，，，]T

b.每步消去过程总选取按模最大的元素作为主元

矩阵A的阶数:

n=20

计算条件数使用p-范数,p=1

选主元方式（0:

自动;1:

手动），r=1

第1步消元，请输入第1列所选元素所处的行数：

第2步消元，请输入第2列所选元素所处的行数：

…（实际选择时，第k步选择主元处于第k+1行）

最终计算得

x=[1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1]T

（4）A分别为幻方矩阵，Hilbert矩阵，pascal矩阵和随机矩阵

将计算结果列于下表：

11阶幻方矩阵

10阶Hilbert矩阵

选主元方式

模最大

模最小

模最大

模最小

x10

10阶Pascal方矩阵

10阶随机矩阵

选主元方式

模最大

模最小

模最大

模最小

x10

简要分析

计算

（1）表明：

对于同一矩阵，不同范数定义的条件数是不同的；Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。

计算

（2）表明：

通过比较每次选取模最大或模最小的元素的选主元方式，可以发现，在本题给定的问题中，选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确。

因为这样做可以避免使用较小的数做为除数，以免发生结果数量级加大，使大数吃掉小数，产生舍入误差。

计算（3）表明：

首先，n=20得到与计算

（2）,即n=10时一样的结论，即选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确；其次，与计算

（2）（Cond（A10×10,1）=×103）比较，Cond（A20×20,1）=×106显着增大，且计算（3）的误差也远大于计算

（2），即，矩阵的条件数越大，产生的误差也越大。

计算（4）表明：

Gauss消去法在消去过程中，主元的选择与算法的稳定性有密切的联系。

一般来说，选取绝对值大的元素作为主元比绝对值小的元素作为主元时的计算结果更加精确。

但这并不是绝对的，一些特殊的方阵，如Pascal方矩阵，则恰恰是选择模最小的元素作为主元时计算结果最精确（选模最小的元素只是一个表象，这种选主元方法优于其他选主元方法的本质是这种选择方法能使消去过程不产生浮点数，而全是整数运算，只有在回代过程中才有可能会产生浮点数）。

在系数矩阵性质未知，或者说对于绝大多数的系数矩阵来说，选择模最大的元素作为主元是一种比较稳定和精确的方法。

实验2.病态的线性方程组的求解

问题提出：

理论的分析表明，求解病态的线性方程组是困难的。

实际情况是否如此，会出现怎样的现象呢？

实验内容：

考虑方程组Hx=b的求解，其中系数矩阵H为Hilbert矩阵，

这是一个着名的病态问题。

通过首先给定解（例如取为各个分量均为1）再计算出右端b的办法给出确定的问题。

实验要求：

（1）选择问题的维数为6，分别用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何将计算结果与问题的解比较，结论如何

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

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

程序清单

clear

n=6;%Hilbert矩阵的阶数

H=hilb（n）;

cond_H=cond（H）;

b=H*ones（n,1）;

D=diag（diag（H））;

U=-triu（H,1）;

L=-tril（H,-1）;

w=;

tol=;

maxk=10000;

Ab=[Hb];

%Gauss消去法

fori=1:

n-1

[aii,ip]=max（abs（Ab（i:

n,i）））;

ip=ip+i-1;

Ab（[iip],:

）=Ab（[ipi],:

）;

aii=Ab（i,i）;

ifabs（aii）<

disp（'Singularmatrix,stop!

'）;

pause;%如果主元绝对值小于限值，则认为A奇异，终止

end

fork=i+1:

Ab（k,i:

n+1）=Ab（k,i:

n+1）-（Ab（k,i）/aii）*Ab（i,i:

n+1）;

end;

x（n）=Ab（n,n+1）/Ab（n,n）;

fori=n-1:

-1:

x（i）=（Ab（i,n+1）-Ab（i,i+1:

n）*x（i+1:

n））/Ab（i,i）;

end

x_Gauss=x

%Jacobi迭代法

BJ=D\（L+U）;

fJ=D\b;

x（:

1）=*ones（n,1）;

x（:

2）=BJ*x（:

1）+fJ;

k=2;

while（norm（x（:

k）-x（:

k-1））>=tol）&&（k

k=k+1;

x（:

k）=BJ*x（:

k-1）+fJ;

end

k_J=k

x_J=x（:

k）

%GS迭代法

BG=（D-L）\U;

fG=（D-L）\b;

x（:

1）=*ones（n,1）;

x（:

2）=BG*x（:

1）+fG;

k=2;

while（norm（x（:

k）-x（:

k-1））>=tol）&&（k

k=k+1;

x（:

k）=BG*x（:

k-1）+fG;

end;

k_G=k

x_G=x（:

k）

%SOR迭代法

lw=（D-w*L）\（（1-w）*D+w*U）;

fS=（D-w*L）\b*w;

x（:

1）=*ones（n,1）;

x（:

2）=lw*x（:

1）+fS;

k=2;

while（norm（x（:

k）-x（:

k-1））>=tol）&&（k

k=k+1;

x（:

k）=lw*x（:

k-1）+fS;

end

k_SOR=k

x_SOR=x（:

k）

运行结果及简要分析

（1）6阶Hilbert矩阵

计算方法

Gauss

J法

GS法

SOR法

迭代矩阵的1-范数

迭代矩阵的谱半径

迭代次数

489

6665

9951

Inf

NaN

Gauss消元法得到了较为精确的解。

以作为收敛的标准：

J法迭代矩阵的谱半径大于1，迭代不收敛；GS法和SOR法迭代矩阵的谱半径都略小于1，迭代是收敛的，但收敛速度非常慢，在很多次迭代之后仍与精确解有一定误差，GS收敛速度略快于SOR法。

（2）n阶Hilbert矩阵

计算从6阶到25阶的Hilbert矩阵，为观察收敛速度，以作为收敛的标准。

计算结果如下（仅列出有代表性的6、7、10、11、12、13、21、25阶计算结果）

计算方法

Gauss

J法

GS法

SOR法

阶

489

751

607

Inf

NaN

阶

435

593

2537

Inf

NaN

阶

349

1760

1878

Inf

NaN

X10

NaN

阶

332

1606

1589

Inf

NaN

X10

NaN

X11

NaN

阶

318

1447

1367

Inf

NaN

X10

NaN

X11

NaN

X12

NaN

阶

306

1307

1197

Inf

NaN

X10

NaN

X11

NaN

X12

NaN

X13

NaN

阶

252

2586

2893

NaN

X10

NaN

X11

NaN

X12

NaN

X13

NaN

X14

NaN

X15

NaN

X16

NaN

X17

NaN

X18

NaN

X19

NaN

X20

NaN

X21

NaN

阶

237

2253

2240

NaN

X10

NaN

X11

NaN

X12

NaN

X13

NaN

X14

NaN

X15

NaN

X16

NaN

X17

NaN

X18

NaN

X19

NaN

X20

NaN

X21

NaN

X22

NaN

X23

NaN

X24

NaN

X25

NaN

得到的主要结论如下：

（1）Gauss消去法：

Gauss消去法求得的解与精确解的误差随Hilbert矩阵阶数的增加而增加，Hilbert矩阵阶数不大于11时，误差较小（小于1%），对于要求不要高的工程问题，这样的误差可以接受。

当阶数大于11时，误差迅速增加，当阶数为13时，误差已经超过100%，一般来说这样的近似解不可接受。

当阶数为25时，误差达到38900%。

即低阶Hilbert矩阵可用Gauss消元法求解。

（2）Jacobi迭代方法：

无论Hilbert矩阵为多少阶，Jacobi迭代矩阵的谱半径都大于1，迭代不稳定、不收敛。

（3）GS迭代方法：

GS迭代矩阵的谱半径略小于1，迭代收敛，但收敛速度非常缓慢。

（4）SOR迭代方法：

取w=，SOR迭代矩阵的谱半径略小于1，迭代收敛，但收敛速度亦非常缓慢。

在本实验中，多数情况SOR法较GS迭代更慢一些（收敛步数更大）。

从上面的结果可以看出：

Hilbert矩阵阶数较小时，可用Gauss消元法直接求解，解的精度比迭代法高；随着阶数增加，Gauss消元法的精度迅速下降，解变得不可靠，这时，有一些迭代法，如GS法和SOR法仍可继续求得比较精确的解。

另外，在三种迭代法中，GS和SOR法相对Jacobi法更有优势，但这两种方法的迭代矩阵谱半径已经非常接近1（病态问题），收敛速度都很慢。

（3）病态问题的求解

求解病态问题，主要的方法是对原方程进行预处理，以降低系数矩阵的条件数。

例如选择非奇异矩阵，使方程组化为等价方程组，原方程的解。

原则上应使矩阵的条件数比有所改善。

一般和可选择为三角形矩阵或对角矩阵。

理论上最好选择对角阵和，满足：

。

对于1—100阶的Hilbert矩阵，假设，从设为由Hilbert矩阵对角元素组成的对角阵，设是中，使得最小的一种。

的结果如下图所示：

虽然这里取的并不是最优的D矩阵，但在上图中，可观察到经过预处理后的DHD的条件数相比原Hilbert矩阵减小了，在一定程度上改善了原Hilbert矩阵的性质。

1515

展开阅读全文