数值分析希尔伯特病态线性方程组.docx
《数值分析希尔伯特病态线性方程组.docx》由会员分享，可在线阅读，更多相关《数值分析希尔伯特病态线性方程组.docx（13页珍藏版）》请在冰点文库上搜索。
数值分析希尔伯特病态线性方程组
病态线性代数方程组的求解
理论的分析表明，求解病态的线性代数方程组是困难的。
考虑方程组Hx=b的求解，其中H为Hilbert矩阵，H=（hij）n⨯n，hij=1，i,j=1,2,...,ni+j-1
1.估计Hilbert矩阵2-条件数与阶数的关系；
2.选择问题的不同维数，分别用Gauss消去法，Jacobi迭代，GS迭代和SOR迭代求解，比较结果；
3.讨论病态问题求解的算法。
解：
1、
取Hilbert矩阵阶数最高分别为n=20和n=100。
采用Hilbert矩阵的2-条件数作为实验的比较对象，画出的曲线如下图所示：
lg（cond（Hn））n
lg（cond（Hn））~n关系图
lg（cond（Hn））n
从图中可以看出，在n≤13之前，图像接近直线，在n＞13之后，图像趋于平缓，在一定的范围内上下波动。
为了比较图像的线性部分，作出一条直线与已知曲线进行比较。
比较直线的关系式为：
lg（cond（Hn））=1.519n-1.833，结果下图所示。
lg（cond（Hn））~n
关系图
lg（cond（Hn））n
从图2中可以看出，当n较小时，lg（cond（Hn））~n之间近似满足线性关系。
当n继续增大到100时，lg（cond（Hn））~n关系图下图所示：
lg（cond（Hn））~n关系图
lg（cond（Hn））n
从图中可以看出，图像的走势符合在n=20时的猜想，在n大于一定的值之后，图像趋于平缓，且在一定范围内震荡，同时又有一定上升趋势，但上升速度很慢。
2、
选择不同的阶数n，设方程组的精确解为xz=（1,1,…,1）T
进行计算，用四种方法解x_Guass1、x_Jacobi1、x_GS1、x_SOR1对比表如下
表所示。
Gauss消去法求解：
选择问题的阶数为3~8时，用Gauss消去法求得的解与精确解一致，当阶数为9~14时，解开始出现偏差，而且n越大，偏差越大。
用迭代法求解：
取初始向量为1.2（1,1,…,1）T.
无论n为多少阶，用Jacobi迭代方法迭代出现发散的不稳定现象，无法求解；用GS迭代方法迭代不发散，能求得解，但收敛非常缓慢，当迭代次数取得相当大（20000次）时解仍在精确解附近波动；
取w=1.5，用SOR迭代方法迭代不发散，能求得解，收敛速度较GS迭代快一些，但仍非常缓慢。
选择问题的阶数为21时：
用Gauss消去法求得的解与精确解相差很大，相差10的量级。
用迭代法求解时，取初始向量为1.2（1,1,…,1）T.
用Jacobi迭代方法迭代发散，无法求解；
用GS迭代方法迭代不发散，能求得解，但收敛非常缓慢，迭代20000次后，算得的值与精确值1相差很大。
用SOR迭代方法迭代不发散，能求得解，收敛速度还行。
从上面的结果可以看出当病态问题的阶数升高时作为直接法的Gauss消去法又能求解变成不能求解。
用Jacobi迭代方法迭代发散，无法求解；而GS和SOR迭代法在阶数升高时仍能求解，选择合适的w，可使SOR迭代法收敛速度加快。
但在阶数较低时直接法能求得精确解而迭代法却总存在一定的误差。
可见直接法与迭代法各有各的优势与不足。
3、
关于病态问题的求解，主要的方法是对原方程作某些预处理，以降低系数矩阵的条件数。
例如选择非奇异矩阵P,Q∈Rn⨯n，使方程组Ax=b化为等价方程组
（PAQ）y=Pb，原方程的解x=Qy.原则上应使矩阵PAQ的条件数比A有所改善。
一般P和Q可选择为三角形矩阵或对角矩阵。
理论上最好选择对角阵P和Q，满足：
cond（PAQ）=mincond（PAQ）。
以Hilbert矩阵为例。
lg（cond（hn）/cond（Hn））~n关系图
lg（cond（hn））n
-1-1上图中的hn=DnHnDn（Dn为Hn的对角线元素开方构成的对角矩阵），n
为希尔伯特矩阵的阶数。
我们观察到，随希尔伯特矩阵阶数的增大，函数lg[cond（hn）/cond（Hn）]在[-3,1.5]区间波动，主要集中在[-1.5,0.5]区间。
我们知道当cond（hn）≤cond（Hn）时，有lg[cond（hn）/cond（Hn）]≤0，在上图中，我们可以很容易的观察到，对于大部分n，函数值都是小于或者等于零的，这说明Hn经过预处理后的hn的条件数较小，在一定程度上改善了原希尔伯特矩阵的性质。
由于希尔伯特矩阵病态的性质，对于对原希尔伯特矩阵进行预处理后的新希尔伯特矩阵的条件数在一定的区间呈波动的变化规律。
从整体上观察，对于大多数的n，进行预处理后的希尔伯特矩阵的条件数有明显的降低，这就说明预处理改善了大多数阶数的希尔伯特矩阵的性质。
第1小题Matlab编程如下：
t1=20
i1=1:
t1
forn=1:
t1
H=hilb（n）
K1（n）=cond（H）
K（n）=log10（K1（n））
end
l1=1.519*i1-1.833
plot（i1,K（i1）,i1,l1）
xlabel（'n','fontsize',20）
ylabel（'lg（cond（Hn））','fontsize',20）
title（'lg（cond（Hn））~n关系图','fontsize',24）set（gca,'fontsize',16）
figure
t2=100
i2=1:
t2
forn=1:
t2
H=hilb（n）
K1（n）=cond（H）
K（n）=log10（K1（n））
end
plot（i2,K（i2））
xlabel（'n','fontsize',20）
ylabel（'lg（cond（Hn））','fontsize',20）
title（'lg（cond（Hn））~n关系图','fontsize',24）set（gca,'fontsize',16）
第2小题Matlab编程如下：
clear
form=3:
25
n=m%Hilbert矩阵的阶数
n1=n+1
H=hilb（n）
K=cond（H）
u=1.2
x0=u*ones（n,1）
xz=ones（n,1）
b=H*xz
a=H
D=diag（diag（a））
U=-triu（a,1）
L=-tril（a,-1）
w=1.5
d=1.0e-6
max=1000
c=[a,b]
%Gauss消去法
fork=1:
n
bmax=0;%选主元，并存放在变量bmax中
fori=k:
n
ifbmax-abs（c（i,k））<0
bmax=abs（c（i,k））;
l=i;%记下绝对值最大者所在位置
end
end
ifbmax<1.0e-20
'stop';pause;%如果主元绝对值小于1.0e-20，则A奇异,计算终止end
%进行行交换
ifl~=k%如果l=k,则不需要交换
forj=k:
n1
t=c（l,j）;c（l,j）=c（k,j）;c（k,j）=t;
end
end
%进行消元
t=1/c（k,k）;
forj=k+1:
n1
c（k,j）=c（k,j）*t;c（k+1:
n,j）=c（k+1:
n,j）-c（k+1:
n,k）*c（k,j）;end
end
%回代法求解上三角形方程组，解存放在矩阵C的第n+1列中fork=2:
n
i=n1-k;
forj=i+1:
n
c（i,n1）=c（i,n1）-c（i,j）*c（j,n1）;
end
end
xGuass=c（1:
n,n1）
%Jacobi迭代法
BJ=D\（L+U）
fJ=D\b
x2=x0
xJacobi=BJ*x2+fJ
kJacobi=1
while（norm（xJacobi-x2）>=d）&（kJacobix2=xJacobi
xJacobi=BJ*x2+fJ
kJacobi=kJacobi+1
end
%GS迭代法
BG=（D-L）\U
fG=（D-L）\b
x3=x0
xGS=BJ*x3+fG
kGS=1
while（norm（xGS-x3）>=d）&（kGSx3=xGS
xGS=BG*x3+fG
kGS=kGS+1
end
%SOR迭代法
lw=（D-w*L）\（（1-w）*D+w*U）
fS=（D-w*L）\b*w
x4=x0
xSOR=lw*x4+fS
kSOR=1
while（norm（xSOR-x4）>=d）&（kSORx4=xSOR
xSOR=lw*x4+fS
kSOR=kSOR+1
end
forp=1:
m
xGuass1（p,m）=xGuass（p,1）
xJacobi1（p,m）=xJacobi（p,1）
xGS1（p,m）=xGS（p,1）
xSOR1（p,m）=xSOR（p,1）
x（p,4*m-3）=xGuass1（p,m）
x（p,4*m-2）=xJacobi1（p,m）
x（p,4*m-1）=xGS1（p,m）
x（p,4*m）=xSOR1（p,m）
end
end
第3小题Matlab编程如下：
clear
k=500
i=1:
k
forn=1:
k
H=hilb（n）
K（n）=cond（H）
D1=diag（sqrt（diag（H）））
D=inv（D1）
Hy=D*H*D
Ky1（n）=cond（Hy）
Ky（n）=log10（Ky1（n）/K（n））
end
plot（i,Ky）
xlabel（'n','fontsize',20）
ylabel（'lg（cond（hn））','fontsize',20）
title（'lg（cond（hn）/cond（Hn））~n关系图','fontsize',24）set（gca,'fontsize',16）
升级会员 