线性方程组的直接法

Hilbert矩阵Hn=[hi,j]∈Rn*n,其元素hi,j=1/(i+j-1),分别对n=2,3,4…...(1)计算cond(Hn)∞，分析条件数作为n的函数如何让变化？（2）令x=（1，1…..1）T∈Rn，计算bn=Hn*x，然后用Gauss消去法或cholesky方法解方程组bn=Hn*x，解出的解为x,计算剩余变量rn=bn－Hn*x和误差向量；▽x=x－x。（3）分析对每个n，x分量的有效数字如何随n变化，此变化与条件数有何联系。当n为多大时绝对误差达到100%（即X连一位有效数字也没有了）。解：希尔伯特矩阵是著名的病态矩阵。n阶Hilbert矩阵为)12/(1)1/(1/1)1/(13/12/1/12/11nnnnnHn其条件数Cond(Hn)∞如下表所示n23456Cond(Hn)∞27.0000e+0007.4800e+0012.8375e+0049.4366e+0052.9070e+007Input(‘inputn=’);H=hilb(n)Cond(H,inf)随着n的增大，矩阵条件数迅速增加。猜测：希尔伯特矩阵条件数以指数规律增长。即，设矩阵阶数为n，有Cond(Hn)∞≈exp(an+b)用数据拟合的方法验证。问题分析：由于选择拟合函数为指数函数，直接列出超定方程组将是非线性的方程组。为了便于计算，对表中的条件数做对数变换，问题转化为线性拟合问题ln[Cond(Hn)∞]=an+b，(n=2，3，4，5，6)实际操作时使用MATLAB的多项式拟合命令。线性拟合图形如下线性函数的两个系数分别为a=3.3216，b=–3.3474故指数拟合函数为：Cond(Hn)∞≈exp(3.3216n–3.3474)拟合函数的残差向量为r1r2r3r4r50-9.0000e+0018.8000e+0024.0420e+003-5.9917e+005MATLAB程序段如下C=[];fork=2:6H=hilb(k);C=[C,cond(H,inf)];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))2.x=（1，1…..1）T∈Rn，.，令b=Hnx，用高斯消去法求解方程组Hnx=b解出x。。函数定义部分functionx=gauss(A,b)%x=gauss(A,b)n=length(A);a=[A,b];fork=1:n-1maxa=max(abs(a(k:n,k)));ifmaxa==0return;endfori=k:nifabs(a(i,k))==maxay=a(i,k:n+1);a(i,k:n+1)=a(k,k:n+1);a(k,k:n+1)=y;break;endendfori=k+1:nl(i,k)=a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)-l(i,k).*a(k,k+1:n+1);endend%回代ifa(n,n)==0returnendx(n)=a(n,n+1)/a(n,n);fori=n-1:-1:1x(i)=(a(i,n+1)-sum(a(i,i+1:n).*x(i+1:n)))/a(i,i);endclear;n=input('inputn=')；E=ones(n,1);H=hilb(n);b=H*E;x=gauss(H,b)err=norm(b-(x*H)',inf)err=norm(x-E’,inf)err=norm((x-E’)/E,inf)