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

数值分析课程实验报告
题目：病态线性方程组的求解
理论分析表明，数值求解病态线性方程组很困难。

考虑求解如下的线性方程组的求解
Hx = b ，期中H 是Hilbert 矩阵，()ij n n H h ⨯=，1
1
ij h i j =+-，i ，j = 1,2，…，n
1. 估计矩阵的2条件数和阶数的关系
2. 对不同的n ，取(1,1,
,1)n
x =∈
，分别用Gauss 消去，Jacobi 迭代，Gauss-seidel 迭
代，SOR 迭代和共轭梯度法求解，比较结果。

3. 结合计算结果，试讨论病态线性方程组的求解。

解答过程
1.估计矩阵的2-条件数和阶数的关系
矩阵的2-条件数定义为：1
222
()Cond A A A
-=⨯，将Hilbert 矩阵带入有：
1222
()Cond H H H -=⨯
调用自编的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。

C o n d (H )
n
图3.1~10阶Hilbert 矩阵2-条件数的线性拟合
由拟合结果知，Cond （H ）2与n 的关系为：
2lg () 1.65183 1.47863Cond H n
=-+
其线性相关系数r=0.99985，可见二者具有较好的线性关系。

对上式稍作变形得：
1.651831.478632()10n Cond 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，对于迭代法设定终止计算精度为10
10ε-=，所得
计算结果以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量级），但仔细分析上一节
的计算结果不难发现：从n=2到n=5再到n=10，几种迭代法的迭代次数迅速增大，且计算精度显著降低；而从n=10到n=20再到n=50，迭代次数和计算精度却趋于稳定，由于矩阵的条件数一直呈指数增大，因此这种计算结果的趋势明显是不合理的。

究其原因，我认为和一开始图1所示的2-条件数的计算失真是类似的。

当矩阵的条件数非常大时，inv指令的求逆结果不可靠，从而造成图1中2-条件数趋于平稳的那个阶段，而在GS迭代以及SOR迭代中也需要对矩阵进行求逆运算，因此很有可能又是因为inv对高条件数矩阵求逆的偏差造成最终计算结果的失真。

因此，我预计实际上当n=20时，GS法和SOR法的迭代次数以及计算误差要远远高于n=10的值，同理n=50时的结果又要远远大于n=20的值。

不难看出，造成上述两个特点的本质因素还是因为病态矩阵的高条件数，无论收敛性、收敛速度还是计算精度等，都与矩阵的条件数直接相关，因此改善病态矩阵的求解的根本方法还是在于降低矩阵的条件数。

因此我尝试着对Hilbert矩阵进行了一些预处理，并发现在某些预处理方法下条件数的阶数被有效减小，因此预处理是求解病态方程组的一个十分有效且必要的环节。

此外，若给定条件数并横向比较上一节中的不同求解方法，可以看出CG法的表现是最令人满意的，当矩阵阶数n不太大时（不超过10），CG法都能以非常少的迭代次数（只需要几次）以及较高的计算精度完成计算任务，相比之下，GS法以及SOR法的上千万次迭代次数显得非常吃力。

此外，事实上n较小时Gauss消元法的计算结果也比较理想，具有很快的计算速度和较高的计算精度。

我们知道，在传统意义上Gauss消元法的最大弊端在于n非常大时计算量大得令人难以接受（n三次方量级），而迭代法的单步计算量则为n 平方量级，但在求解病态问题时，n只能比较小，且迭代法的迭代次数却非常大，所以造成迭代法的计算速度反而远慢于Gauss消元法。

综上，求解病态方程组时首先应该对系数矩阵进行预处理，尽可能降低它的条件数；其次可以选择共轭梯度（CG）法对其进行求解。

如果阶数n不大，Gauss消元法也是比较有效的。