数值分析课程实验报告数值分析课程实验报告题目：病态线性方程组的求解题目：病态线性方程组的求解理论分析表明，数值求解病态线性方程组很困难考虑求解如下的线性方程组的求解Hx = b，期中 H 是 Hilbert 矩阵，，，i，j = 1,2，…，n()ijn nHh1 1ijhij1.估计矩阵的 2 条件数和阶数的关系2.对不同的 n，取，分别用 Gauss 消去，Jacobi 迭代，Gauss-seidel 迭(1,1,,1)nx K¡代，SOR 迭代和共轭梯度法求解，比较结果 3.结合计算结果，试讨论病态线性方程组的求解解答过程解答过程1.估计矩阵的估计矩阵的 2-条件数和阶数的关系条件数和阶数的关系矩阵的 2-条件数定义为：，将 Hilbert 矩阵带入有：1 222( )Cond AAA1 222()Cond HHH调用自编的 Hilbert_Cond 函数对其进行计算，取阶数 n = 50，可得从 1 阶到 50 阶的 2-条件 数，以五位有效数字输出，其中前 10 项见表 1表 1.前十阶 Hilbert 矩阵的 2-条件数阶数123452-条件数119.281524.061.5514e+0044.7661e+005阶数6789102-条件数1.4951e+0074.7537e+0081.5258e+0104.9315e+0111.6025e+013从表 1 可以看出，随着阶数每递增 1，Hilbert 矩阵的 2-条件数都至少增加一个数量级， 但难以观察出明显的相依规律。
故考虑将这些数据点绘制在以 n 为横轴、Cond（H）2 为 纵轴的对数坐标系中（编程用 Hilbert_Cond 函数同时完成了这个功能） ，生成结果如图 1 图 1.不同阶数下 Hilbert 矩阵的 2-条件数分布由图可见，当维数较小时，在 y-对数坐标系中 Cond（H）2 与 n 有良好的线性关系； 但 n 超过 10 后，线性趋势开始波动，n 超过 14 后更是几乎一直趋于平稳事实上，从 n = 12 开始，系统便已经开始提出警告：“Warning: Matrix is close to singular or badly scaled.Results may be inaccurate.”也就是说，当 n 较大时，H 矩阵已经接近奇异，计算结 果可能是不准确的通过查阅相关资料，我找到了造成这种现象的原因：在 matlab 中，用 inv 函数求条件数过大的矩阵的逆矩阵将是不可靠的而调用系统自带的专门对 Hilbert 矩 阵求逆的 invhilb（n）函数则不存在这个问题，生成结果如图 2图 2. 修正后的不同阶数下 Hilbert 矩阵的 2-条件数分布简便起见，取 n 不大于 10 的前十项进行线性拟合，结果如图 3。
0123456789101111010010001000010000010000001E71E81E91E101E111E121E13Cond(H)Linear Fit of Sheet1 Cond(H)Cond(H)nEquationy = a + b*xWeightNo WeightingResidual Sum of Squares0.05264Pearson's r0.99985Adj. R-Square0.99967ValueStandard ErroCond(H)Intercept-1.651830.05541Cond(H)Slope1.478630.00893图 3.1~10 阶 Hilbert 矩阵 2-条件数的线性拟合由拟合结果知，Cond（H）2 与 n 的关系为：2lg()1.65183 1.47863Cond Hn 其线性相关系数 r=0.99985，可见二者具有较好的线性关系对上式稍作变形得：1.65183 1.47863 2()10nCond H这就是对 Hilbert 矩阵的-2 条件数与其阶数 n 的关系估计可见 Hilbert 矩阵的 2-条件数会 随其阶数 n 的增加呈指数增大趋势，因此当 n 较大时 Hilbert 矩阵将是严重病态的，甚至导 致 matlab 中 inv 求逆运算失真。
2.对不同的对不同的 n，采用各种方法求解方程，采用各种方法求解方程调用自编的 Calc 主函数（其中包括的 Hilbert 函数以及 b_函数可创建出对应阶数的 H 矩阵以及向量 b，Gauss_Cal 函数、Jacobi_Cal 函数、Gauss_Seidel_Cal 函数、SOR_Cal 函 数（该函数自动寻找最优松弛因子，然后以最优因子进行求解）以及 CG_Cal 函数则可完成各自方法的求解） ，分别取 n = 2,5,10,20,50，对于迭代法设定终止计算精度为，1010所得计算结果以 16 位有效数字输出，分别见表 2~表 6表 2.n=2 的计算结果求解方法求解方法迭代矩阵谱半径迭代矩阵谱半径 p；；是否收敛是否收敛迭代迭代次数次数N解向量解向量 x误差误差 eGAUSS消元消元——10.9999999999999998.00593208497344e-016Jacobi 迭迭代代0.866025403784439收敛1691.000000000016 1.0000000000485.05958147990869e-011GS 迭代迭代0.75收敛771.00000000015982 0.9999999997602732.88116190608079e-010SOR 迭迭代代0.35（）1.35b收敛261.00000000001783 0.9999999999822262.51745780938445e-011CG 迭代迭代—收敛2113.27844930319752e-015从表 2 可看出，n=2 时，四种迭代法都能够收敛，迭代次数最大为 e+2 量级（J 法） ， 最小仅要 2 次（CG 法） ，并且五种解法都能给出非常精确的结果，最大误差为 e-10 量级 （GS 法） 。
表 3.n=5 的计算结果求解方法求解方法迭代矩阵谱半径迭代矩阵谱半径 p；；是否收敛是否收敛迭代次数迭代次数 N误差误差 eGAUSS 消元消元——3.47269899767798e-012Jacobi 迭代迭代3.44414219116596不收敛——GS 迭代迭代0.999957671222957收敛2318152.3623228567773e-006SOR 迭代迭代0.999190149180609（）1.83b收敛215008.86896685436631e-008CG 迭代迭代—收敛75.84409588812087e-011从表 3 可以看出，n=5 时，J 法已经不收敛，其余解法依然收敛（但值得注意的是 GS 法以及 SOR 法的谱半径以及非常接近 1，达到了收敛的边缘） ，最大迭代次数达到 e+6 量 级（GS 法） ，最小需要 7 次（CG 法） ，计算精度依然较高，最大误差为 e-6 量级SOR 法比 GS 法节省了较多的迭代次数表 4.n=10 的计算结果求解方法求解方法迭代矩阵谱半径迭代矩阵谱半径 p；；是否收敛是否收敛迭代次数迭代次数 N误差误差 eGAUSS 消元消元——0.00070763269657929Jacobi 迭代迭代7.77981513192998不收敛——GS 迭代迭代0.999999999997045收敛30868430.00213101923167426SOR 迭代迭代0.999999999871931（）1.93b收敛103038010.00177481054386451CG 迭代迭代—收敛80.00116732746950012从表 4 可以看出，n=10 时，各种解法的收敛性与 n=5 时相同（GS 法与 SOR 法的谱半 径进一步趋近 1） ，最大迭代次数达 e+8 量级（SOR 法） ，最小需要 8 次（CG 法） ，计算精 度已经较低，最大误差达到 e-3 量级。
此时有一个异常现象：SOR 法的迭代次数不但不比 GS 法少，反而超出好几倍但根据计算出的两种算法下的迭代矩阵谱半径，可以看出 SOR 法为 0.999999999871931，而 GS 法为 0.999999999997045，在最优松弛因子之下的 SOR 法确实具有更小的迭代矩阵谱半径，因此应当更快收敛经检查，计算程序的编制应 该没有错误，但计算结果确实与理论不符，希望老师指点迷津！表 5.n=20 的计算结果求解方法求解方法迭代矩阵谱半径迭代矩阵谱半径 p；；是否收敛是否收敛迭代次数迭代次数 N误差误差 eGAUSS 消元消元——105.281872512979Jacobi 迭代迭代16.4920989837926不收敛——GS 迭代迭代1收敛125450460.00166739398859172SOR 迭代迭代1（）0.51b收敛131181710.00111174295388826CG 迭代迭代—收敛100.00220319859317489从表 5 可以看出，n=20 时，各种解法的收敛性依然没变（但 GS 法以及 SOR 法的谱半 径在 matlab 的十六位显示精度下已经等于 1） ，最大迭代次数 e+8 量级（SOR 法） ，最小需 要 10 次（CG 法） ，Gauss 消元法的计算结果已失去意义（误差 e+2 量级） ，迭代法的计算 误差均为 e-3 量级。
且此时 SOR 的迭代次数依然高于 GS表 6.n=50 的计算结果求解方法求解方法迭代矩阵谱半径迭代矩阵谱半径 p；；是否收敛是否收敛迭代次数迭代次数 N误差误差 eGAUSS 消元消元——1697.4503939611Jacobi 迭代迭代42.6768950976645不收敛——GS 迭代迭代1收敛181452240.00294769067781031SOR 迭代迭代1（）0.06b收敛39574210.00150790137625883CG 迭代迭代—收敛150.00187865108201776从表 6 可以看出，n=50 时，计算结果的特点与 n=20 相似，不同之处在于 Gauss 消元 法的误差进一步增大（e+3 量级） 3. 综合讨论综合讨论根据上一部分的计算结果，可以总结出求解病态矩阵的两个特点：1.收敛性差，收敛速度慢收敛性差，收敛速度慢 对于本例，当阶数 n 大于 2 以后 J 法就不收敛了，这是收敛性差的一个体现GS 法和 SOR 法虽然一直保持收敛（事实上这是由于 Hilbert 矩阵是正定对称的，所以决定了 GS 法 以及 SOR 法一定收敛） ，但随着 n 的增大它们的谱半径迅速趋近于 1（n=20 时 matlab 的 16 位显示精度已经无法显示出它们的谱半径与 1 的差别） ，根据理论我们知道，趋近于 1 的谱 半径决定了计算的收敛速度将会非常慢，之前计算结果中的迭代次数 N 也很好地验证了这 一点。
2.计算精度低，阶数较高时误差惊人计算精度低，阶数较高时误差惊人 根据理论，线性方程组系数矩阵的条件数的直观意义就是初始数据扰动的相对误差的 放大倍数由于数值计算时的舍入误差等因素，初始数据的微小扰动难以避免，因此解病 态方程组时，这种微小的初始相对误差将被病态矩阵巨大的条件数所放大，最后造成计算 结果的巨大偏差对于本例，n=20 时 Gauss 消元法的计算误差竟然达到 e+2 量级，使得计 算结果完全失去意义，而其他几种迭代法的计算误差看似能够接受（e-3 量级） ，但仔细分析上一节的计算结果不难发现：从 n=2 到 n=5 再到 n=10，几种迭代法的迭代次数迅速增大， 且计算精度显著降低；而从 n=10 到 n=20 再到 n=50，迭代次数和计算精度却趋于稳定，由 于矩阵的条件。