数值求解Hilbert病态线性方程组.docx
《数值求解Hilbert病态线性方程组.docx》由会员分享，可在线阅读，更多相关《数值求解Hilbert病态线性方程组.docx（16页珍藏版）》请在冰豆网上搜索。
数值求解Hilbert病态线性方程组
病态线性代数方程组的求解
理论的分析表明，求解病态的线性代数方程组是困难的。
考虑方程组Hx=
b的求解，其中H为Hilbert矩阵，
，
，
1.估计Hilbert矩阵2-条件数与阶数的关系；
2.选择问题的不同维数，分别用Gauss消去法，Jacobi迭代，GS迭代和SOR迭代求解，比较结果；
3.讨论病态问题求解的算法。
解：
1、
取Hilbert矩阵阶数最高分别为n=20和n=100。
采用Hilbert矩阵的2-条件数作为实验的比较对象，画出的曲线如下图所示：
从图中可以看出，在n≤13之前，图像接近直线，在n＞13之后，图像趋于平缓，在一定的范围内上下波动。
为了比较图像的线性部分，作出一条直线与已知曲线进行比较。
比较直线的关系式为：
，结果下图所示。
从图2中可以看出，当n较小时，
之间近似满足线性关系。
当n继续增大到100时，
关系图下图所示：
从图中可以看出，图像的走势符合在n=20时的猜想，在n大于一定的值之后，图像趋于平缓，且在一定范围内震荡，同时又有一定上升趋势，但上升速度很慢。
2、
选择不同的阶数n，设方程组的精确解为xz=（1,1,…,1）T
进行计算，用四种方法解x_Guass1、x_Jacobi1、x_GS1、x_SOR1对比表如下表所示。
四种方法解x_Guass、x_Jacobi、x_GS、x_SOR对比表
n
x_Guass
x_Jacobi
x_GS
x_SOR
3
1
-8.22E+307
0.99999
1
1
-1.552e+308
1
0.99999
1
-Inf
0.99997
1
4
1
1.0172e+308
1
0.99999
1
Inf
0.99965
1.0001
1
Inf
1.0008
0.99974
1
Inf
0.99948
1.0002
5
1
-6.07E+307
0.99879
1.0001
1
-1.2685e+308
1.0223
0.99904
1
-1.6592e+308
0.90511
1.004
1
-Inf
1.1421
0.99411
1
-Inf
0.93095
1.0028
6
1
-6.64E+307
0.99949
1
1
-1.4322e+308
1.0059
0.99954
1
-Inf
0.99518
1.0011
1
-Inf
0.95149
1.0009
1
-Inf
1.1021
0.99621
1
-Inf
0.94536
1.0023
7
1
-9.20E+307
1.0007
1
1
-Inf
0.98192
1.0002
1
-Inf
1.1022
0.9979
1
-Inf
0.80824
1.0062
1
-Inf
1.062
0.99464
1
-Inf
1.1448
0.99898
1
-Inf
0.89956
1.0021
8
1
1.4378e+308
1.0019
0.99997
1
Inf
0.95875
1.001
1
Inf
1.1962
0.99408
1
Inf
0.73024
1.0114
1
Inf
0.95449
0.99464
1
Inf
1.1824
0.99806
1
Inf
1.1466
0.99668
1
Inf
0.82804
1.0042
9
1
5.09E+307
1.0031
0.99992
1
1.1719e+308
0.93946
1.0019
1
1.6249e+308
1.2489
0.99029
1
Inf
0.75446
1.0154
0.99999
Inf
0.83824
0.99549
1
Inf
1.0985
0.99876
0.99998
Inf
1.2247
0.99405
1
Inf
1.1148
0.99969
1
Inf
0.77561
1.0046
10
1
-7.82E+307
1.0036
0.99994
1
-Inf
0.9376
1.0011
1
-Inf
1.2142
0.99555
1
-Inf
0.88982
1.0045
0.99992
-Inf
0.77999
1.0011
1.0002
-Inf
0.95825
1.0003
0.99964
-Inf
1.1609
0.99748
1.0003
-Inf
1.2207
0.99776
0.99982
-Inf
1.0817
0.99927
1
-Inf
0.75019
1.003
11
1
7.97E+307
1.0029
0.99995
1
Inf
0.95857
1.0007
0.99998
Inf
1.0913
0.99809
1.0002
Inf
1.0909
0.9998
0.99894
Inf
0.80291
1.003
1.0037
Inf
0.82849
1.001
0.99187
Inf
1.0236
0.99936
1.0111
Inf
1.1886
0.99776
0.99081
Inf
1.2128
0.99799
1.0042
Inf
1.0601
0.99961
0.99916
Inf
0.73672
1.0028
12
1
8.11E+307
1.0012
0.99997
1
Inf
0.99592
1.0003
0.99993
Inf
0.91904
1.0002
1.0009
Inf
1.2974
0.99621
0.99316
Inf
0.88753
1.004
1.0298
Inf
0.74309
1.0014
0.91765
Inf
0.87451
1.001
1.1481
Inf
1.0837
0.99834
0.82734
Inf
1.2215
0.99782
1.1258
Inf
1.2148
0.99803
0.9479
Inf
1.0431
0.99993
1.0094
Inf
0.71515
1.0029
13
1
4.90E+307
0.99898
0.99999
0.99997
1.1863e+308
1.0388
0.9999
1.0012
1.6981e+308
0.7414
1.0022
0.98043
Inf
1.4694
0.99325
1.1771
Inf
1.003
1.0045
0.033234
Inf
0.70616
1.0015
4.3909
Inf
0.74649
1.0023
-6.8988
Inf
0.95394
0.99913
13.35
Inf
1.1601
0.99821
-11.81
Inf
1.2635
0.99739
9.4534
Inf
1.2194
0.99831
-2.2127
Inf
1.0198
1.0002
1.5352
Inf
0.67595
1.0032
14
1
5.94E+307
0.99677
1
0.99999
1.4525e+308
1.0799
0.99951
1.0004
Inf
0.5837
1.0041
0.99479
Inf
1.5915
0.99068
1.0373
Inf
1.1284
1.0046
0.85719
Inf
0.71157
1.0015
1.2465
Inf
0.65091
1.0035
1.176
Inf
0.82424
1
-0.85802
Inf
1.0613
0.9989
5.3287
Inf
1.2421
0.99733
-4.4541
Inf
1.302
0.99748
5.0263
Inf
1.2166
0.99841
-0.64129
Inf
0.98559
1.0006
1.2863
Inf
0.62196
1.0035
21
1
Inf
0.98949
1
0.99978
Inf
1.1624
0.99929
1.0084
Inf
0.55916
1.0017
0.86586
Inf
0.93832
1.0014
2.1174
Inf
1.5357
0.9961
-4.2693
Inf
1.4454
0.99857
15.159
Inf
1.0507
0.99966
-18.28
Inf
0.69866
1.0017
6.2524
Inf
0.52654
1.0021
13.831
Inf
0.53833
1.0021
10.008
Inf
0.68185
1.0014
-49.838
Inf
0.89247
1.0004
45.871
Inf
1.1126
0.99942
-27.1
Inf
1.2984
0.99861
54.555
Inf
1.4197
0.99808
-34.976
Inf
1.4585
0.99793
-31.788
Inf
1.4057
0.99818
24.491
Inf
1.2593
0.99885
41.386
Inf
1.0215
0.99993
-43.821
Inf
0.69776
1.0014
13.526
Inf
0.29537
1.0032
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、
关于病态问题的求解，主要的方法是对原方程作某些预处理，以降低系数矩阵的条件数。
例如选择非奇异矩阵
，使方程组
化为等价方程组
，原方程的解
.原则上应使矩阵
的条件数比
有所改善。
一般
和
可选择为三角形矩阵或对角矩阵。
理论上最好选择对角阵
和
，满足：
。
以Hilbert矩阵为例。
上图中的
（Dn为Hn的对角线元素开方构成的对角矩阵），n为希尔伯特矩阵的阶数。
我们观察到，随希尔伯特矩阵阶数的增大，函数
在[-3,1.5]区间波动，主要集中在[-1.5,0.5]区间。
我们知道当
时，有
，在上图中，我们可以很容易的观察到，对于大部分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）
升级会员 