数值求解Hilbert病态线性方程组Word文件下载.docx
《数值求解Hilbert病态线性方程组Word文件下载.docx》由会员分享，可在线阅读，更多相关《数值求解Hilbert病态线性方程组Word文件下载.docx（16页珍藏版）》请在冰豆网上搜索。
1.0223
0.99904
-1.6592e+308
0.90511
1.004
1.1421
0.99411
0.93095
1.0028
6
-6.64E+307
0.99949
-1.4322e+308
1.0059
0.99954
0.99518
1.0011
0.95149
1.0009
1.1021
0.99621
0.94536
1.0023
7
-9.20E+307
1.0007
0.98192
1.1022
0.9979
0.80824
1.0062
1.062
0.99464
1.1448
0.99898
0.89956
1.0021
8
1.4378e+308
1.0019
0.95875
1.001
1.1962
0.99408
0.73024
1.0114
0.95449
1.1824
0.99806
1.1466
0.99668
0.82804
1.0042
9
5.09E+307
1.0031
0.99992
1.1719e+308
0.93946
1.6249e+308
1.2489
0.99029
0.75446
1.0154
0.83824
0.99549
1.0985
0.99876
0.99998
1.2247
0.99405
1.1148
0.99969
0.77561
1.0046
10
-7.82E+307
1.0036
0.99994
0.9376
1.2142
0.99555
0.88982
1.0045
0.77999
0.95825
1.0003
0.99964
1.1609
0.99748
1.2207
0.99776
0.99982
1.0817
0.99927
0.75019
1.003
11
7.97E+307
1.0029
0.99995
0.95857
1.0913
0.99809
1.0909
0.9998
0.99894
0.80291
1.0037
0.82849
0.99187
1.0236
0.99936
1.0111
1.1886
0.99081
1.2128
0.99799
1.0601
0.99961
0.99916
0.73672
12
8.11E+307
1.0012
0.99592
0.99993
0.91904
1.2974
0.99316
0.88753
1.0298
0.74309
1.0014
0.91765
0.87451
1.1481
1.0837
0.99834
0.82734
1.2215
0.99782
1.1258
1.2148
0.99803
0.9479
1.0431
1.0094
0.71515
13
4.90E+307
1.1863e+308
1.0388
0.9999
1.6981e+308
0.7414
1.0022
0.98043
1.4694
0.99325
1.1771
0.033234
0.70616
1.0015
4.3909
0.74649
-6.8988
0.95394
0.99913
13.35
1.1601
0.99821
-11.81
1.2635
0.99739
9.4534
1.2194
0.99831
-2.2127
1.0198
1.5352
0.67595
1.0032
14
5.94E+307
0.99677
1.4525e+308
1.0799
0.99951
1.0004
0.5837
1.0041
0.99479
1.5915
0.99068
1.0373
1.1284
0.85719
0.71157
1.2465
0.65091
1.0035
1.176
0.82424
-0.85802
1.0613
0.9989
5.3287
1.2421
0.99733
-4.4541
1.302
5.0263
1.2166
0.99841
-0.64129
0.98559
1.0006
1.2863
0.62196
21
0.98949
0.99978
1.1624
0.99929
1.0084
0.55916
1.0017
0.86586
0.93832
2.1174
1.5357
0.9961
-4.2693
1.4454
0.99857
15.159
1.0507
0.99966
-18.28
0.69866
6.2524
0.52654
13.831
0.53833
10.008
0.68185
-49.838
0.89247
45.871
1.1126
0.99942
-27.1
1.2984
0.99861
54.555
1.4197
0.99808
-34.976
1.4585
0.99793
-31.788
1.4057
0.99818
24.491
1.2593
0.99885
41.386
1.0215
-43.821
0.69776
13.526
0.29537
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消去法又能求解变成不能求解。
而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:
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））'
title（'
lg（cond（Hn））~n关系图'
24）
set（gca,'
16）
figure
t2=100
i2=1:
t2
plot（i2,K（i2））
第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:
bmax=0;
%选主元，并存放在变量bmax中
fori=k:
ifbmax-abs（c（i,k））<
bmax=abs（c（i,k））;
l=i;
%记下绝对值最大者所在位置
end
ifbmax<
1.0e-20
'
stop'
;
pause;
%如果主元绝对值小于1.0e-20，则A奇异,计算终止
%进行行交换
ifl~=k%如果l=k,则不需要交换
forj=k:
n1
t=c（l,j）;
c（l,j）=c（k,j）;
c（k,j）=t;
%进行消元
t=1/c（k,k）;
forj=k+1:
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）;
%回代法求解上三角形方程组，解存放在矩阵C的第n+1列中
fork=2:
i=n1-k;
forj=i+1:
c（i,n1）=c（i,n1）-c（i,j）*c（j,n1）;
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）&
（kJacobi<
max）
x2=xJacobi
xJacobi=BJ*x2+fJ
kJacobi=kJacobi+1
%GS迭代法
BG=（D-L）\U
fG=（D-L）\b
x3=x0
xGS=BJ*x3+fG
kGS=1
while（norm（xGS-x3）>
（kGS<
x3=xGS
xGS=BG*x3+fG
kGS=kGS+1
%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）>
（kSOR<
x4=xSOR
xSOR=lw*x4+fS
kSOR=kSOR+1
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）
第3小题Matlab编程如下：
k=500
i=1:
k
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）
lg（cond（hn））'
lg（cond（hn）/cond（Hn））~n关系图'
升级会员 