数值代数第五次作业.doc

第一题：程序：1. gongetidu.mx=gongetidu(A,b,s,kmax)输入值：n*n的矩阵A ，n*1的矩阵b，以及控制迭代停止的s和最大允许的迭代次数kmax输出值：gongetidu法求解问题A*x=b的解x2. xldSOR.mx,w=xldSOR(A,b)输入值：n*n的矩阵A ，n*1的矩阵b输出值：超松弛迭代法迭代用法求解问题A*x=b的解x,以及最佳松弛因子w3.call1.m第一题的程序问题分析：1.方程分析将Dirichlet问题转化为如下差分形式：1+h24uij-14ui+1,j+ui-1,j+ui,j+1+ui,j-1=h24fi,j i,j=1,2,n-1ui,0=ui,0 ,ui,n=uj,n i=0,1,nu0,j=u0,j , un,j=un,j i=0,1,n将 ui,j ( i,j=1,2,n-1) 按照自然顺序排序并记向量x ，大小为 1 ,(n-1)2系数矩阵A为A= SBBSBSBBS 矩阵大小为(n-1)2, (n-1)2 , 其中B=-I/4 ,S是对角元均为 1+h2/4 ,次对角元均为-1/4向量b ,讨论两种情况：（1） 对于 i,j=2n-2 即 uij不在边界上 有 1+h24uij-14ui+1,j+ui-1,j+ui,j+1+ui,j-1=h24fi,j 即对应的b的分量为 h24fi,j=h24sin(i,j)（2）对于 uij 在边界上，例如 i=1, j=2n-2 ,有 ui-1,j=u0,j=p0,j=j2 ，有 1+h24uij-14ui+1,j+ui,j+1+ui,j-1=h24fi,j+14*j2即对应的b的分量为 h24fi,j+14*j2为此，可以用for循环同理定义矩阵b，令b=h24fi,j若i=1 , 对应的b的分量+14*j2若i=n-1 , 对应的b的分量+14*(j2+n2)若j=1 , 对应的b的分量+14*i2若j=n-1 , 对应的b的分量+14*(i+n2)再将矩阵b按照自然顺序排序得到向量b由此计算 A*x=b最后再将1*（n-1）2的向量x变化为（n-1）*（n-1）的矩阵。2.算法分析：1. 实共轭梯度法采用实共轭梯度法解对称正定方程组A*x=b，在|rk|很小或者迭代次数达到最大允许的迭代次数kmax来终止迭代，即令sqrt(rr)s|b|2或kkmax2. SOR迭代法：SOR迭代法需要确定最佳松弛因子wopt。选用不同的w值1:0.01:2，用相同的初始向量迭代400次，考虑 x400-x399的无穷范数，无穷范数最小的认为其收敛速度最快，作为最佳松弛因子wopt .再用这个w ,令xk=（D-wL）-1 (1-w)D+wU xk-1+w（D -wL）-1 b; 直到 xk- xk-1 的无穷范数0.0001, 将此时的xk 作为解。这样的解满足有四位有效数字。实验结果1. 用实共轭梯度法求解x=gongetidu(A,b,10(-6),1000)即在|rk|小于10(-6)*|b|或者迭代次数达到最大允许的迭代次数kmax=1000来终止迭代。得到的解x及r=A*x-b的前若干项为x208.4135.1259527104.129590.849945786.534676487.629956292.4930221100.260195110.413275122.587999b-Ax6.99E-06-8.77E-064.24E-05-3.40E-06-2.35E-05-2.29E-05-5.08E-05-2.61E-05-1.25E-05-9.38E-07残量r(x)的二范数|A*x-b|2=|r(x)|2= 8.7224e-004系统计算时间为Elapsed time is 0.026339 seconds.结果分析：实共轭梯度法求解速度非常快，这与书中的结论“只要线性方程组的系数矩阵是十分良态的，共轭梯度法就会收敛的非常快”相符合。可以通过增加迭代次数增加精度。将x重新排列成矩阵U，并画图得：2.用上节课的SOR迭代法进行求解。其中通过对不同的w分别迭代400次选取迭代速度最快的最为最佳松弛因子wopt。对wopt迭代指导精确到0.0001终止。得到x和b-Ax的前若干项为：x208.4005135.1266104.1390.8506286.5348787.6295492.49222100.2599110.4141122.5888b-A*x2.91E-064.42E-065.69E-066.71E-067.47E-067.99E-068.26E-068.30E-068.13E-067.76E-06残量r(x)的二范数|A*x-b|2=|r(x)|2= 1.5601e-004系统计算时间为Elapsed time is 9.345953 seconds.最佳松弛因子为1.5530，精切到小数点后两位结果分析：由于需要选取SOR迭代法的最佳松弛因子，迭代法需要的时间较长，但是精度较好。可以通过以下方式改进：先确定松弛因子的大概范围再精细计算以减少计算量3.同时用共轭梯度法和SOR迭代法进行对比求解为了方便比较，令两种迭代法的终止条件都设定为|rk|10(-10)*|b|或者kmax=2000, 比较两种方法在结果要求精度较高的情况下的计算效果。运算时间残差二范数|r|2共轭梯度0.0301866.3743e-008SOR 2.9351068.9492e-008结果分析：1）相较于SOR迭代法，共轭梯度法的运算速度非常快，精度高，是一种实用的迭代方法。2）SOR迭代法在计算最佳松弛因子时需要较长时间，松弛因子的计算精度越高，运算量越大。在迭代的最终结果要求并不非常高的时候这种方法并不适合。第二题：4.call.m第二题的程序5. myHilbert.mH=myHilbert(n)输入n，输出n阶hilbert方阵。实验结果：对hilbert矩阵 ai,j=1(i+j-1) ,bi=13j=1nai,j 求解Ax=b首先，易知真实解为1，1，1/31）对于不同的矩阵大小（n=10,20,200）,分别用共轭梯度法计算解，将计算解与真实解进行比较，并计算（以下列出n=10,20,30,40,50的结果）n=10n=20n=30n=40n=50计算时间0.0005000.0005200.0004930.0004980.000606|r|28.25e-052.56e-059.97e-050.0002263212.26e-05max|x-1/3|4.8258e-051.174441e-054.99781e-0511.7e-058.9917748e-06X0.33655 0.33169 0.32958 0.32718 0.33486 0.31701 0.34435 0.35320 0.36108 0.32240 0.33581 0.32708 0.32945 0.33358 0.34085 0.34480 0.32340 0.31984 0.31930 0.34339 0.34603 0.32598 0.31903 0.31504 0.33986 0.34277 0.33034 0.32217 0.31608 0.33523 0.33709 0.33453 0.32666 0.31961 0.33128 0.33017 0.33782 0.33126 0.32408 0.32849 0.32273 0.34001 0.33537 0.32867 0.32683 0.31515 0.34114 0.33877 0.33295 0.32612 0.34132 0.34138 0.33673 0.32613 0.34070 0.34323 0.33993 0.32667 0.33940 0.34437 0.34252 0.32758 0.33755 0.34488 0.34453 0.32871 0.33527 0.34484 0.34601 0.32997 0.33264 0.34432 0.34699 0.33128 0.32974 0.34339 0.34752 0.33257 0.32663 0.34211 0.34766 0.33381 0.32336 0.34054 0.34745 0.33497 0.31998 0.33871 0.34692 0.33602 0.33668 0.34613 0.33694 0.33447 0.34510 0.33775 0.33212 0.34386 0.33841 0.32966 0.34244 0.33895 0.32710 0.34086 0.33935 0.32447 0.33915 0.33963 0.32178 0.33733 0.33978 0.31906 0.33541 0.33981 0.31630 0.33340 0.33973 0.31352 0.33133 0.33954 0.32919 0.33925 0.32701 0.33886 0.32479 0.33838 0.32254 0.33782 0.32027 0.33718 0.31798 0.33647 0.31568 0.33569 0.31336 0.33484 0.31105 0.33394 0.30874 0.33299 0.33198 0.33093 0.32984 0.32872 0.32755 0.32636 0.32513 0.32389 0.32261 0.32132 2）随着矩阵阶层增大（n=10,20,200），其运算速度，计算精度比较：矩阵大小nmax|x-1/3|残差二范数|r|2运行时间t104.83E-058.25E-050.000551201.17E-052.56E-050.000545305.00E-059.97E-050.000663400.0001170.0002260.000513508.99E-062.26E-050.000613601.64E-053.94E-050.000833702.62E-056.09E-050.000679803.83E-058.73E-050.000718905.28E-050.0001180.0011211006.94E-050.0001530.0009691108.81E-050.0001930.0012791200.0001090.0002360.0008561300.0001310.0002830.0009081400.0001550.0003330.0010171500.0001810.0003860.0009861600.0002080.0004430.0009791700.0002360.0005030.0011171800.0002660.0005650.0011111900.0002970.000630.001072000.0003290.0006970.001192结果分析：1）随着矩阵阶层的增加，误差增大，运算时间增加，但是增加范围可接受。2）由上面的结果可用共轭梯度法求解这样的方程是合适的，计算速度快，精度较高；对比前几章的计算结果可得，相比于Gauss变换，选主元高斯变换，平方根法，改进的平方根法等方法，共和梯度法计算速度更快更准确。习题三：6.call3.m第三题的程序7. gongetidu.m；xldjacobi.m ; xldGS.mx=gongetidu(A,b,s,kmax); x= xldjacobi.(A,b,s,kmax) x=xldGS(A,b,s,kmax)输入值：n*n的矩阵A ，n*1的矩阵b，以及控制迭代停止的s和kmax输出值：用该方法求解问题A*x=b的解x问题分析：要求用共轭梯度法，Jacobi迭代法，G_S迭代法，三种方法求解方程并对比结果。为了方便比较，将三种迭代法的终止迭代条件都设置为残差二范数|r|10(-10)*|b|或者最大允许的迭代次数kmax=2000，三种方法结果如下：jacabi迭代法的解x1b-A*x1G-S迭代法的解x2b-A*x2共轭梯度法的解x3b-A*x31.000075-0.000131.000078-0.000251-3.55E-15-2.000070.00010713-2.000060.000157-2-3.55E-152.9998740.000152732.9998890.0001753-7.11E-15-1.99998-4.22E-05-1.999993.21E-05-2-7.11E-150.999