Hilbert代数方程组的数值解法.docx
《Hilbert代数方程组的数值解法.docx》由会员分享，可在线阅读，更多相关《Hilbert代数方程组的数值解法.docx（16页珍藏版）》请在冰豆网上搜索。
Hilbert代数方程组的数值解法
Hilbert代数方程组的数值解法
1Hilbert矩阵的条件数和矩阵的阶数的关系
编制Matlab程序：
clear,clc
formatlong
%formatshort
forn=1:
14;
A=hilb（n）;
%A=rand（n）;
condA（n）=cond（A,inf）;
end
plot（log10（condA））
得出图形
图1：
Hilbert矩阵和阶次的关系
由图可见a.Hilbert矩阵的2-条件数的对数与阶次近似正比;
b.阶次超过12后，计算困难，舍入误差会导致计算不准
结论：
Hilbert矩阵的2-条件数随阶次指数增长，关系大概是：
Log10（Cond（A））=1.33791720780254（n-1）
2.各种求解方法的对比
2.1直接法：
Gauss消去方法（程序见附录，下同）
在双精度型变量下，即eps=2.220446049250313e-016，对于直接法Gauss消去方法，求解Hilbert矩阵方程组，在阶次n=13时，解得:
X=[10.999971.00130.978811.1906-0.0354414.6171-7.396714.089-12.5429.9162-2.38171.5624]
出现明显错误，可见Gauss消去方法对病态问题比较敏感。
求解阶次不高。
2.2Jacobi迭代法
很遗憾，jakobi迭代法的迭代矩阵的谱半径随阶次线性上升（程序见附录2），普遍大于1，计算结果发送，无法迭代出满意的结果。
需放弃
图2：
Jacobi迭代矩阵的谱半径
2.3Gauss-Seidel迭代与SOR迭代求解
先分析SOR迭代矩阵的谱半径和松弛因子w的关系
图3：
SOR迭代矩阵的谱半径和松弛因子w的关系
可以发现SOR迭代和G-S迭代的谱半径都普遍接近于1，因此松弛因子的选择对计算结果的影响不大。
我们选择w=1，即G-S迭代来做数值实验：
选定e=0.000001；
迭代次数：
15610次
下降曲线：
图4：
G-S迭代下降曲线
取w=1.5时，
迭代次数需58798次
前100步下降曲线
图5：
SOR（w=1.5）下降曲线
w=0.5迭代次数16290
图6：
SOR（w=0.5）下降曲线
2.4SOR按照寻优算法的迭代求解
寻优得出w=0.1时，谱半径相对较小，迭代次数只需3689次，下降曲线也非常快，如图7
图7：
SOR（w=0.1）下降曲线
2.5共轭梯度法迭代求解
奇迹发生了，迭代只需要9次就完成了，下降也非常快，见图
图8：
共轭梯度法下降曲线
2.6各算法综合对比
综上所述，在双精度的数值精度下和e=0.00001下，在笔记本可接受的计算时间里，做各种算法的对比如下：
Gauss
Jacobi
G-S
SOR
CG
能计算到n阶
12
2
1200
2000
5500
迭代次数
/
/
最多
中等
较少
计算时间
/
慢
较快
超快
稳定性
较差
最差
一般
较好
很好
3.讨论病态方程组的求解
3.1对于Hilbert病态方程组的求解
Gauss直接法和Jacobi迭代法都比较差，G-S法和SOR虽能求解，但由于其条件数都接近于1，收敛太慢，计算时间太长，因此实际可计算的阶次不高。
3.2采用高精度计算
我们分别采用单精度和双精度型变量，以Gauss消去方法为求解办法，做数值实验，
x2=gauss（A,b）
x3=gauss（single（A）,single（b））
当阶次n=6时，得
x2=
0.99999999999908
1.00000000002630
0.99999999982172
1.00000000046435
0.99999999948691
1.00000000020235
x3=
1.0003130
0.9906988
1.0645127
0.8293791
1.1905806
0.9242350
可以发现，x3已经明显偏离了正解，得出结论，采用高精度计算所得解有更多有效位。
3.3采用预处理方法
预处理共轭梯度法见附录，采用后，可以发现，计算4000阶的问题，速度加快很多，迭代次数也减少了。
附录1Gauss消去方法
functionx=Gauss（a,b）
%Gauss消去法
[n,m]=size（a）;
nb=length（b）;
det=1;%存储行列式值
x=zeros（n,1）;
fork=1:
n-1
fori=k+1:
n
ifa（k,k）==0
return
end
m=a（i,k）/a（k,k）;
forj=k+1:
n
a（i,j）=a（i,j）-m*a（k,j）;
end
b（i）=b（i）-m*b（k）;
end
det=det*a（k,k）;
end
det=det*a（n,n）;
fork=n:
-1:
1%回代
forj=k+1:
n
b（k）=b（k）-a（k,j）*x（j）;
end
x（k）=b（k）/a（k,k）;
end
附录2Jacobi迭代的收敛性分析
clear,clc
formatlong
%formatshort
forn=1:
25;
A=hilb（n）;
D=diag（diag（A））;
B=inv（D）*（D-A）;
x=ones（n,1）;
e=0.000001;
e=eps;
b=A*x;
f=inv（D）*b;
%%谱半径大于等于1就不收敛
p（n）=max（abs（eig（B）））;
end
plot（p）
xlabel（'阶次'）
ylabel（'谱半径'）
title（'Jacobi迭代收敛性分析'）
附录3G-S迭代
function[x,k,Ve]=gauss_seidel（A,b,e,M）
[NAr,NA]=size（A）;
ifNAr~=NA,error（'Aisnotasquarematrix'）;end
%检查A，b的大小是否匹配
[Nbr,Nb]=size（b）;
if（（Nbr~=NA）&&（1~=Nb））,error（'Aandbhavenon-compatibledimensions'）;end
%e不能小于计算机的最小精度
ifex=zeros（size（b））;%初始解设置为与b同型的零向量
k=0;%迭代次数的记数变量，初始量设为0
r=1+e;%前后项之间的误差
%构造D、B、f矩阵
D=diag（diag（A））;
U=-triu（A,1）;
L=-tril（A,-1）;
DsubL=（D-L）;
B=inv（DsubL）*U;
f=inv（DsubL）*b;
%为了画图，返回误差向量
Ve=zeros（100,1）;
%谱半径大于等于1就不收敛
p=max（abs（eig（B）））;
ifp>=1
sprintf（'G-Siterationmethodnotconvergent'）
%return
end
while（r>=e）&&（k~=M）
x0=x;
x=B*x0+f;
k=k+1;
r=max（abs（x-x0））;
%
%ifk<=size（Ve,1）
%Ve（k,1）=r;
%end
end
附录4SOR迭代
function[x,n]=SOR（A,b,w,e,M）
x0=zeros（size（b））;%初始解设置为与b同型的零向量
%if（w<=0||w>=2）
%error;
%return;
%end
D=diag（diag（A））;%求A的对角矩阵
L=-tril（A,-1）;%求A的下三角阵
U=-triu（A,1）;%求A的上三角阵
B=inv（D-L*w）*（（1-w）*D+w*U）;
f=w*inv（（D-L*w））*b;
n=0;%迭代次数
r=1+e;
whiler>=e
x=B*x0+f;
n=n+1;
r=max（abs（x-x0））;
%Ve（n）=r;
x0=x;
if（n>=M）
disp（'Warning:
迭代次数太多，可能不收敛！
'）;
return;
end
end
%figure
%stem（Ve（1:
100））
%title（'SOR迭代误差曲线'）;
%xlabel（'迭代次数'）;
%ylabel（'误差'）;...
附录5CG方法
function[x,n]=conjgrad（A,b,x0）
r1=b-A*x0;
p=r1;
n=0;
fori=1:
rank（A）%以下过程可参考算法流程
if（dot（p,A*p）<1.0e-50）%循环结束条件
break;
end
alpha=dot（r1,r1）/dot（p,A*p）;
x=x0+alpha*p;
r2=r1-alpha*A*p;
if（r2<1.0e-50）%循环结束条件
break;
end
belta=dot（r2,r2）/dot（r1,r1）;
p=r2+belta*p;
n=n+1;
end
附录6PCG方法
function[x,n]=preconjgrad（A,b,M,eps）
x0=zeros（size（b））;%初始解设置为与b同型的零向量
ifnargin==4
eps=1.0e-6;
end
r1=b-A*x0;
iM=inv（M）;
z1=iM*r1;
p=z1;
n=0;
tol=1;
whiletol>=eps
alpha=dot（r1,z1）/dot（p,A*p）;
x=x0+alpha*p;
r2=r1-alpha*A*p;
z2=iM*r2;
belta=dot（r2,z2）/dot（r1,z1）;
p=z2+belta*p;
n=n+1;
tol=norm（x-x0）;
x0=x;%更新迭代值
r1=r2;
z1=z2;
end
升级会员 