计算方法大作业第一次.docx
《计算方法大作业第一次.docx》由会员分享，可在线阅读，更多相关《计算方法大作业第一次.docx（15页珍藏版）》请在冰豆网上搜索。
计算方法大作业第一次
数值计算第一次大作业
实验目的以Hilbert矩阵为例，研究处理病态问题可能遇到的困难。
内容Hilbert矩阵的定义是
它是一个对称正定矩阵，而且
随着n的增加迅速增加，其逆矩阵
，这里
1）画出
之间的曲线（可以用任何的一种范数）。
你能猜出
之间有何种关系吗？
提出你的猜想并想法验证。
用行范数
forn=1:
50
fori=1:
n
forj=1:
n
A（i,j）=1/（i+j-1）;
B（i,j）=factorial（n+i-1）*factorial（n+j-1）/（（i+j-1）*（factorial（i-1）*factorial（j-1））^2*factorial（n-i）*factorial（n-j））;
end
end
result1=0;
forj=1:
n
result1=result1+A（1,j）;
end
result1=log（result1）;
result2=0;
fori=1:
n
forj=1:
n
result2=B（i,j）+result2;
end
result（i）=log（result2）;
end
m=max（result）;
x（n）=result1+m;
end
plot（[1:
50],x）
对于更大的n值，由于Hilbert逆矩阵中的元素过大，溢出，故在此取50以内的n。
图1
关系曲线图
猜想
之间存在线性关系
验证：
设
在以上程序基础上，再添加
>>;>>y=x';
>>l=1:
40;
>>k=l';
>>p=polyfit（k,y,1）%一次多项式拟合
p=
3.5446-3.0931
%P=polyfit（k,y,2）%二次多项式拟合
p=
-0.00083.5778-3.3253
%P=polyfit（k,y,3）%三次多项式拟合
0.0000-0.00333.6198-3.4777
%P=polyfit（k,y,4）%四次多项式拟合
-0.00000.0002-0.00823.6654-3.5815
%P=polyfit（k,y,5）%五次多项式拟合
p=
0.0000-0.00000.0007-0.01563.7107-3.6542
从上式可以看出，高次项系数相对于一次项和常数项系数要小很多，
所以取
2）设
是
的对角线元素开方构成的矩阵。
，不难看出
依然是对称矩阵，而且对角线元素都是1。
把
变成
的技术称为预处理。
画出
之间的曲线（可以用任何一种范数）。
你能对于预处理得出什么印象？
本小题用2范数
>>clear
>>n=500;
>>c=[];
>>fork=2:
n
H=hilb（k）;
D=diag（sqrt（diag（H）））;
D1=inv（D）;
H1=D1*H*D1;
c=[c,cond（H1）/cond（H）];
end
C=log（c）;
k=2:
n;
plot（k,C,'r-'）
图2
随n的变化曲线图
从图中给出了函数
的变化曲线。
我们观察到随着Hilbert矩阵阶数的增大，函数值在[-6,4]区间波动，主要集中在[-3,1]区间。
我们知道在
时，有
，在上图中，我们可以容易观察到，对于大部分
，函数值都是小于或者等于零的，这说明
经过预处理后的
地条件数较小，由于条件数愈大，方程组的病态愈严重，也就愈难得到方程组比较准确地解，所以预处理在一定程度上改善了原Hilbert矩阵的特性。
3）对于
，给定不同的右端项
。
分别用
以及
，求解
，比较计算结果。
取n=4，b=[1;2;3;4]
>>b=[1,2,3,4]';
>>H=hilb（4）;
>>H1=inv（H）;
>>x1=H1*b
x1=
1.0e+003*-0.06400.9000-2.52001.8200
>>D=diag（sqrt（diag（H）））;
>>D1=inv（D）;
>>H2=D1*H*D1;
>>H3=inv（H2）;
>>x2=D1*H3*D1*b
x2=
1.0e+003*-0.06400.9000-2.52001.8200
x3=H\b
x3=
1.0e+003*-0.06400.9000-2.52001.8200
同理，当n=5时，b=[1;2;3;4;5]
x1=
1.0e+004*0.0125-0.28801.4490-2.46401.3230
x2=
1.0e+004*0.0125-0.28801.4490-2.46401.3230
x3=
1.0e+004*0.0125-0.28801.4490-2.46401.3230
当n=6时，b=[1;2;3;4;5;6]
x1=
1.0e+005*-0.00220.0735-0.57121.6632-2.01600.8593
x2=
1.0e+005*-0.00220.0735-0.57121.6632-2.01600.8593
x3=
1.0e+005*-0.00220.0735-0.57121.6632-2.01600.8593
当n=7时，b=[1;2;3;4;5;6;7]
x1=
1.0e+006*0.0003-0.01610.1777-0.77281.5593-1.46360.5165
x2=
1.0e+006*0.0003-0.01610.1777-0.77281.5593-1.46360.5165
x3=
1.0e+006*0.0003-0.01610.1777-0.77281.5593-1.46360.5165
当n=8时，b=[1;2;3;4;5;6;7;8]
x1=
1.0e+007*-0.00010.0032-0.04690.2818-0.83161.2757-0.97540.2934
x2=
1.0e+007*-0.00010.0032-0.04690.2818-0.83161.2757-0.97540.2934
x3=
1.0e+007*-0.00010.0032-0.04690.2818-0.83161.2757-0.97540.2934
当n=9时，b=[1;2;3;4;5;6;7;8;9]
x1=
1.0e+007*0.0001-0.00580.1095-0.86493.4685-7.66849.4594-6.09521.5972
x2=
1.0e+007*0.0001-0.00580.1095-0.86493.4685-7.66849.4594-6.09521.5972
x3=
1.0e+007*0.0001-0.00580.1095-0.86493.4685-7.66859.4595-6.09521.5972
当n=10时，b=[1;2;3;4;5;6;7;8;9;10]
x1=
1.0e+008*-0.00000.0010-0.02330.2330-1.21073.5942-6.32246.5105-3.62280.8406
x2=
1.0e+008*-0.00000.0010-0.02330.2330-1.21073.5943-6.32276.5108-3.62290.8406
x3=
1.0e+008*-0.00000.0010-0.02330.2330-1.21083.5947-6.32336.5114-3.62320.8407
当n=11时，b=[1;2;3;4;5;6;7;8;9;10;11]
x1=
1.0e+009*0.0000-0.00020.0046-0.05640.3674-1.39913.2758-4.77254.2140-2.06290.4294
x2=
1.0e+009*0.0000-0.00020.0046-0.05670.3687-1.40383.2861-4.78674.2259-2.06850.4305
x3=
1.0e+009*0.0000-0.00020.0046-0.05670.3687-1.40373.2858-4.78624.2255-2.06820.4305
当n=12时，b=[1;2;3;4;5;6;7;8;9;10;11;12]
Warning:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=2.692153e-017.
x1=
1.0e+010*-0.00000.0000-0.00070.0110-0.08860.4225-1.26862.4585-3.06912.3821-1.04520.1980
x2=
1.0e+010*-0.00000.0000-0.00090.0133-0.10550.4975-1.47882.8408-3.51902.7127-1.18310.2229
x3=
1.0e+010*-0.00000.0000-0.00080.0123-0.09780.4630-1.38172.6636-3.30992.5586-1.11870.2113
由此可见，当n较小时，三种方法得出的结果基本相同，随着n的增大，三种方法得出的结果的偏差也越来越大。
4）取不同的
并以
的第一列为右端向量
，用高斯-塞德尔迭代法求解
，观察其收敛性。
最后你能对于有关Hilbert矩阵的计算得出哪些结论。
输入：
A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，max为最大迭代次数，w为误差精度
输出：
x为求得的方程组的解构成的列向量，n为迭代次数
取n=4，则
>>A=hilb（4）;
>>b=A（:
1）；
>>X=[0;0;0;0];
>>max=50;
>>w=10^-6;
>>gaussseidel（A,b,X,max,w）
迭代次数为
n=
2
方程组的解为
x=
1.0000
-0.0000
0.0000
-0.0000
取n=8，则
>>X=[0;0;0;0;0;0;0;0];
>>max=50;
>>w=10^-6;
>>gaussseidel（A,b,X,max,w）
迭代次数为
n=
2
方程组的解为
x=
1.0000
-0.0000
0.0000
-0.0000
0.0000
-0.0000
-0.0000
-0.0000
>>A=hilb（1000）;
>>b=A（:
1）;
>>X=zeros（1000,1）;
>>max=50;
>>w=10^-6;
>>gaussseidel（A,b,X,max,w）
迭代次数为
n=
2
迭代结果为第一行为1，其余为0的向量。
取元素全为一的向量作为起始迭代向量时，有：
>>A=hilb（4）;
>>b=A（:
1）;
>>X=ones（4,1）;
>>max=50;
>>w=10^-6;
>>gaussseidel（A,b,X,max,w）
在最大迭代次数内不收敛!
最大迭代次数后的结果为
x=
1.0206
-0.0457
-0.0884
0.1310
>>A=hilb（8）;
>>b=A（:
1）;
>>X=ones（8,1）;
>>max=50;
>>w=10^-6;
>>gaussseidel（A,b,X,max,w）
在最大迭代次数内不收敛!
最大迭代次数后的结果为
x=
0.9510
0.2777
-0.1542
-0.2055
-0.1283
-0.0147
0.1016
0.2091
>>A=hilb（200）;
>>b=A（:
1）;
>>X=ones（200,1）;
>>max=50;
>>w=10^-6;
>>gaussseidel（A,b,X,max,w）
在最大迭代次数内不收敛!
最大迭代次数后的结果为
x=
0.6054
1.3942
0.2032
-0.4150……………（由于元素太多故不一一列出）
通过推导，不难发现，起始的迭代向量为零时，对于不同的n值，均只需迭代两次就可以得出答案，且结果均是第一个元素为1，其余元素为0。
这是因为右端向量b取的是Hilbert矩阵的第一列。
如果b取其它向量，则可以知道用迭代法求解时的求解误差比较大，不收敛。
当起始的迭代向量不为零时，可以得到不同的答案，如上题中的取元素全为一的向量时，结果是发散。
由此可以知道对于不同的初始条件，迭代结果也不同。
所以用高斯-赛德尔迭代法解此方程组并不是对任意的初始向量和右端向量都收敛。
Hilbert矩阵是病态的，对于对原Hilbert矩阵进行预处理后的新Hilbert矩阵的条件数在一定范围内呈波动的变化规律。
从整体上观察，对于大多数的n，进行预处理后的Hilbert矩阵地条件数有明显的降低，这就说明预处理改善了大多数Hilbert矩阵的性质。
用迭代法求解方程时，如果系数矩阵为Hilbert矩阵，则求解结果的误差较大，根据迭代法基本定理可知用高斯-赛德尔迭代法解系数矩阵为Hilbert矩阵时是不收敛的。
当n越大时，Hn的病态越严重。
附录
function[n,x]=gaussseidel（A,b,X,nm,w）
%用高斯-赛德尔迭代法求解方程组Ax=b
%输入：
A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，max为最大迭代次数，w为误差精度
%输出：
x为求得的方程组的解构成的列向量，n为迭代次数
n=1;
m=length（A）;
I=eye（m）;%生成m*m阶的单位矩阵
D=diag（diag（A））;%令A=D-L-U,计算矩阵D
L=tril（-A）+D;%令A=D-L-U,计算矩阵L
U=triu（-A）+D;%令A=D-L-U,计算矩阵U
M=inv（D-L）*U;%计算迭代矩阵
g=inv（I-inv（D）*L）*（inv（D）*b）;%计算迭代格式中的常数项
%下面是迭代过程
whilen<=max
x=M*X+g;%用迭代格式进行迭代
ifnorm（x-X,2）disp（'迭代次数为'）;n
disp（'方程组的解为'）;x
return;
%上面：
达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=x;n=n+1;
end
%下面：
如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp（'在最大迭代次数内不收敛!
'）;
disp（'最大迭代次数后的结果为'）;x
升级会员 