清华大学数值分析实验报告.doc

.数值分析实验报告 一、 实验3.1题目：考虑线性方程组，编制一个能自动选取主元，又能手动选取主元的求解线性代数方程组的Gauss消去过程。（1）取矩阵，则方程有解。取计算矩阵的条件数。分别用顺序Gauss消元、列主元Gauss消元和完全选主元Gauss消元方法求解，结果如何？（2）现选择程序中手动选取主元的功能，每步消去过程都选取模最小或按模尽可能小的元素作为主元进行消元，观察并记录计算结果，若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。（4）选取其他你感兴趣的问题或者随机生成的矩阵，计算其条件数，重复上述实验，观察记录并分析实验的结果。1. 算法介绍首先，分析各种算法消去过程的计算公式，顺序高斯消去法：第k步消去中，设增广矩阵中的元素（若等于零则可以判定系数矩阵为奇异矩阵，停止计算），则对k行以下各行计算，分别用乘以增广矩阵的第行并加到第行，则可将增广矩阵中第列中以下的元素消为零；重复此方法，从第1步进行到第n-1步，则可以得到最终的增广矩阵，即；列主元高斯消去法：第k步消去中，在增广矩阵中的子方阵中，选取使得，当时，对中第行与第行交换，然后按照和顺序消去法相同的步骤进行。重复此方法，从第1步进行第n-1步，就可以得到最终的增广矩阵，即；完全主元高斯消去法：第k步消去中，在增广矩阵中对应的子方阵中，选取使得，若或，则对中第行与第行、第列与第列交换，然后按照和顺序消去法相同的步骤进行即可。重复此方法，从第1步进行到第n-1步，就可以得到最终的增广矩阵，即；接下来，分析回代过程求解的公式，容易看出，对上述任一种消元法，均有以下计算公式：2. 实验程序的设计一、输入实验要求及初始条件；二、计算系数矩阵A的条件数及方程组的理论解；三、对各不同方法编程计算，并输出最终计算结果。3. 计算结果及分析（1）先计算系数矩阵的条件数，结果如下，可知系数矩阵的条件数较大，故此问题属于病态问题， b或A的扰动都可能引起解的较大误差；采用顺序高斯消去法，计算结果为：最终解为x=(1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000001, 0.999999999999998, 1.000000000000004, 0.999999999999993, 1.000000000000012, 0.999999999999979, 1.000000000000028)T使用无穷范数衡量误差，得到=2.842170943040401e-14，可以发现，采用顺序高斯消元法求得的解与精确解之间误差较小。通过进一步观察，可以发现，按照顺序高斯消去法计算时，其选取的主元值和矩阵中其他元素大小相近，因此顺序高斯消去法方式并没有对结果造成特别大的影响。 若采用列主元高斯消元法，则结果为：最终解为x=(1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000)T同样使用无穷范数衡量误差，有=0； 若使用完全主元高斯消元法，则结果为最终解x=(1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000000)T同样使用无穷范数衡量误差，有=0；（2）若每步都选取模最小或尽可能小的元素为主元，则计算结果为最终解x=(1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000001 0.999999999999998 1.000000000000004 0.999999999999993 1.000000000000012 0.999999999999979 1.000000000000028)T使用无穷范数衡量误差，有为2.842170943040401e-14；而完全主元消去法的误差为=0。从（1）和（2）的实验结果可以发现，列主元消去法和完全主元消去法都得到了精确解，而顺序高斯消去法和以模尽量小的元素为主元的消去法没有得到精确解。在后两种消去法中，由于程序计算时的舍入误差，对最终结果产生了一定的影响，但由于方程组的维度较低，并且元素之间相差不大，所以误差仍比较小。为进一步分析，计算上述4种方法每步选取的主元数值，并列表进行比较，结果如下：第n次消元顺序列主元完全主元模最小16.000000000000000886.00000000000000024.666666666666667884.66666666666666734.285714285714286884.28571428571428644.133333333333333884.13333333333333354.064516129032258884.06451612903225864.031746031746032884.03174603174603274.015748031496063884.01574803149606384.007843137254902884.00784313725490294.003913894324853884.003913894324853104.0019550342130990.0156173706054690.0156173706054694.001955034213099从上表可以发现，对这个方程组而言，顺序高斯消去选取的主元恰好事模尽量小的元素，而由于列主元和完全主元选取的元素为8，与4在数量级上差别小，所以计算过程中的累积误差也较小，最终4种方法的输出结果均较为精确。在这里，具体解释一下顺序法与模最小法的计算结果完全一致的原因。该矩阵在消元过程中，每次选取主元的一列只有两个非零元素，对角线上的元素为4左右，而其正下方的元素为8，该列其余位置的元素均为0。在这样的情况下，默认的主元也就是该列最小的主元，因此两种方法所得到的计算结果是一致的。理论上说，完全高斯消去法的误差最小，其次是列主元高斯消去法，而选取模最小的元素作为主元时的误差最大，但是由于方程组的特殊性（元素相差不大并且维度不高），这个理论现象在这里并没有充分体现出来。（3）时，重复上述实验过程，各种方法的计算结果如下所示，在这里，仍采用无穷范数衡量绝对误差。顺序高斯消去法列主元高斯消去完全主元高斯消去选取模最小或尽可能小元素作为主元消去X1.000000000000000 1.000000000000000 1.000000000000000 1.000000000000001 0.999999999999998 1.000000000000004 0.999999999999993 1.000000000000014 0.999999999999972 1.000000000000057 0.9999999999998861.000000000000227 0.999999999999547 1.000000000000902 0.999999999998209 1.000000000003524 0.999999999993179 1.000000000012732 0.999999999978173 1.00000000002910211111111111111111111111111111111111111111.000000000000000 1.000000000000000 1.000000000000000 1.000000000000001 0.999999999999998 1.000000000000004 0.999999999999993 1.000000000000014 0.999999999999972 1.000000000000057 0.9999999999998861.000000000000227 0.999999999999547 1.000000000000902 0.999999999998209 1.000000000003524 0.999999999993179 1.000000000012732 0.999999999978173 1.0000000000291022.910205409989430e-11002.910205409989430e-11 可以看出，此时列主元和完全主元的计算结果仍为精确值，而顺序高斯消去和模尽可能小方法仍然产生了一定的误差，并且两者的误差一致。与n=10时候的误差比相比，n=20时的误差增长了大约1000倍，这是由于计算过程中舍入误差的不断累积所致。所以，如果进一步增加矩阵的维数，应该可以看出更明显的现象。（4）不同矩阵维度下的误差如下，在这里，为方便起见，选取2-条件数对不同维度的系数矩阵进行比较。维度条件数顺序消去列主元完全主元模尽量小1.7e+32.84e-14002.84e-141.8e+62.91e-11002.91e-115.7e+79.31e-10009.31e-101.8e+92.98e-08002.98e-081.9e+123.05e-05003.05e-053.8e+163.28e+043.88e-123.88e-123.28e+048.5e+163.52e+134.2e-34.2e-33.52e+13 从上表可以看出，随着维度的增加，不同方法对计算误差的影响逐渐体现，并且增长较快，这是由于舍入误差逐步累计而造成的。不过，方法二与方法三在维度小于40的情况下都得到了精确解，这两种方法的累计误差远比方法一和方法四慢；同样地，出于与前面相同的原因，方法一与方法四的计算结果保持一致，方法二与方法三的计算结果保持一致。4. 结论本文矩阵中的元素差别不大，模最大和模最小的元素并没有数量级上的差异，因此，不同的主元选取方式对计算结果的影响在维度较低的情况下并不明显，四种方法都足够精确。对比四种方法，可以发现采用列主元高斯消去或者完全主元高斯消去法，可以尽量抑制误差，算法最为精确。不过，对于低阶的矩阵来说，四种方法求解出来的结果误差均较小。另外，由于完全选主元方法在选主元的过程中计算量较大，而且可以发现列主元法已经可以达到很高的精确程度，因而在实际计算中可以选用列主元法进行计算。附录：程序代码clearclc;format long;%方法选择n=input(矩阵A阶数:n=);disp(选取求解方式);disp(1 顺序Gauss消元法，2 列主元Gauss消元法，3 完全选主元Gauss消元法，4 模最小或近可能小的元素作为主元);a=input(求解方式序号：);%赋值A和bA=zeros(n,n);b=zeros(n,1);for i=1:n A(i,i)=6; if i1 A(i,i-1)=8; end if i1，所以这是一个病态问题。由上表可以看出，四种方法的求解都存在一定的误差。下面分析误差的来源：LU分解方法的误差存在主要是由于Hilbert矩阵各元素由分数形式转换为小数形式时，不能除尽情况下会出现舍入误差，在进行LU分解时也存在这个问题，所以最后得到的结果不是方程的精确解 ，但结果显示该方法的误差非常小；Jacobi迭代矩阵的谱半径为4.30853，故此迭代法不收敛；GS迭代法在迭代次数为541次时得到了方程的近似解，其误差约为0.05 ，比较大。GS迭代矩阵的谱半径为0.999998，很接近1，所以GS迭代法收敛速度较慢；SOR迭代法在迭代次数为100次时误差约为0.08，误差较大。SOR迭代矩阵的谱半径为0.999999，也很接近1，所以时SOR迭代法收敛速度不是很快，但是相比于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.9978292219453491.0375262031068390.8969732619760151.0203451363750361.0690711669325761.0511799950366120.9968147571853640.9263432373255361.0129389722756340.9397138368551710.9882618050730811.0646370905351541.0836333450939741.0450601771155140.9706030247784690.880212649657655迭代次数-356100谱半径-6.0421310.999999999208776-时，算法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.9974420733067511.0190699093584090.9922782477867390.9564418583132370.9864203333613531.0213016119565911.0387010268066081.0359427734985331.0166937631494220.9857164549462500.9471812875006971.0157760397865720.9664291470644830.9286748681579100.9969315484827271.0667378039135371.0977924305964681.0880304408550691.0481106208111920.9899194185724240.9228408137041420.853252417221922迭代次数-1019100谱半径-8.6496410.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.532393125791018 -7.355047567109081 7.380667063930484-1.129041418095142 0.425748747257065 1.733284233971601 0.817952344733362-不收敛1.0043857406415901.0463460678775540.9071783477077290.9057634559490530.9725218027884571.0437314453679031.0915351694487641.1100900207039441.1031296846797681.0771686511460561.0385147362651760.9922599908320410.9421513904780030.8907853666840650.839876442493220迭代次数-262100谱半径-6.0421311.0000000000000008.355047567109082-0.160123557506780分析以上结果可以发现，随着n值的增加，Gauss消去法误差逐渐增大，而且误差增大的速度很快，在维数小于等于10情况下，Gauss消去法得到的结果误差较小；但当维数达到15时，计算结果误差已经达到精确解的很多倍；J法迭代不收敛，无论n如何取值，其谱半径始终大于1，因而J法不收敛，所以J迭代法不能用于Hilbert矩阵的求解；对于GS迭代法和SOR迭代法，两种方法均收敛，GS迭代法是SOR迭代法松弛因子取值为1的特例，SOR方法受到取值的影响，会有不同的收敛情况。可以得出GS迭代矩阵的谱半径小于1但是很接近1，收敛速度很慢。虽然随着维数的增大，所需迭代的次数逐渐减少，但是当维数达到15的时候，GS法已经不再收敛。因此可以得出结论，GS迭代方法在Hilbert矩阵维数较低时，能够在一定程度上满足迭代求解的需求，不过迭代的速度很慢。另外，随着矩阵维数的增加， SOR法的误差水平基本稳定，而且误差在可以接受的范围之内。经过比较可以得出结论，如果求解较低维度的Hibert矩阵问题，Gauss消去法、GS迭代法和SOR迭代法均可使用，且Gauss消去法的结果精确度较高；如果需要求解较高维度的Hibert矩阵问题，只有采用SOR迭代法。（3）系数矩阵的条件数较大时，为病态方程。由实验可知，Gauss法在解上述方程时，结果存在很大的误差。而对于收敛的迭代法，可以通过选取最优松弛因子的方法来求解，虽然迭代次数相对较多，但是结果较为精确。总体来看，对于一般病态方程组的求解，可以采用以下方式：1. 低维度下采用Gauss消去法直接求解是可行的；Jacobi迭代方法不适宜于求解病态问题；GS迭代方法可以解决维数较低的病态问题，但其谱半径非常趋近于1，导致迭代算法收敛速度很慢，维数较大的时候，GS法也不再收敛；SOR方法较适合于求解病态问题，特别是矩阵维数较高的时候，其优势更为明显。2. 采用高精度的运算，如选用双倍或更多倍字长的运算，可以提高收敛速度；3. 可以对原方程组作某些预处理，从而有效降低系数矩阵的条件数。4. 实验结论（1）对Hibert矩阵问题，其条件数会随着维度的增加迅速增加，病态性会越来越明显；在维度较低的时候，Gauss消去法、GS迭代法和SOR迭代法均可使用，且可以优先使用Gauss消去法；如果需要求解较高维度的Hibert矩阵问题，只有SOR迭代法能够求解。（2）SOR方法比较适合于求解病态问题，特别是矩阵维数较高的时候，其优点更为明显。从本次实验可以看出，随着矩阵维数的增大，SOR方法所需的迭代次数减少，而且误差基本稳定，是解决病态问题的适宜方法。附录：程序代码clear allclc;format long;%矩阵赋值n=input(矩阵H的阶数:n=);for i=1:n for j=1:n H(i,j)=1/(i+j-1); endendb=H*ones(n,1);disp(H矩阵为：);Hdisp(向量b：);b%方法选择disp(选取求解方式);disp(1 Gauss消去法，2 J迭代法，3 GS迭代法，4 SOR迭代法);a=input(求解方式序号：);%Gauss消去法if a=1; H1=H;b1=b; for k=1:n if H1(k,k)=0 disp(主元为零，Gauss消去法无法进行); break end fprintf(第%d次消元所选取的主元是：%gn,k,H1(k,k) for p=k+1:n m5=-H1(p,k)/H1(k,k); H1(p,k:n)=H1(p,k:n)+m5*H1(k,k:n); b1(p)=b1(p)+m5*b1(k); end endx1(n)=b1(n)/H1(n,n);for k=n-1:-1:1 for v=k+1:n b1(k)=b1(k)-H1(k,v)*x1(v); end x1(k)=b1(k)/H1(k,k);end disp(Gauss消去法解为:); disp(x1); disp(解与精确解之差的无穷范数); norm(x1-a),inf)endD=diag(diag(H);L=-tril(H,-1);U=-triu(H,1); %J迭代法if a=2; %给定初始x0 ini=input(初始值设定:x0=); x0(:,1)=ini*diag(ones(n); disp(初始解向量为:); x0 xj(:,1)=x0(:,1); B=(D(-1)*(L+U); f=(D(-1)*b; fprintf(J法B矩阵谱半径为：%gn,vrho(B); if vrho(B)1; for m2=1:5000 xj(:,m2+1)=B*xj(:,m2)+fj; if norm(xj(:,m2+1)-xj(:,m2),inf)0.0001 break end end disp(J法计算结果为:); xj(:,m2+1) disp(解与精确解之差的无穷范数); norm(xj(:,m2+1)-diag(ones(n),inf) disp(J迭代法迭代次数:); m2 else disp(由于B矩阵谱半径大于1，因而J法不收敛); endend %GS迭代法if a=3; %给定初始x0 ini=input(初始值设定:x0=); x0(:,1)=ini*diag(ones(n); disp(初始解向量为:); x0 xG(:,1)=x0(:,1); G=inv(D-L)*U; fG=inv(D-L)*b; fprintf(GS法G矩阵谱半径为：%gn,vrho(G); if vrho(G)1 for m3=1:5000 xG(:,m3+1)=G*xG(:,m3)+fG; if norm(xG(:,m3+1)-xG(:,m3),inf)0.0001 break; end end disp(GS迭代法计算结果:); xG(:,m3+1) disp(解与精确解之差的无穷范数); norm(xG(:,m3+1)-diag(ones(n),inf) disp(GS迭代法迭代次数:); m3 else disp(由于G矩阵谱半径大于1，因而GS法不收敛); end end %SOR迭代法 if a=4; %给定初始x0 ini=input(初始值设定:x0=); x0(:,1)=ini*diag(ones(n); disp(初始解向量为:); x0 A=H; for i=1:n b(i)=sum(A(i,:); end x_star=ones(n,1); format long w=input(松弛因子:w=); Lw=inv(D-w*L)*(1-w)*D+w*U); f=w*inv(D-w*L)*b; disp(迭代矩阵的谱半径：) p=vrho(Lw) time_max=100;%迭代次数 x=zeros(n,1);%迭代初值 for i=1:time_max x=Lw*x+f; end disp(SOR迭代法得到的解为); x disp(解与精确解之差的无穷范数)