计算方法Hilbert矩阵病态问题

Solution （1）
猜想：()()
ln cond Hn n 成线性比例关系。

验证：猜想部分正确。

根据不同n 建立对应的Hilbert 矩阵，计算该矩阵的2-范数的条件数，并绘制曲线。

编程实现见附录1中solution(1)。

实际编程中，分别取N=5，=10N ，=20N ，N=30，可以发现()()
ln cond Hn n 曲线先是按线性规律变化；当n 达到大约12以后，()()
ln cond Hn 的取值在40~50之间产生振荡。

各结果如下：
N=5时 N=10时
ln(cond(Hn))-n
n
l n (c o n d (H n ))
ln(cond(Hn))-n
n
l n (c o n d (H n ))
ln(cond(Hn))-n
n
l n (c o n d (H n )
)
ln(cond(Hn))-n
n
l n (c o n d (H n ))
N=20时N=30时
Solution （2）
猜想：通过对Hilbert 矩阵的预处理，能改善Hilbert 的病态性质。

验证：根据不同n 建立对应的Hilbert 矩阵，计算该矩阵的2-范数的条件数，其次计算经过预处理后的H^矩阵的2-范数的条件数，并绘制()()()
ln ^/cond Hn cond Hn n 曲线。

编程实现见附录1中的solution(2)。

实际编程中，分别取N=5，=12N ，=50N ，N=100，可以发现所得到的变化曲线
()()()ln ^/cond Hn cond Hn n 在12n ≤时，()()()ln ^/cond Hn cond Hn 的取值较为
平缓的减小，且均小于0。

而随着Hilbert 矩阵阶数的增大，()()()
ln ^/cond Hn cond Hn 的取值在[]-7,3区间中振荡，主要集中在[]-31，，且多数取值()()()
ln ^/0cond Hn cond Hn ≤。

则对多数的Hilbert 矩阵，其条件数在经过预处理后，有()()^cond Hn cond Hn ≤，即预处理在一定程度上改善了原Hilbert 矩阵的病态性质。

各结果如下：
N=5时 N=12时
ln(cond(Hn)/cond(Hn))-n
n
l n (c o n d (H n )/c o n d (H n ))
ln(cond(Hn)/cond(Hn))-n
n
l n (c o n d (H n )/c o n d (H n ))
N=50时N=100时
ln(cond(Hn)/cond(Hn))-n
n
l n (c o n d (H n )/c o n d (H n ))
ln(cond(Hn)/cond(Hn))-n
n
l n (c o n d (H n )/c o n d (H n ))
取定=8n ，取不同的右端项b ：b1=ones(N,1);b2=[1 1 1 1 1 1 1 0.999999]';b3=[1 1 1 1 1 1 1 1.000001]';分别按照题示3种方法求解方程解。

编程实现见附录1中的solution(3)。

对各结果列表比较如下：
分析比较各结果，对方法一-11n
x =H b ，方法二^1
--1
^
n x =H D b ，^
1
2-x =D x ，方法三
3\n x =H b 有：
（1） 对相同的右端项b ，方法一、三计算的结果较为接近，但也有所差别。

（2） 对不同右端项，微小扰动下，方法一、二、三的解均有所变化，方法一的变化
最明显，方法三变化较小，而方法二的解基本保持稳定。

综上，对=8n 时，说明了Hilbert 矩阵的病态性质：(i) Hilbert 矩阵求逆不精确，从而利用方法一较难得到方程组的精确解；(ii)以Hilbert 矩阵为系数矩阵的线性方程组的解对数据误差有较大的敏感性，方法一、三的结果均有较明显的差距。

而经过预处理后，Hilbert 矩阵的病态性质有所改善，利用方法二计算得到的解基本保持稳定。

取定右端项为Hilbert 矩阵的第一列，=(:,1)b H ，对不同的=5n ，=8n ，=12n ，=15n ，用高斯-塞德尔迭代法求解方程解。

编程实现见附录1中的solution(4)；其中利用高斯-塞德尔迭代法的函数，由function [n,x]=gaussseidel(A,b,X,delta,max)给出，见附录2。

对各结果列表比较如下：
表2 不同右端项下3种方法求解结果比较
无直接联系。

附录1 Solution（1）
% solution1
clear all
N=30; % N分别取5,10,20，30
for i=1:N
H=hilb(i); % 定义Hilebert矩阵
h=cond(H); % 求对应的2-范数条件数
y(i)=log(h); % 求ln(cond(Hn))
end
x=1:N;
plot(x,y);
title('ln(cond(Hn))-n');
xlabel('n'),ylabel('ln(cond(Hn))')
Solution（2）
% solution 2
clear all
N=30;% N分别取5,10,20，30
for i=1:N
H=hilb(i);% 定义Hilebert矩阵
h=cond(H);% 求对应的2-范数条件数
D=diag(sqrt(diag(H)));% 求D矩阵
Dinv=inv(D);% 求D的逆矩阵
HH=Dinv*H*Dinv;% 求H^矩阵
hh=cond(HH);% 求H^矩阵对应的2-范数条件数
y(i)=log(hh/h);% 求ln(cond(Hn^)/cond(Hn)) end
x=1:N;
plot(x,y);
title('ln(cond(Hn)/cond(Hn))-n');
xlabel('n'),ylabel('ln(cond(Hn^)/cond(Hn))')
% solution 2
clear all
N=8;
H=hilb(N);
Hinv=invhilb(N);
D=diag(sqrt(diag(H)));
Dinv=inv(D);
HH=Dinv*H*Dinv;
b1=ones(N,1);
b2=[1 1 1 1 1 1 1 0.999999]'; b3=[1 1 1 1 1 1 1 1.000001]';
x11=Hinv*b1
x1p=HH*Dinv*b1;
x12=Dinv*x1p
x13=H\b1
x21=Hinv*b2
x2p=HH*Dinv*b2;
x22=Dinv*x2p
x23=H\b2
x31=Hinv*b3
x3p=HH*Dinv*b3;
x32=Dinv*x3p
x33=H\b3
% solution 4
clear all
N=15;
H=hilb(N);
b=H(:,1);
P=ones(N,1);
delta=0.000001;
max=1000000;
x=gaussseidel(H,b,P,delta,max);
附录2
高斯-塞德尔迭代法
function [n,x]=gaussseidel(A,b,X,delta,max)
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，% max为最大迭代次数，delta为误差精度
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1;
m=length(A);
I=eye(m);
D=diag(diag(A));
L=tril(-A)+D;
U=triu(-A)+D;
G=inv(D-L)*U; %计算迭代矩阵
f=inv(D-L)*b; %计算迭代矩阵的常数项
while n<=max
x=G*X+f; %迭代格式
if norm(x-X,2)<delta
%若达到精度要求就结束程序，输出迭代次数和方程组的解
disp('迭代次数为');n
disp('方程组的解为');x
return;
end
X=x;n=n+1;
end
%下面：如果达到迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）disp('在最大迭代次数内不收敛!');
disp('最大迭代次数后的结果为');x
end。