数值分析第一次上机实验练习报告.docx
《数值分析第一次上机实验练习报告.docx》由会员分享，可在线阅读，更多相关《数值分析第一次上机实验练习报告.docx（20页珍藏版）》请在冰豆网上搜索。
数值分析第一次上机实验练习报告
数值分析第一次上机实验练习报告
——线性方程组求解相关问题
材博141李根2014310937
一、问题的描述
设
是Hilbert矩阵，即
对n=2,3,4…（建议算到15左右）
（a）取
，及
.再用Gauss消去法和Cholesky分解方法来求解
，看看误差有多大.
（b）计算条件数：
cond（Hn）2.
（c）使用某种正则化方法改善（a）中的结果.
二、方法描述与方案设计
我们使用MATLAB来完成此次上机实验。
MATLAB是一款强大的科学计算软件，内置有多种函数可以方便地调用，同时软件界面非常友好，便于学习和使用。
我们的实验方案设计如下：
1、Hilbert矩阵的生成与条件数的计算
MATLAB内置有Hilbert矩阵的生成函数hilb（）和矩阵条件数的计算函数cond（），直接调用即可。
2、Gauss消去法的MATLAB实现
对线性方程组
，将
和
写为增广矩阵
，令之为
。
顺序消元是对增广矩阵进行初等行变换，经n-1步消元将其化为
，其中
位一个上三角矩阵。
过程如下：
经第k-1步消元（
），对应增广矩阵化为
假设
，令
前k-1行不动，将第k行乘以
加到第i行，进行第k步消元，即
这样经n-1步后，得到了矩阵如下
，其中
为一个上三角矩阵，当方程有唯一解时，
可逆，而经初等变换后的系数矩阵行列式不变，即
所以
。
解方程组
，从最后一个方程开始，依次向前一方程代入，即可求解，此过程为回代，计算公式为
通过查询MATLAB应用的相关书籍，我们发现，MATLAB解线性方程组常用的方法有两种，一是在方程Ax=b两边同时乘以A的逆矩阵A-1，二是利用预先定义好的矩阵“左除”：
x=A\b，而后者的算法即为Gauss消元法，因此对于某个给定的Hilbert矩阵
，我们进行矩阵左除运算所得的结果即为Gauss消去法的结果。
3、Cholesky分解法的MATLAB实现
Cholesky分解法又称平方根法，相较于Gauss消去法，Cholesky分解法求解线性方程组时具有较小的计算量，但要求系数矩阵对称且正定。
对系数矩阵
作Cholesky分解即寻找下三角矩阵L，使得
这样的L对角元全为正数且唯一。
分解时，由上式首先可计算出第一列的元素，当我们计算出前j-1列的元素时，第j列的对角元素由下式算出
从而计算出第j列元素
然后我们可依次求解方程组
、
，最终求解，即
MATLAB内置有矩阵Cholesky分解的函数chol（），但需要注意，函数给出的矩阵
为上三角阵，和课本上的定义略有不同。
得到
后，再利用Gauss消去法求解两个新的线性方程组即可。
4、Тихонов正则化方法的MATLAB实现
所谓Тихонов正则化方法其实是一种近似方法，其核心是通过引入了一个小的纯量矩阵来改善原病态矩阵的性质，使其条件数显著降低从而得到更精确的解。
对本实验而言，我们令
，
从而将问题转化为
，该方程仍然由Gauss消去法求解。
需要注意的是，由于引入了纯量矩阵
，现在的解并不是原方程组的精确解，而是一个近似解，因此我们应该合理选区参数
的值使得结果落在一个我们可以接受的范围之内。
此外，原方法中应取矩阵的共轭转置，但由于本实验中Hilbert矩阵为实对称矩阵，因此我们简单地取转置即可。
三、计算结果及分析
1、使用Gauss消去法和Cholesky分解法求解原方程组
我们使用yn1和yn2来表示当系数矩阵为n阶Hilbert矩阵时，利用Gauss消去法和Cholesky分解法得到的原方程组
的解，为定量分析误差，我们引入sigma1和sigma2作为量化标准，它们是求得的解与bn相减得到的残差向量的模。
我们编写的程序如下所示：
>>clearall;
clc;
formatlong;%使用双精度显示
n=input（'n='）;%输入矩阵阶数
Xn=linspace（1,1,n）;%生成行向量之后转置为列向量
Hn=hilb（n）;%生成n阶希尔伯特矩阵
bn=Hn*Xn';%生成原方程右端列向量
Rn=chol（Hn）;%对Hilbert矩阵进行Cholesky分解Hn=Rn'*Rn
zn=Rn'\bn;%Cholesky方法的中间解
yn1=Hn\bn%高斯消元法的解
yn2=Rn\zn%Cholesky方法的解
errn1=（yn1-Xn'）.^2,errn2=（yn2-Xn'）.^2%计算残差平方向量
sigma1=sqrt（sum（errn1））,sigma2=sqrt（sum（errn2））%计算相对Xn的标准偏差
下面的一系列表格给出了两种方法的结果（数据由MATLAB给出，具16位有效数字）：
求解方法
Gauss
Cholesky
所得结果
（n=2）
1.000000000000000
0.999999999999999
1.000000000000000
0.999999999999999
残差向量的模
8.950904182623619e-016
8.950904182623619e-016
求解方法
Gauss
Cholesky
所得结果
（n=3）
0.999999999999999
1.000000000000007
0.999999999999993
0.999999999999999
1.000000000000007
0.999999999999993
残差向量的模
1.017354490604556e-014
1.018383803549894e-014
求解方法
Gauss
Cholesky
所得结果
（n=4）
0.999999999999992
1.000000000000093
0.999999999999767
1.000000000000155
0.999999999999984
1.000000000000183
0.999999999999555
1.000000000000291
残差向量的模
2.948931463388727e-013
5.624264157788779e-013
求解方法
Gauss
Cholesky
所得结果
（n=5）
1.000000000000021
0.999999999999566
1.000000000001974
0.999999999996921
1.000000000001539
0.999999999999944
1.000000000000944
0.999999999996231
1.000000000005375
0.999999999997485
残差向量的模
3.991872242690353e-012
7.092940334727893e-012
求解方法
Gauss
Cholesky
所得结果
（n=6）
0.999999999999127
1.000000000024739
0.999999999833385
1.000000000431972
0.999999999524352
1.000000000187069
0.999999999999122
1.000000000024819
0.999999999833114
1.000000000432145
0.999999999524615
1.000000000186823
残差向量的模
6.900785737295254e-010
6.900074426813605e-010
求解方法
Gauss
Cholesky
所得结果
（n=7）
0.999999999991674
1.000000000335732
0.999999996741135
1.000000012741888
0.999999976535428
1.000000020349350
0.999999993298420
0.999999999992490
1.000000000300196
0.999999997104746
1.000000011264819
0.999999979335571
1.000000017864962
0.999999994131788
残差向量的模
3.439014525977920e-008
3.026516485826501e-008
求解方法
Gauss
Cholesky
所得结果
（n=8）
0.999999999984100
1.000000000892975
0.999999987934077
1.000000067024900
0.999999815794504
1.000000265031910
0.999999808767942
1.000000054585126
0.999999999967524
1.000000001728951
0.999999977534558
1.000000121122128
0.999999674865752
1.000000458968002
0.999999674006787
1.000000091828723
残差向量的模
3.851772775183292e-007
6.680166553264907e-007
求解方法
Gauss
Cholesky
所得结果
（n=9）
0.999999999860744
1.000000009421568
0.999999842880911
1.000001108720852
0.999995970955548
1.000008163929252
0.999990683903648
1.000005596586910
0.999998623655625
0.999999999745327
1.000000017473990
0.999999705372702
1.000002097525213
0.999992322101513
1.000015652254346
0.999982046535176
1.000010833229504
0.999997325588041
残差向量的模
1.428785278708662e-005
2.748214245580457e-005
求解方法
Gauss
Cholesky
所得结果
（n=10）
0.999999998795572
1.000000101747845
0.999997870728245
1.000019082978195
0.999910046468220
1.000244832219840
0.999601711970687
1.000382070330071
0.999800706958791
1.000043578494796
0.999999998774002
1.000000105231645
0.999997769587892
1.000020199312176
0.999903949505314
1.000263370806271
0.999568815579411
1.000415924440761
0.999781991088109
1.000047876491003
残差向量的模
6.439185125919934e-004
6.983810322355236e-004
求解方法
Gauss
Cholesky
所得结果
（n=11）
0.999999997920812
1.000000208162879
0.999994794393711
1.000056428985217
0.999672623147925
1.001124781412886
0.997600118498217
1.003213391149767
0.997373403022044
1.001197774967262
0.999766477566859
0.999999995169251
1.000000508362225
0.999986748262932
1.000148802389948
0.999110102316951
1.003139443198704
0.993143674377092
1.009372318631840
0.992196146947349
1.003618426774329
0.999283830342860
残差向量的模
0.005084195451562
0.014833271077307
求解方法
Gauss
Cholesky
所得结果
（n=12）
0.999999974044749
1.000003334581326
0.999893851846554
1.001461831770340
0.989181965552700
.0479********
0.865468311256050
1.245121567095965
0.710937852345233
1.212819000001615
0.911095819203608
1.016087165180365
0.999999965880967
1.000004316449056
0.999864380709050
1.001846945341169
0.986462994827836
.0594********
0.834278310085044
1.300003329532084
0.648246004872842
1.257647708011406
0.892864938877564
1.019305272304730
残差向量的模
0.466488328530673
0.568465973687995
求解方法
Gauss
Cholesky
所得结果
（n=13）
0.999999985819868
1.000001939754577
0.999932541558656
1.001038048033260
0.991231595703128
.0453********
0.847760592541553
1.342408467694683
0.479726245670882
1.527141311509344
0.658804929786682
1.127592826303361
0.979029832575293
1.000000100800966
0.999983771063206
1.000637205639667
0.989244282315779
.0976********
0.465977282773442
2.876103771702808
-3.375641729215241
7.847802305139874
-6.108210067961484
5.693767548924378
-0.784713030312417
1.297431751763950
残差向量的模
0.907882051103852
12.070004065947840
求解方法
Gauss
Cholesky
所得结果
（n=14）
0.999999482439617
1.000068029614560
0.997857975573164
.027*********
0.833143160876506
1.406101562833689
1.692627958831394
-6.860065364290957
26.482050799360898
-45.462787656255678
53.342818634605649
-35.207104190739898
15.139********5737
-1.391730391291961
残差向量的模
84.421272602931595
求解方法
Gauss
Cholesky
所得结果
（n=15）
0.999999988484888
0.999998820001038
1.000153035493171
0.995678337603479
.0552********
0.606823831028672
2.717791261290424
-3.811419278271013
9.689857910094865
-8.657463514524082
6.344376120567810
1.692954637371676
.0796********
2.820985686699631
0.624725651254457
残差向量的模
15.395524776711373
值得注意的是，当n>13时，MATLAB认为矩阵此时非正定，无法进行Cholesky分解，说明此时矩阵病态十分严重。
我们用残差向量的模来估计计算误差，然后做出误差对数与矩阵阶数的关系图如下：
由图可知，n较小时两种方法的误差基本相同，但当n较大时，Cholesky方法由于受矩阵病态性的影响，误差急剧增大，甚至认为矩阵非正定而无法给出计算结果。
2、计算各阶Hilbert矩阵的条件数
计算一系列Hilbert矩阵条件数的程序如下：
>>formatlong,n=input（'n='）,
fori=1:
n
C（i）=cond（hilb（i））;
end
C
我们得到的结果如下表所示（保留三位有效数字）：
矩阵阶数
条件数
矩阵阶数
条件数
矩阵阶数
条件数
1
1.00
6
1.50E7
11
5.23E14
2
19.3
7
4.75E8
12
1.68E16
3
524
8
1.53E10
13
1.76E18
4
1.55E4
9
4.93E11
14
3.08E17
5
4.77E5
10
1.60E13
15
4.43E17
将矩阵条件数对矩阵阶数作图，如下所示：
由图可知，当n较小时，矩阵条件数的对数随矩阵阶数n的增大而线性增加，而当n较大时，则增长幅度明显下降，并开始出现平台。
3、利用Тихонов正则化方法求解原方程组的研究
我们使用如下程序来实现Тихонов正则化方法：
>>formatlong;
n=input（'n='）;%输入矩阵阶数
a=input（'a='）;%输入正则化参数
Xn=linspace（1,1,n）;%生成行向量之后转置为列向量
Hn=hilb（n）;%生成n阶希尔伯特矩阵
Bn=Hn'*Hn*Xn';%生成新的右端列向量
In=eye（n）;%生成n阶单位阵
An=a.*In+Hn'*Hn;%生成正则化矩阵
yn3=An\Bn%求解正则化方程An*yn3=Bn的解
errn3=（yn3-Xn'）.^2%计算残差平方向量
sigma3=sqrt（sum（errn3））%计算新解相对Xn的标准偏差
在（a）中，我们已经看到，对于n>10的情形，由于Hilbert矩阵病态过于严重，矩阵的条件数过大，我们得到的解甚至无法给出3位有效数字；而现在经过Тихонов正则化处理之后，情况得到了明显的改善，为节约篇幅仅试举n=15为例，此时MATLAB给出的解为：
yn3=
0.999136503929967
1.006390641270175
0.997170658144827
0.991092760424955
0.992635470673801
0.997976335849701
1.003810718709688
1.008287686563388
.010*********
.010*********
1.008128508970063
1.003684743467315
0.997462622419310
0.989750418555657
0.980817121011882
sigma3=
0.032438344028382
显然，此时的结果大大改观，在n=15时仍能给出两位有效数字。
为更清楚地说明问题，我们计算了各正则化An矩阵的条件数，所用程序如下：
>>formatlong;
n=input（'n='）;
a=input（'a='）;
fori=1:
n
Hi=hilb（i）;
Ii=eye（i）;
Ai=a.*Ii+Hi'*Hi;
C（i）=cond（Ai）;
end
C
程序给出的结果如下表所示（仍保留3位有效数字，其中a=10-7）：
矩阵阶数
条件数
矩阵阶数
条件数
矩阵阶数
条件数
1
1.00
6
2.62E7
11
3.15E7
2
372
7
2.76E7
12
3.22E7
3
2.71E5
8
2.88E7
13
3.29E7
4
2.06E7
9
2.98E7
14
3.35E7
5
2.46E7
10
3.07E7
15
3.41E7
我们继续使用图来表示条件数的对数与矩阵阶数之间的关系：
显然，当使用了正则化方法之后，新矩阵明显优于原Hilbert矩阵：
其条件数在n≥4时进入一个稳定的平台，基本保持不变，因此我们在n较大的时候仍能得到稳定的解。
问题看上去已经得到了圆满解决，但其实并未结束：
我们在使用Тихонов正则化方法时，使用了近似，引入了新的参量a，那么a取何值时能够取得最好的效果？
换言之，引入a能够很好地消除舍入误差的影响，但它同时又给我们的计算带来了截断误差，那么a为多少时能够在这两方面得到平衡？
于是，我们继续编写程序进行研究，研究对象是15阶的Hilbert矩阵：
>>formatlong;
n=input（'n='）;
fori=1:
n
a=10^（-i）;
H=hilb（15）;
X=linspace（1,1,15）;
B=H'*H*X';
I=eye（15）;
Ai=a.*I+H'*H;
yi=Ai\B;
erri=（yi-X'）.^2;
S（i）=sqrt（sum（erri））;
end
S
程序给出的结果如下表所示（保留3位有效数字）：
a的负对数
残差向量模
a的负对数
残差向量模
a的负对数
残差向量模
1
1.105
6
6.64E-2
11
3.93E-3
2
0.669
7
3.24E-2
12
2.02E-3
3
0.323
8
1.89E-2
13
6.44E-3
4
0.213
9
1.20E-2
14
7.68E-2
5
9.82E-2
10
5.04E-3
15
0.684
之后作图如下：
由图可见，随着a的减小，误差先减后增，当a=10-12时，残差
升级会员 