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

数值分析课程实验报告

题目：病态线性方程组的求解

理论分析表明，数值求解病态线性方程组很困难。考虑求解如下的线性方程组的求解

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可以看出，随着阶数每递增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的计算结果

从表2可看出，n=2时，四种迭代法都能够收敛，迭代次数最大为e+2量级（J法），最小仅要2次（CG法），并且五种解法都能给出非常精确的结果，最大误差为e-10量级（GS 法）。

表3.n=5的计算结果

从表3可以看出，n=5时，J法已经不收敛，其余解法依然收敛（但值得注意的是GS 法以及SOR法的谱半径以及非常接近1，达到了收敛的边缘），最大迭代次数达到e+6量级（GS法），最小需要7次（CG法），计算精度依然较高，最大误差为e-6量级。SOR法比

GS法节省了较多的迭代次数。

表4.n=10的计算结果

从表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的计算结果

从表5可以看出，n=20时，各种解法的收敛性依然没变（但GS法以及SOR法的谱半径在matlab的十六位显示精度下已经等于1），最大迭代次数e+8量级（SOR法），最小需要10次（CG法），Gauss消元法的计算结果已失去意义（误差e+2量级），迭代法的计算误差均为e-3量级。且此时SOR的迭代次数依然高于GS。

表6.n=50的计算结果

从表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量级），但仔细分析上一节