hilb病态矩阵

《hilb病态矩阵》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

许多线性代数方程组理论上是可以理解的，但是实际计算中由于受到舍入误差的影响无法得到精确解。此类问题成为病态问题。一、希尔伯特矩阵的来源及矩阵的正定性设f(x)∈C[0，1]，求n次多项式P(x)=a0+a1x+a2x2+……+anxn使得取最小值。其中，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(an+b)用数据拟合的方法验证。问题分析：由于选择拟合函数为指数函数，直接列出超定方程组将是非线性的方程组。为了便于计算，对表中的条件数做对数变换，问题转化为线性拟合问题ln[Cond(Hn)]=an+b，(n=2，3，4，5，6)实际操作时使用MATLAB的多项式拟合命令。线性拟合图形如下线性函数的两个系数分别为a=3.3935，b=–3.8811故指数拟合函数为：Cond(Hn)≈exp(3.3935n–3.8811)拟合函数的残差向量为r1r2r3r4r59.9860e-001-2.0231e+001-6.8994e+002-5.7825e+0035.9012e+005MATLAB程序段如下C=[];fork=2:6H=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('inputn=');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)%CGmethod2009_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-013 61150451.0000e-0080.00593.7577e-01082125049.9999e-0090.00319.3203e-007103715571.0000e-0080.00421.4717e-004实验程序：?typehilblab00n=input('inputn=');E=ones(n,1);H=hilb(n);b=H*E;x=zeros(n,1);err=1;k=0;whileerr>1.0e-8k=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的对角线元素开方构成的对角矩阵，令，不难看出该矩阵仍然是对称正定矩阵，而且其对角元素全为1。将Hn变换为的技术称为预处理（Preconditioning）。试分析函数随n变化的规律，绘制该函数的曲线。三、如果只用左预处理，即D取希尔伯特矩阵对角元构成对角矩阵，用D的逆阵左乘H后条件数如何。