计算方法大作业第一次

数值计算第一次大作业

实验目的 以Hilbert 矩阵为例，研究处理病态问题可能遇到的困难。

内容 Hilbert 矩阵的定义是

,()

11/21/31/1/21/31/41/(1)1/31/41/51/(2)1/1/(1)1/(2)1/(21)n i j H h n n n n n n n =⎡⎤⎢⎥+⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥++-⎣⎦

它是一个对称正定矩阵，而且()n cond H 随着n 的增加迅速增加，其逆矩阵1,()n i j H α-=，

这里

,2(1)(1)!(1)!(1)[(1)!(1)!]()!()!

i j i j n i n j i j i j n i n j α+-+-+-=+----- 1) 画出ln(())~n cond H n 之间的曲线（可以用任何的一种范数）。你能猜出

ln(())~n cond H 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)*facto rial(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 ln(())~n cond H n 关系曲线图

猜想ln(())~n cond H n 之间存在线性关系

验证：设ln(()n cond H an b ∞=+

在以上程序基础上，再添加

>>;>> 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

从上式可以看出，高次项系数相对于一次项和常数项系数要小很多，

所以取ln(() 3.5446 3.0931n cond H n ∞=-

2）设D 是n H 的对角线元素开方构成的矩阵。11

ˆn n H D H D --=，不难看出ˆn H 依然是对称矩阵，而且对角线元素都是1。把n H 变成ˆn

H 的技术称为预处理。画出ˆln(()/())n n

cond H cond H n 之间的曲线（可以用任何一种范数）。你能对于预处理得出什么印象？

本小题用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 ˆln(()/())n n

cond H cond H 随n 的变化曲线图 从图中给出了函数ˆln(()/())n n

cond H cond H 的变化曲线。我们观察到随着Hilbert 矩阵阶数的增大，函数值在[-6,4]区间波动，主要集中在[-3,1]区间。我们知

道在ˆ()()n n cond H cond H ≤时，有ˆln(()/())0n n

cond H cond H ≤，在上图中，我们可以容易观察到，对于大部分n ，函数值都是小于或者等于零的，这说明n H 经过预处理

后的ˆn

H 地条件数较小，由于条件数愈大，方程组的病态愈严重，也就愈难得到方程组比较准确地解，所以预处理在一定程度上改善了原Hilbert 矩阵的特性。

3）对于412n ≤≤，给定不同的右端项b 。分别用111112ˆˆ,,n n x H b x H D b x D x

----===以及3\n x H b =，求解n H x b =，比较计算结果。 取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