hilb病态矩阵.doc

许多线性代数方程组理论上是可以理解的，但是实际计算中由于受到舍入误差的影响无法得到精确解。此类问题成为病态问题。 一、希尔伯特矩阵的来源及矩阵的正定性 设f(x)∈C[0，1]，求n次多项式 P(x) = a0 + a1x + a2 x2 + …… + an x n 使得 取最小值。其中，P(x)称为函数f(x)的最佳平方逼近多项式。由极值必要条件，对 中的变量求导数，得 令其为零，得方程组 令，将方程组写为矩阵形式 方程组的系数矩阵就是著名的Hilbert矩阵。 多项式可表示为内积形式 所以 积分，得 由此可知，希尔伯特矩阵是对称正定矩阵。 二、希尔伯特矩阵的病态性 希尔伯特矩阵是著名的病态矩阵。n阶Hilbert矩阵为 其条件数 Cond(Hn) 如下表所示 n 2 3 4 5 6 Cond(Hn) 1.9281e+001 5.2406e+002 1.5514e+004 4.7661e+005 1.4951e+007 随着n的增大，矩阵条件数迅速增加。 猜测：希尔伯特矩阵条件数以指数规律增长。即，设矩阵阶数为n，有 Cond(Hn) ≈ exp( a n + b ) 用数据拟合的方法验证。 问题分析：由于选择拟合函数为指数函数，直接列出超定方程组将是非线性的方程组。为了便于计算，对表中的条件数做对数变换，问题转化为线性拟合问题 ln[Cond(Hn)] = a n + b，( n = 2，3，4，5，6) 实际操作时使用MATLAB的多项式拟合命令。线性拟合图形如下 线性函数的两个系数分别为 a = 3.3935，b = – 3.8811 故指数拟合函数为： Cond(Hn) ≈ exp(3.3935 n –3.8811 ) 拟合函数的残差向量为 r1 r2 r3 r4 r5 9.9860e-001 -2.0231e+001 -6.8994e+002 -5.7825e+003 5.9012e+005 MATLAB程序段如下 C=[]; for k=2:6 H=hilb(k);%k阶希尔伯特矩阵 C=[C,cond(H)]; end LC=log(C);n=2:6; P=polyfit(n,LC,1) plot(n,LC,o,n,polyval(P,n)) residure=C-exp(polyval(P,n)) 三、希尔伯特矩阵方程组求解实验 实验1：CG法实验 设E为n阶全“1”列向量，令b=HnE，用共轭梯度法求解方程组 Hnx = b 并与MATLAB的反斜杠命令求解方程组的结果作对比。 N 迭代次数 残差平方和 误差 MATLAB直接解误差 4 4 1.2094e-008 6.3891e-009 7.1987e-013 6 4 1.1845e-007 0.0055 3.7577e-010 8 6 4.4100e-009 0.0018 9.3203e-007 10 6 3.2342e-008 0.0034 1.4717e-004 12 5 5.0641e-007 0.0056 0.2385 14 5 3.3218e-007 0.0079 14.3112 16 5 7.0032e-007 0.0102 8.2762 实验程序： n=input(input n=); E=ones(n,1); H=hilb(n);b=H*E; [x,Res,k]=cg(H,b); k,Res x1=H\b; error=norm(x-E,inf) error1=norm(x1-E,inf) function [sol,res,k] = CG (A,b) %CG method 2009_10-27 sol = zeros(size(b)); r = b - A * sol; p = r ; ro1 = r * r; k = 0; while (ro1 1.0e-12) ro = ro1; Ap = A * p; alpha = ro / ( Ap * p ) ; sol = sol + alpha * p ; r = r - alpha * Ap; ro1= r*r; beta = ro1 / ro ; p = r + beta* p; k=k+1; end res=sqrt(ro1); 实验2：Seidel 迭代法 N 迭代次数 相邻两次差 误差 MATLAB直接解误差 4 9721 9.9939e-009 1.0291e-005 7.1987e-013 6 115045 1.0000e-008 0.0059 3.7577e-010 8 212504 9.9999e-009 0.0031 9.3203e-007 10 371557 1.0000e-008 0.