2022年清华大学数值分析实验报告.docx

《2022年清华大学数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《2022年清华大学数值分析实验报告.docx（39页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、2022年清华大学数值分析实验报告 数值分析试验报告 一、 试验3.1 题目： 考虑线性方程组，编制一个能自动选取主元，又能手动选取主元的求解线性代数方程组的Gauss消去过程。 （1）取矩阵，则方程有解。取计算矩阵的条件数。分别用依次Gauss消元、列主元Gauss消元和完全选主元Gauss消元方法求解，结果如何？ （2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，视察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析试验的结果。 （3）取矩阵阶数n=20或者更大，重复上述试验过程，视察记录并分析不同的问题及消去过程

2、中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。 （4）选取其他你感爱好的问题或者随机生成的矩阵，计算其条件数，重复上述试验，视察记录并分析试验的结果。 1. 算法介绍 首先，分析各种算法消去过程的计算公式， 依次高斯消去法： 第k步消去中，设增广矩阵中的元素（若等于零则可以判定系数矩阵为奇异矩阵，停止计算），则对k行以下各行计算，分别用乘以增广矩阵的第行并加到第行，则可将增广矩阵中第列中以下的元素消为零；重复此方法，从第1步进行到第n-1步，则可以得到最终的增广矩阵，即； 列主元高斯消去法： 第k步消去中，在增广矩阵中的子方阵中，选取使得，当时，对中第行与第行交换，然后

3、根据和依次消去法相同的步骤进行。重复此方法，从第1步进行第n-1步，就可以得到最终的增广矩阵，即； 完全主元高斯消去法： 第k步消去中，在增广矩阵中对应的子方阵中，选取使得，若或，则对中第行与第行、第列与第列交换，然后根据和依次消去法相同的步骤进行即可。重复此方法，从第1步进行到第n-1步，就可以得到最终的增广矩阵，即； 接下来，分析回代过程求解的公式，简单看出，对上述任一种消元法，均有以下计算公式： 2. 试验程序的设计 一、输入试验要求及初始条件； 二、计算系数矩阵A的条件数及方程组的理论解； 三、对各不同方法编程计算，并输出最终计算结果。 3. 计算结果及分析 （1） 先计算系数矩阵的条

4、件数，结果如下， 可知系数矩阵的条件数较大，故此问题属于病态问题， b或A的扰动都可能引起解的较大误差； 采纳依次高斯消去法，计算结果为： 最终解为x=T 运用无穷范数衡量误差，得到=2.842170943040401e-14，可以发觉，采纳依次高斯消元法求得的解与精确解之间误差较小。通过进一步视察，可以发觉，根据依次高斯消去法计算时，其选取的主元值和矩阵中其他元素大小相近，因此依次高斯消去法方式并没有对结果造成特殊大的影响。 若采纳列主元高斯消元法，则结果为： 最终解为x=T 同样运用无穷范数衡量误差，有=0； 若运用完全主元高斯消元法，则结果为 最终解x=T 同样运用无穷范数衡量误差，有=

5、0； （2） 若每步都选取模最小或尽可能小的元素为主元，则计算结果为 最终解x=T 运用无穷范数衡量误差，有为2.842170943040401e-14；而完全主元消去法的误差为=0。 从（1）和（2）的试验结果可以发觉，列主元消去法和完全主元消去法都得到了精确解，而依次高斯消去法和以模尽量小的元素为主元的消去法没有得到精确解。在后两种消去法中，由于程序计算时的舍入误差，对最终结果产生了肯定的影响，但由于方程组的维度较低，并且元素之间相差不大，所以误差仍比较小。 为进一步分析，计算上述4种方法每步选取的主元数值，并列表进行比较，结果如下： 第n次消元 依次 列主元 完全主元 模最小 1 6.0

6、00000000000000 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.0078431

7、37254902 8 8 4.007843137254902 9 4.003913894324853 8 8 4.003913894324853 10 4.001955034213099 0.015617370605469 0.015617370605469 4.001955034213099 从上表可以发觉，对这个方程组而言，依次高斯消去选取的主元恰好事模尽量小的元素，而由于列主元和完全主元选取的元素为8，与4在数量级上差别小，所以计算过程中的累积误差也较小，最终4种方法的输出结果均较为精确。 在这里，详细说明一下依次法与模最小法的计算结果完全一样的缘由。该矩阵在消元过程中，每次选取主元的一

8、列只有两个非零元素，对角线上的元素为4左右，而其正下方的元素为8，该列其余位置的元素均为0。在这样的状况下，默认的主元也就是该列最小的主元，因此两种方法所得到的计算结果是一样的。 理论上说，完全高斯消去法的误差最小，其次是列主元高斯消去法，而选取模最小的元素作为主元时的误差最大，但是由于方程组的特别性（元素相差不大并且维度不高），这个理论现象在这里并没有充分体现出来。 （3） 时，重复上述试验过程，各种方法的计算结果如下所示，在这里，仍采纳无穷范数衡量肯定误差。 依次高斯消去法 列主元高斯消去 完全主元高斯消去 选取模最小或尽可能小元素作为主元消去 X 1.000000000000000 1.

9、000000000000000 1.000000000000000 1.000000000000001 0.999999999999998 1.000000000000004 0.999999999999993 1.000000000000014 0.999999999999972 1.000000000000057 0.999999999999886 1.000000000000227 0.999999999999547 1.000000000000902 0.999999999998209 1.000000000003524 0.999999999993179 1.000000000012

10、732 0.999999999978173 1.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.000000000000000 1.000000000000000 1.000000000000000 1.000000000000001 0.999999999999998 1.000000000000004 0.999999999999993 1.000000000000014 0.999999999999972 1.000000000000057

11、0.999999999999886 1.000000000000227 0.999999999999547 1.000000000000902 0.999999999998209 1.000000000003524 0.999999999993179 1.000000000012732 0.999999999978173 1.000000000029102 2.910205409989430e-11 0 0 2.910205409989430e-11 可以看出，此时列主元和完全主元的计算结果仍为精确值，而依次高斯消去和模尽可能小方法仍旧产生了肯定的误差，并且两者的误差一样。与n=10时候的误差

12、比相比，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

13、+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. 结论 本文矩阵中的元素差别不大，模最大和模最小的元素并没有数量级上的差异，因此，不同的主元选取方式对计算结果的影响在维度

14、较低的状况下并不明显，四种方法都足够精确。 对比四种方法，可以发觉采纳列主元高斯消去或者完全主元高斯消去法，可以尽量抑制误差，算法最为精确。不过，对于低阶的矩阵来说，四种方法求解出来的结果误差均较小。 另外，由于完全选主元方法在选主元的过程中计算量较大，而且可以发觉列主元法已经可以达到很高的精确程度，因而在实际计算中可以选用列主元法进行计算。 附录：程序代码 clear clc; format long; %方法选择 n=input; disp; disp; a=input; %赋值A和b A=zeros; b=zeros; for i=1:n A=6; if i>1 A=8; end

15、if i<n A=1; end end for i=1:n for j=1:n b=b+A; end end disp; A disp; b %求条件数及理论解 disp; X= fprintf); fprintf); fprintf); %依次Gauss消元法 if a=1 A1=A;b1=b; for k=1:n if A1=0 disp; break end fprintf) %disp; %A1 for p=k+1:n l=A1/A1; A1=A1-l*A1; b1=b1-l*b1; end end x1=b1/A1; for k=n-1:-1:1 for w=k+1:n b1=

16、b1-A1*x1; end x1=b1/A1; end disp; disp; disp; norm end %列主元Gauss消元法 if a=2 A2=A;b2=b; for k=1:n max_i,max_j=find=max); if max_i=k A2_change=A2; A2=A2; A2=A2_change; b2_change=b2; b2=b2; b2=b2_change; end if A2=0 disp; break end fprintf) %disp; %A2 for p=k+1:n l=A2/A2; A2=A2-l*A2; b2=b2-l*b2; end end

17、 x2=b2/A2; for k=n-1:-1:1 for w=k+1:n b2=b2-A2*x2; end x2=b2/A2; end disp; disp; disp; norm end %完全选主元Gauss消元法 if a=3 A3=A;b3=b; for k=1:n VV=eye; max_i,max_j=find=max); if numel=0 max_i,max_j=find=-max); end W=eye; W+k-1,max_i+k-1)=0; W=0; W+k-1,k)=1; W+k-1)=1; V=eye; V=0; V+k-1,max_j+k-1)=0; V+k-1

18、)=1; V+k-1,k)=1; A3=W*A3*V; b3=W*b3; VV=VV*V; if A3=0 disp; break end fprintf) %disp; %A3 for p=k+1:n l=A3/A3; A3=A3-l*A3; b3=b3-l*b3; end end x3=b3/A3; for k=n-1:-1:1 for w=k+1:n b3=b3-A3*x3; end x3=b3/A3; end disp; disp; disp; norm end %模最小或近可能小的元素作为主元 if a=4 A4=A;b4=b; for k=1:n AA=A4; AA=NaN; mi

19、n_i,j=find=min); if numel=0 min_i,j=find=-min); end W=eye; W+k-1,min_i+k-1)=0; W=0; W+k-1,k)=1; W+k-1)=1; A4=W*A4; b4=W*b4; if A4=0 disp; break end fprintf) %A4 for p=k+1:n l=A4/A4; A4=A4-l*A4; b4=b4-l*b4; end end x4=b4/A4; for k=n-1:-1:1 for w=k+1:n b4=b4-A4*x4; end x4=b4/A4; end disp; disp; disp;

20、norm end 二、试验3.3 题目： 考虑方程组的解，其中系数矩阵H为Hilbert矩阵： 这是一个闻名的病态问题。通过首先给定解（例如取为各个重量均为1）再计算出右端的方法给出确定的问题。 （1）选择问题的维数为6，分别用Gauss消去法（即LU分解）、J迭代法、GS迭代法和SOR迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何。 （2）逐步增大问题的维数，仍用上述的方法来解它们，计算的结果如何？计算的结果说明的什么？ （3）探讨病态问题求解的算法。 1. 算法设计 对随意线性方程组，分析各种方法的计算公式如下， （1）Gauss消去法： 首先对系数矩阵进行LU分解

21、，有，则原方程转化为，令，则原方程可以分为两步回代求解： 详细方法这里不再赘述。 （2）J迭代法： 首先分解，再构造迭代矩阵，其中 ，进行迭代计算，直到误差满意要求。 （3）GS迭代法： 首先分解，再构造迭代矩阵 ，其中 ，进行迭代计算，直到误差满意要求。 （4）SOR迭代法： 首先分解，再构造迭代矩阵 ，其中，进行迭代计算，直到误差满意要求。 2. 试验过程 一、依据维度n确定矩阵H的各个元素和b的各个重量值； 二、选择计算方法（ Gauss消去法，J迭代法，GS迭代法，SOR迭代法），对迭代法设定初值，此外SOR方法还须要设定松弛因子； 三、进行计算，直至满意误差要求（对迭代法，设定相邻两

22、次迭代结果之差的无穷范数小于0.0001；对SOR方法，设定为输出迭代100次之后的结果及误差值），输出试验结果。 3. 计算结果及分析 （1）时，问题可以详细定义为 计算结果如下， Gauss消去法 第1次消元所选取的主元是：1 第2次消元所选取的主元是：0.0833333 第3次消元所选取的主元是：0.00555556 第4次消元所选取的主元是：0.000357143 第5次消元所选取的主元是：2.26757e-05 第6次消元所选取的主元是：1.43155e-06 解得X=T 运用无穷范数衡量误差，可得=4.254160357319847e-10； J迭代法 设定迭代初值为零，计算得到

23、J法的迭代矩阵B的谱半径为4.308531，所以J法不收敛； GS迭代法 设定迭代初值为零，计算得到GS法的迭代矩阵G的谱半径为：0.9999981，故GS法收敛，经过541次迭代计算后，结果为X=T 运用无穷范数衡量误差，有=0.047045569989162； SOR迭代法 设定迭代初值为零向量，并设定，计算得到SOR法迭代矩阵谱半径为0.999999433815223，经过100次迭代后的计算结果为 X=T； 运用无穷范数衡量误差，有=0.082326357506473； 对SOR方法，可变，变更值，计算结果可以列表如下 迭代次数 100 100 100 100 迭代矩阵的谱半径 0.9

24、99999433815223 0.999998867083155 0.999996830135013 0.999982309342386 X 1.003653917714694 0.974666041209353 1.011814573842440 1.042837929171827 1.017190220902681 0.945462001336268 1.014676015634604 0.896636864424096 1.090444578936265 1.107070542628148 1.006315452225331 0.873244842279255 1.02802221550

25、5147 0.790604920509843 1.267167365524072 1.061689730857891 0.990084054872602 0.846005956774467 1.051857392323966 0.653408758549156 1.486449891152510 0.783650360698119 1.349665420488270 0.664202350634588 0.054537998663732 0.126755157720745 0.267167365524072 0.486449891152510 可以发觉，松弛因子的取值对迭代速度造成了不同的影响

26、，上述四种方法中，松弛因子=0.5时，收敛相对较快。 综上，四种算法的结果列表如下： 算法 Gauss消去法 Jacobi法 GS法 SOR法（取） 迭代次数 - 不收敛 541 100 迭代矩阵的谱半径 - 4.30853 0.999998 0.999999433815223 X 0.999999999999228 1.000000000021937 0.999999999851792 1.000000000385369 0.999999999574584 1.000000000167680 - 1.001178105812706 0.999144082651860 0.9689290939

27、84902 1.047045569989162 1.027323158370281 0.954352032784608 1.003380614145078 0.962420297458423 1.031857023134559 1.061814901289881 1.014037815827164 0.917673642493527 4.254160357319847e-10 - 0.047045569989162 0.082326357506473 计算可得，矩阵H的条件数为>>1，所以这是一个病态问题。由上表可以看出，四种方法的求解都存在肯定的误差。下面分析误差的来源： LU分

28、解方法的误差存在主要是由于Hilbert矩阵各元素由分数形式转换为小数形式时，不能除尽状况下会出现舍入误差，在进行LU分解时也存在这个问题，所以最终得到的结果不是方程的精确解 ，但结果显示该方法的误差特别小； Jacobi迭代矩阵的谱半径为4.30853，故此迭代法不收敛； GS迭代法在迭代次数为541次时得到了方程的近似解，其误差约为0.05 ，比较大。GS迭代矩阵的谱半径为0.999998，很接近1，所以GS迭代法收敛速度较慢； SOR迭代法在迭代次数为100次时误差约为0.08，误差较大。SOR迭代矩阵的谱半径为0.999999，也很接近1，所以时SOR迭代法收敛速度不是很快，但是相比于

29、GS法，在迭代速度方面已经有了明显的提高；另外，对不同的，SOR方法的迭代速度会相应有改变，假如选用最佳松弛因子，可以实现更快的收敛； （2） 考虑不同维度的状况，时， 算法 Gauss消去 J法 GS法 SOR法（w=0.5） 计算结果 0.999999999966269 1.000000001809060 0.999999976372676 1.000000127868103 0.999999655764116 1.000000487042164 0.999999653427125 1.000000097774747 - 0.997829221945349 1.03752620310683

30、9 0.896973261976015 1.020345136375036 1.069071166932576 1.051179995036612 0.996814757185364 0.926343237325536 1.012938972275634 0.939713836855171 0.988261805073081 1.064637090535154 1.083633345093974 1.045060177115514 0.970603024778469 0.880212649657655 迭代次数 - - 356 100 谱半径 - 6.04213 1 0.99999999920

31、8776 - 时， 算法 Gauss消去法 Jacobi法 GS法 SOR法（w=0.5） 计算结果 0.999999994751197 1.000000546746354 0.999985868343700 1.000157549468631 0.999063537004329 1.003286333127805 0.992855789229370 1.009726486881556 0.991930155925812 1.003729850349020 0.999263885025643 - 0.997442073306751 1.019069909358409 0.99227824778

32、6739 0.956441858313237 0.986420333361353 1.021301611956591 1.038701026806608 1.035942773498533 1.016693763149422 0.985716454946250 0.947181287500697 1.015776039786572 0.966429147064483 0.928674868157910 0.996931548482727 1.066737803913537 1.097792430596468 1.088030440855069 1.048110620811192 0.98991

33、9418572424 0.922840813704142 0.853252417221922 迭代次数 - - 1019 100 谱半径 - 8.64964 1 0.999999999999966 - 时 算法 Gauss消去法 Jacobi法 GS法 SOR法（w=0.5） 计算结果 0.999999968723799 1.000002417094896 0.999994922439769 0.998640261957706 1.025668111139297 0.781933485305194 2.066840925345890 -2.279036697492128 7.532393125

34、791018 -7.355047567109081 7.380667063930484 -1.129041418095142 0.425748747257065 1.733284233971601 0.817952344733362 - 不收敛 1.004385740641590 1.046346067877554 0.907178347707729 0.905763455949053 0.972521802788457 1.043731445367903 1.091535169448764 1.110090020703944 1.103129684679768 1.0771686511460

35、56 1.038514736265176 0.992259990832041 0.942151390478003 0.890785366684065 0.839876442493220 迭代次数 - - 262 100 谱半径 - 6.04213 1 1.000000000000000 8.355047567109082 - - 0.160123557506780 分析以上结果可以发觉，随着n值的增加，Gauss消去法误差渐渐增大，而且误差增大的速度很快，在维数小于等于10状况下，Gauss消去法得到的结果误差较小；但当维数达到15时，计算结果误差已经达到精确解的许多倍； J法迭代不收敛，无论

36、n如何取值，其谱半径始终大于1，因而J法不收敛，所以J迭代法不能用于Hilbert矩阵的求解； 对于GS迭代法和SOR迭代法，两种方法均收敛，GS迭代法是SOR迭代法松弛因子取值为1的特例，SOR方法受到取值的影响，会有不同的收敛状况。可以得出GS迭代矩阵的谱半径小于1但是很接近1，收敛速度很慢。虽然随着维数的增大，所需迭代的次数渐渐削减，但是当维数达到15的时候，GS法已经不再收敛。因此可以得出结论，GS迭代方法在Hilbert矩阵维数较低时，能够在肯定程度上满意迭代求解的需求，不过迭代的速度很慢。另外，随着矩阵维数的增加， SOR法的误差水平基本稳定，而且误差在可以接受的范围之内。 经过比

37、较可以得出结论，假如求解较低维度的Hibert矩阵问题，Gauss消去法、GS迭代法和SOR迭代法均可运用，且Gauss消去法的结果精确度较高；假如须要求解较高维度的Hibert矩阵问题，只有采纳SOR迭代法。 （3） 系数矩阵的条件数较大时，为病态方程。由试验可知，Gauss法在解上述方程时，结果存在很大的误差。而对于收敛的迭代法，可以通过选取最优松弛因子的方法来求解，虽然迭代次数相对较多，但是结果较为精确。 总体来看，对于一般病态方程组的求解，可以采纳以下方式： 1. 低维度下采纳Gauss消去法干脆求解是可行的； Jacobi迭代方法不相宜于求解病态问题； GS迭代方法可以解决维数较低的

38、病态问题，但其谱半径特别趋近于1，导致迭代算法收敛速度很慢，维数较大的时候，GS法也不再收敛； SOR方法较适合于求解病态问题，特殊是矩阵维数较高的时候，其优势更为明显。 2. 采纳高精度的运算，如选用双倍或更多倍字长的运算，可以提高收敛速度； 3. 可以对原方程组作某些预处理，从而有效降低系数矩阵的条件数。 4. 试验结论 （1）对Hibert矩阵问题，其条件数会随着维度的增加快速增加，病态性会越来越明显；在维度较低的时候，Gauss消去法、GS迭代法和SOR迭代法均可运用，且可以优先运用Gauss消去法；假如须要求解较高维度的Hibert矩阵问题，只有SOR迭代法能够求解。 （2）SOR方

39、法比较适合于求解病态问题，特殊是矩阵维数较高的时候，其优点更为明显。从本次试验可以看出，随着矩阵维数的增大，SOR方法所需的迭代次数削减，而且误差基本稳定，是解决病态问题的相宜方法。 附录：程序代码 clear all clc; format long; %矩阵赋值 n=input; for i=1:n for j=1:n H=1/; end end b=H*ones; disp; H disp; b %方法选择 disp; disp; a=input; %Gauss消去法 if a=1; H1=H;b1=b; for k=1:n if H1=0 disp; break end fprintf

40、) for p=k+1:n m5=-H1/H1; H1=H1+m5*H1; b1=b1+m5*b1; end end x1=b1/H1; for k=n-1:-1:1 for v=k+1:n b1=b1-H1*x1; end x1=b1/H1; end disp; disp; disp; norm,inf) end D=diag); L=-tril; U=-triu; %J迭代法 if a=2; %给定初始x0 ini=input; x0=ini*diag); disp; x0 xj=x0; B=)*; f=)*b; fprintf); if vrho<1; for m2=1:5000 xj=B*xj+fj; if norm-xj),inf)<0.0001 break end end disp; xj disp; norm-diag),inf) disp;