数值分析实验报告——Hilbert矩阵的求解.docVIP

数值分析课程实验报告 题目：病态线性方程组的求解 理论分析表明，数值求解病态线性方程组很困难。考虑求解如下的线性方程组的求解Hx = b，期中H是Hilbert矩阵，，，i，j = 1,2，…，n 1. 估计矩阵的2条件数和阶数的关系 2. 对不同的n，取，分别用Gauss消去，Jacobi迭代，Gauss-seidel迭代，SOR迭代和共轭梯度法求解，比较结果。 3. 结合计算结果，试讨论病态线性方程组的求解。 解答过程 1.估计矩阵的2-条件数和阶数的关系 矩阵的2-条件数定义为：，将Hilbert矩阵带入有： 调用自编的Hilbert_Cond函数对其进行计算，取阶数n = 50，可得从1阶到50阶的2-条件数，以五位有效数字输出，其中前10项见表1。 表1.前十阶Hilbert矩阵的2-条件数 阶数 1 2 3 4 5 2-条件数 1 19.281 524.06 1.5514e+004 4.7661e+005 阶数 6 7 8 9 10 2-条件数 1.4951e+007 4.7537e+008 1.5258e+010 4.9315e+011 1.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。 图3.1~10阶Hilbert矩阵2-条件数的线性拟合 由拟合结果知，Cond（H）2与n的关系为： 其线性相关系数r=0.99985，可见二者具有较好的线性关系。对上式稍作变形得： 这就是对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，对于迭代法设定终止计算精度为，所得计算结果以16位有效数字输出，分别见表2~表6。  表2.n=2的计算结果 求解方法 迭代矩阵谱半径p； 是否收敛 迭代次数N 解向量x 误差e GAUSS消元 — — 1 0.999999999999999 8.00593208497344e-016 Jacobi迭代 0.866025403784439 收敛 169 1.000000000016 1.000000000048 5.05958147990869e-011 GS迭代 0.75 收敛 77 1.00000000015982 0.999999999760273 2.88116190608079e-010 SOR迭代 0.35 （） 收敛 26 1.00000000001783 0.999999999982226 2.51745780938445e-011 CG迭代 — 收敛 2 1 1 3.27844930319752e-015 从表2可看出，n=2时，四种迭代法都能够收敛，迭代次数最大为e+2量级（J法），最小仅要2次（CG法），并且五种解法都能给出非常精确的结果，最大误差为e-10量级（GS法）。 表3.n=5的计算结果 求解方法 迭代矩阵谱半径p； 是否收