计算方法大作业(第一次).docVIP

数值计算第一次大作业 实验目的 以Hilbert矩阵为例，研究处理病态问题可能遇到的困难。 内容 Hilbert矩阵的定义是 它是一个对称正定矩阵，而且随着n的增加迅速增加，其逆矩阵，这里 画出之间的曲线（可以用任何的一种范数）。你能猜出之间有何种关系吗？提出你的猜想并想法验证。 用行范数 for n=1:50 for i=1:n for j=1:n A(i,j)=1/(i+j-1); B(i,j)=factorial(n+i-1)*factorial(n+j-1)/((i+j-1)*(factorial(i-1)*factorial(j-1))^2*factorial(n-i)*factorial(n-j)); end end result1=0; for j=1:n result1=result1+A(1,j); end result1=log(result1); result2=0; for i=1:n for j=1:n result2=B(i,j)+result2; end result(i)=log(result2); end m=max(result); x(n)=result1+m; end plot([1:50],x) 对于更大的n值，由于Hilbert逆矩阵中的元素过大，溢出，故在此取50以内的n。 图1 关系曲线图 猜想之间存在线性关系 验证：设 在以上程序基础上，再添加 ; y=x; l=1:40; k=l; p=polyfit(k,y,1) %一次多项式拟合 p = 3.5446 -3.0931 % P=polyfit(k,y,2) %二次多项式拟合 p = -0.0008 3.5778 -3.3253 % P=polyfit(k,y,3) %三次多项式拟合 0.0000 -0.0033 3.6198 -3.4777 % P=polyfit(k,y,4) %四次多项式拟合 -0.0000 0.0002 -0.0082 3.6654 -3.5815 % P=polyfit(k,y,5) %五次多项式拟合 p = 0.0000 -0.0000 0.0007 -0.0156 3.7107 -3.6542 从上式可以看出，高次项系数相对于一次项和常数项系数要小很多， 所以取 2）设是的对角线元素开方构成的矩阵。，不难看出依然是对称矩阵，而且对角线元素都是1。把变成的技术称为预处理。画出之间的曲线（可以用任何一种范数）。你能对于预处理得出什么印象？ 本小题用2范数 clear n=500; c=[]; for k=2:n H=hilb(k); D=diag(sqrt(diag(H))); D1=inv(D); H1=D1*H*D1; c=[c,cond(H1)/cond(H)]; end C=log(c); k=2:n; plot(k,C,r-) 图 2 随n的变化曲线图 从图中给出了函数的变化曲线。我们观察到随着Hilbert矩阵阶数的增大，函数值在[-6,4]区间波动，主要集中在[-3,1]区间。我们知道在时，有，在上图中，我们可以容易观察到，对于大部分，函数值都是小于或者等于零的，这说明经过预处理后的地条件数较小，由于条件数愈大，方程组的病态愈严重，也就愈难得到方程组比较准确地解，所以预处理在一定程度上改善了原Hilbert矩阵的特性。 3）对于，给定不同的右端项。分别用以及，求解，比较计算结果。 取n=4，b=[1;2;3;4] b=[1,2,3,4]; H=hilb(4); H1=inv(H); x1=H1*b x1 = 1.0e+003 * -0.0640 0.9000 -2.5200 1.8200 D=diag(sqrt(diag(H))); D1=inv(D); H2=D1*H*D1; H3=inv(H2); x2=D1*H3*D1*b x2 = 1.0e+003 * -0.0640 0.9000 -2.5200 1.8200 x3=H\b x3 = 1.0e+003 * -0.0640 0.9000 -2.5200 1.8200 同理，当n=5时，b=[1;2;3;4;5] x1 = 1.0e+004 * 0.0125 -0.2880 1.4490 -2.4640 1.3230 x2 = 1.0e+004 * 0.0125 -0.2880 1.4490