数值分析第一次上机实验练习报告.docx

1、数值分析第一次上机实验练习报告数值分析第一次上机实验练习报告线性方程组求解相关问题材博141 李根 2014310937一、 问题的描述设是Hilbert矩阵，即对n=2,3,4（建议算到15左右）(a) 取，及. 再用Gauss消去法和Cholesky分解方法来求解，看看误差有多大.(b) 计算条件数：cond(Hn)2.(c) 使用某种正则化方法改善(a)中的结果.二、 方法描述与方案设计我们使用MATLAB来完成此次上机实验。MATLAB是一款强大的科学计算软件，内置有多种函数可以方便地调用，同时软件界面非常友好，便于学习和使用。我们的实验方案设计如下：1、 Hilbert矩阵的生成与条

2、件数的计算MATLAB内置有Hilbert矩阵的生成函数hilb()和矩阵条件数的计算函数cond()，直接调用即可。2、 Gauss消去法的MATLAB实现对线性方程组，将和写为增广矩阵，令之为。顺序消元是对增广矩阵进行初等行变换，经n-1步消元将其化为，其中位一个上三角矩阵。过程如下：经第k-1步消元（），对应增广矩阵化为假设，令前k-1行不动，将第k行乘以加到第i行，进行第k步消元，即这样经n-1步后，得到了矩阵如下，其中为一个上三角矩阵，当方程有唯一解时，可逆，而经初等变换后的系数矩阵行列式不变，即所以。解方程组，从最后一个方程开始，依次向前一方程代入，即可求解，此过程为回代，计算公式

3、为通过查询MATLAB应用的相关书籍，我们发现，MATLAB解线性方程组常用的方法有两种，一是在方程Ax=b两边同时乘以A的逆矩阵A-1，二是利用预先定义好的矩阵“左除”：x=Ab，而后者的算法即为Gauss消元法，因此对于某个给定的Hilbert矩阵，我们进行矩阵左除运算所得的结果即为Gauss消去法的结果。3、 Cholesky分解法的MATLAB实现Cholesky分解法又称平方根法，相较于Gauss消去法，Cholesky分解法求解线性方程组时具有较小的计算量，但要求系数矩阵对称且正定。对系数矩阵作Cholesky分解即寻找下三角矩阵L，使得这样的L对角元全为正数且唯一。分解时，由上式

4、首先可计算出第一列的元素，当我们计算出前j-1列的元素时，第j列的对角元素由下式算出从而计算出第j列元素然后我们可依次求解方程组、，最终求解，即MATLAB内置有矩阵Cholesky分解的函数chol()，但需要注意，函数给出的矩阵为上三角阵，和课本上的定义略有不同。得到后，再利用Gauss消去法求解两个新的线性方程组即可。4、 正则化方法的MATLAB实现所谓正则化方法其实是一种近似方法，其核心是通过引入了一个小的纯量矩阵来改善原病态矩阵的性质，使其条件数显著降低从而得到更精确的解。对本实验而言，我们令，从而将问题转化为，该方程仍然由Gauss消去法求解。需要注意的是，由于引入了纯量矩阵，现

5、在的解并不是原方程组的精确解，而是一个近似解，因此我们应该合理选区参数的值使得结果落在一个我们可以接受的范围之内。此外，原方法中应取矩阵的共轭转置，但由于本实验中Hilbert矩阵为实对称矩阵，因此我们简单地取转置即可。三、 计算结果及分析1、 使用Gauss消去法和Cholesky分解法求解原方程组我们使用yn1和yn2来表示当系数矩阵为n阶Hilbert矩阵时，利用Gauss消去法和Cholesky分解法得到的原方程组的解，为定量分析误差，我们引入sigma1和sigma2作为量化标准，它们是求得的解与bn相减得到的残差向量的模。我们编写的程序如下所示：clear all;clc;form

6、at long; %使用双精度显示n=input(n=); %输入矩阵阶数Xn=linspace(1,1,n); %生成行向量之后转置为列向量Hn=hilb(n); %生成n阶希尔伯特矩阵bn=Hn*Xn; %生成原方程右端列向量Rn=chol(Hn); %对Hilbert矩阵进行Cholesky分解Hn=Rn*Rnzn=Rnbn; %Cholesky方法的中间解yn1=Hnbn %高斯消元法的解yn2=Rnzn %Cholesky方法的解errn1=(yn1-Xn).2,errn2=(yn2-Xn).2 %计算残差平方向量sigma1=sqrt(sum(errn1),sigma2=sqrt(

7、sum(errn2) %计算相对Xn的标准偏差下面的一系列表格给出了两种方法的结果（数据由MATLAB给出，具16位有效数字）：求解方法GaussCholesky所得结果(n=2)1.0000000000000000.9999999999999991.0000000000000000.999999999999999残差向量的模8.950904182623619e-0168.950904182623619e-016求解方法GaussCholesky所得结果(n=3)0.9999999999999991.0000000000000070.9999999999999930.9999999999999

8、991.0000000000000070.999999999999993残差向量的模1.017354490604556e-0141.018383803549894e-014求解方法GaussCholesky所得结果(n=4)0.9999999999999921.0000000000000930.9999999999997671.0000000000001550.9999999999999841.0000000000001830.9999999999995551.000000000000291残差向量的模2.948931463388727e-0135.624264157788779e-013求解

9、方法GaussCholesky所得结果(n=5)1.0000000000000210.9999999999995661.0000000000019740.9999999999969211.0000000000015390.9999999999999441.0000000000009440.9999999999962311.0000000000053750.999999999997485残差向量的模3.991872242690353e-0127.092940334727893e-012求解方法GaussCholesky所得结果(n=6)0.9999999999991271.00000000002

10、47390.9999999998333851.0000000004319720.9999999995243521.0000000001870690.9999999999991221.0000000000248190.9999999998331141.0000000004321450.9999999995246151.000000000186823残差向量的模6.900785737295254e-0106.900074426813605e-010求解方法GaussCholesky所得结果(n=7)0.9999999999916741.0000000003357320.99999999674113

11、51.0000000127418880.9999999765354281.0000000203493500.9999999932984200.9999999999924901.0000000003001960.9999999971047461.0000000112648190.9999999793355711.0000000178649620.999999994131788残差向量的模3.439014525977920e-0083.026516485826501e-008求解方法GaussCholesky所得结果(n=8)0.9999999999841001.0000000008929750.

12、9999999879340771.0000000670249000.9999998157945041.0000002650319100.9999998087679421.0000000545851260.9999999999675241.0000000017289510.9999999775345581.0000001211221280.9999996748657521.0000004589680020.9999996740067871.000000091828723残差向量的模3.851772775183292e-0076.680166553264907e-007求解方法GaussChole

13、sky所得结果(n=9)0.9999999998607441.0000000094215680.9999998428809111.0000011087208520.9999959709555481.0000081639292520.9999906839036481.0000055965869100.9999986236556250.9999999997453271.0000000174739900.9999997053727021.0000020975252130.9999923221015131.0000156522543460.9999820465351761.00001083322950

14、40.999997325588041残差向量的模1.428785278708662e-0052.748214245580457e-005求解方法GaussCholesky所得结果(n=10)0.9999999987955721.0000001017478450.9999978707282451.0000190829781950.9999100464682201.0002448322198400.9996017119706871.0003820703300710.9998007069587911.0000435784947960.9999999987740021.0000001052316450

15、.9999977695878921.0000201993121760.9999039495053141.0002633708062710.9995688155794111.0004159244407610.9997819910881091.000047876491003残差向量的模6.439185125919934e-0046.983810322355236e-004求解方法GaussCholesky所得结果(n=11)0.9999999979208121.0000002081628790.9999947943937111.0000564289852170.9996726231479251.0

16、011247814128860.9976001184982171.0032133911497670.9973734030220441.0011977749672620.9997664775668590.9999999951692511.0000005083622250.9999867482629321.0001488023899480.9991101023169511.0031394431987040.9931436743770921.0093723186318400.9921961469473491.0036184267743290.999283830342860残差向量的模0.005084

17、1954515620.014833271077307求解方法GaussCholesky所得结果(n=12)0.9999999740447491.0000033345813260.9998938518465541.0014618317703400.989181965552700.0479*0.8654683112560501.2451215670959650.7109378523452331.2128190000016150.9110958192036081.0160871651803650.9999999658809671.0000043164490560.9998643807090501.0

18、018469453411690.986462994827836.0594*0.8342783100850441.3000033295320840.6482460048728421.2576477080114060.8928649388775641.019305272304730残差向量的模0.4664883285306730.568465973687995求解方法GaussCholesky所得结果(n=13)0.9999999858198681.0000019397545770.9999325415586561.0010380480332600.991231595703128.0453*0.8

19、477605925415531.3424084676946830.4797262456708821.5271413115093440.6588049297866821.1275928263033610.9790298325752931.0000001008009660.9999837710632061.0006372056396670.989244282315779.0976*0.4659772827734422.876103771702808-3.3756417292152417.847802305139874-6.1082100679614845.693767548924378-0.784

20、7130303124171.297431751763950残差向量的模0.90788205110385212.070004065947840求解方法GaussCholesky所得结果(n=14)0.9999994824396171.0000680296145600.997857975573164.027*0.8331431608765061.4061015628336891.692627958831394-6.86006536429095726.482050799360898-45.46278765625567853.342818634605649-35.20710419073989815.1

21、39*5737-1.391730391291961残差向量的模84.421272602931595求解方法GaussCholesky所得结果(n=15)0.9999999884848880.9999988200010381.0001530354931710.995678337603479.0552*0.6068238310286722.717791261290424-3.8114192782710139.689857910094865-8.6574635145240826.3443761205678101.692954637371676.0796*2.8209856866996310.6247

22、25651254457残差向量的模15.395524776711373值得注意的是，当n13时，MATLAB认为矩阵此时非正定，无法进行Cholesky分解，说明此时矩阵病态十分严重。我们用残差向量的模来估计计算误差，然后做出误差对数与矩阵阶数的关系图如下：由图可知，n较小时两种方法的误差基本相同，但当n较大时，Cholesky方法由于受矩阵病态性的影响，误差急剧增大，甚至认为矩阵非正定而无法给出计算结果。2、 计算各阶Hilbert矩阵的条件数计算一系列Hilbert矩阵条件数的程序如下：formatlong,n=input(n=),fori=1:nC(i)=cond(hilb(i);end

23、C我们得到的结果如下表所示（保留三位有效数字）：矩阵阶数条件数矩阵阶数条件数矩阵阶数条件数11.0061.50E7115.23E14219.374.75E8121.68E16352481.53E10131.76E1841.55E494.93E11143.08E1754.77E5101.60E13154.43E17将矩阵条件数对矩阵阶数作图，如下所示：由图可知，当n较小时，矩阵条件数的对数随矩阵阶数n的增大而线性增加，而当n较大时，则增长幅度明显下降，并开始出现平台。3、 利用正则化方法求解原方程组的研究我们使用如下程序来实现正则化方法：format long;n=input(n=); %输入矩

24、阵阶数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=AnBn %求解正则化方程An*yn3=Bn的解errn3=(yn3-Xn).2 %计算残差平方向量sigma3=sqrt(sum(errn3) %计算新解相对Xn的标准偏差在(a)中，我们已经看到，对于n10的情形，由于Hilbert矩阵病态过于严重，矩阵的条件数过大，我们得到的解甚至无法给出3位有

25、效数字；而现在经过正则化处理之后，情况得到了明显的改善，为节约篇幅仅试举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.980817121

26、011882sigma3 = 0.032438344028382显然，此时的结果大大改观，在n=15时仍能给出两位有效数字。为更清楚地说明问题，我们计算了各正则化An矩阵的条件数，所用程序如下：format long;n=input(n=);a=input(a=);fori=1:nHi=hilb(i);Ii=eye(i);Ai=a.*Ii+Hi*Hi;C(i)=cond(Ai);endC程序给出的结果如下表所示（仍保留3位有效数字，其中a=10-7）：矩阵阶数条件数矩阵阶数条件数矩阵阶数条件数11.0062.62E7113.15E7237272.76E7123.22E732.71E582.88

27、E7133.29E742.06E792.98E7143.35E752.46E7103.07E7153.41E7我们继续使用图来表示条件数的对数与矩阵阶数之间的关系：显然，当使用了正则化方法之后，新矩阵明显优于原Hilbert矩阵：其条件数在n4时进入一个稳定的平台，基本保持不变，因此我们在n较大的时候仍能得到稳定的解。问题看上去已经得到了圆满解决，但其实并未结束：我们在使用正则化方法时，使用了近似，引入了新的参量a，那么a取何值时能够取得最好的效果？换言之，引入a能够很好地消除舍入误差的影响，但它同时又给我们的计算带来了截断误差，那么a为多少时能够在这两方面得到平衡？于是，我们继续编写程序进行

28、研究，研究对象是15阶的Hilbert矩阵：format long;n=input(n=);fori=1:na=10(-i);H=hilb(15);X=linspace(1,1,15);B=H*H*X;I=eye(15);Ai=a.*I+H*H;yi=AiB;erri=(yi-X).2;S(i)=sqrt(sum(erri);endS程序给出的结果如下表所示（保留3位有效数字）：a的负对数残差向量模a的负对数残差向量模a的负对数残差向量模11.10566.64E-2113.93E-320.66973.24E-2122.02E-330.32381.89E-2136.44E-340.21391.20E-2147.68E-259.82E-2105.04E-3150.684之后作图如下：由图可见，随着a的减小，误差先减后增，当a=10-12时，残差