hilb病态矩阵.doc

许多线性代数方程组理论上是可以理解的，但是实际计算中由于受到舍入误差的影响无法得到精确解。此类问题成为病态问题。一、希尔伯特矩阵的来源及矩阵的正定性设f(x)C0，1，求n次多项式P(x) = a0 + a1x + a2 x2 + + an x n使得取最小值。其中，P(x)称为函数f(x)的最佳平方逼近多项式。由极值必要条件，对中的变量求导数，得令其为零，得方程组令，将方程组写为矩阵形式方程组的系数矩阵就是著名的Hilbert矩阵。多项式可表示为内积形式所以积分，得由此可知，希尔伯特矩阵是对称正定矩阵。二、希尔伯特矩阵的病态性希尔伯特矩阵是著名的病态矩阵。n阶Hilbert矩阵为其条件数 Cond(Hn) 如下表所示n23456Cond(Hn)1.9281e+0015.2406e+0021.5514e+0044.7661e+0051.4951e+007随着n的增大，矩阵条件数迅速增加。猜测：希尔伯特矩阵条件数以指数规律增长。即，设矩阵阶数为n，有Cond(Hn) exp( a n + b )用数据拟合的方法验证。问题分析：由于选择拟合函数为指数函数，直接列出超定方程组将是非线性的方程组。为了便于计算，对表中的条件数做对数变换，问题转化为线性拟合问题lnCond(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 )拟合函数的残差向量为r1r2r3r4r59.9860e-001-2.0231e+001-6.8994e+002-5.7825e+0035.9012e+005MATLAB程序段如下C=;for k=2:6 H=hilb(k);%k阶希尔伯特矩阵 C=C,cond(H);endLC=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直接解误差441.2094e-0086.3891e-0097.1987e-013641.1845e-0070.00553.7577e-010864.4100e-0090.00189.3203e-0071063.2342e-0080.00341.4717e-0041255.0641e-0070.00560.23851453.3218e-0070.007914.31121657.0032e-0070.01028.2762实验程序：n=input(input n=);E=ones(n,1);H=hilb(n);b=H*E;x,Res,k=cg(H,b);k,Resx1=Hb;error=norm(x-E,inf)error1=norm(x1-E,inf)function sol,res,k = CG (A,b) %CG method 2009_10-27sol = 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;endres=sqrt(ro1);实验2：Seidel 迭代法N迭代次数相邻两次差误差MATLAB直接解误差497219.9939e-0091.0291e-0057.1987e-01361150451.0000e-0080.00593.7577e-01082125049.9999e-0090.00319.3203e-007103715571.0000e-0080.00421.4717e-004实验程序：?type hilblab00n=input(input n=);E=ones(n,1);H=hilb(n);b=H*E;x=zeros(n,1);err=1;k=0;while err1.0e-8 k=k+1;L=tril(H); xx=x+inv(L)*(b-H*x); err=norm(xx-x,inf); x=xx;endk,errerror=norm(x-E,inf)x1=Hb;error1=norm(x1-E,inf)实验问题：一、 由于线性拟合所得残差向量较大，用实验确定：选择用几次多项式拟合较好。二、 设D是由Hn的对角线元素开方构成的对角矩阵，令，不难看出该矩阵仍然是对称正定矩阵，而