数值分析第一次上机实验练习报告Word格式.docx
- 文档编号:16204009
- 上传时间:2022-11-21
- 格式:DOCX
- 页数:20
- 大小:147.23KB
数值分析第一次上机实验练习报告Word格式.docx
《数值分析第一次上机实验练习报告Word格式.docx》由会员分享,可在线阅读,更多相关《数值分析第一次上机实验练习报告Word格式.docx(20页珍藏版)》请在冰豆网上搜索。
作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
残差向量的模
8.950904182623619e-016
(n=3)
1.000000000000007
0.999999999999993
1.017354490604556e-014
1.018383803549894e-014
(n=4)
0.999999999999992
1.000000000000093
0.999999999999767
1.000000000000155
0.999999999999984
1.000000000000183
0.999999999999555
1.000000000000291
2.948931463388727e-013
5.624264157788779e-013
(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
(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
(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
(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
(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
(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
(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
(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
(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
(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
(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('
),
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、利用Тихонов正则化方法求解原方程组的研究
我们使用如下程序来实现Тихонов正则化方法:
a=input('
a='
%输入正则化参数
Bn=Hn'
*Hn*Xn'
%生成新的右端列向量
In=eye(n);
%生成n阶单位阵
An=a.*In+Hn'
*Hn;
%生成正则化矩阵
yn3=An\Bn%求解正则化方程An*yn3=Bn的解
errn3=(yn3-Xn'
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*********
1.008128508970063
1.003684743467315
0.997462622419310
0.989750418555657
0.980817121011882
sigma3=
0.032438344028382
显然,此时的结果大大改观,在n=15时仍能给出两位有效数字。
为更清楚地说明问题,我们计算了各正则化An矩阵的条件数,所用程序如下:
Hi=hilb(i);
Ii=eye(i);
Ai=a.*Ii+Hi'
*Hi;
C(i)=cond(Ai);
程序给出的结果如下表所示(仍保留3位有效数字,其中a=10-7):
2.62E7
3.15E7
372
2.76E7
3.22E7
2.71E5
2.88E7
3.29E7
2.06E7
2.98E7
3.35E7
2.46E7
3.07E7
3.41E7
我们继续使用图来表示条件数的对数与矩阵阶数之间的关系:
显然,当使用了正则化方法之后,新矩阵明显优于原Hilbert矩阵:
其条件数在n≥4时进入一个稳定的平台,基本保持不变,因此我们在n较大的时候仍能得到稳定的解。
问题看上去已经得到了圆满解决,但其实并未结束:
我们在使用Тихонов正则化方法时,使用了近似,引入了新的参量a,那么a取何值时能够取得最好的效果?
换言之,引入a能够很好地消除舍入误差的影响,但它同时又给我们的计算带来了截断误差,那么a为多少时能够在这两方面得到平衡?
于是,我们继续编写程序进行研究,研究对象是15阶的Hilbert矩阵:
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));
S
程序给出的结果如下表所示(保留3位有效数字):
a的负对数
残差向量模
1.105
6.64E-2
3.93E-3
0.669
3.24E-2
2.02E-3
0.323
1.89E-2
6.44E-3
0.213
1.20E-2
7.68E-2
9.82E-2
5.04E-3
0.684
之后作图如下:
由图可见,随着a的减小,误差先减后增,当a=10-12时,残差
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 第一次 上机 实验 练习 报告