计算方法大作业第一次

1、数值计算第一次大作业实验目的 以Hilbert矩阵为例，研究处理病态问题可能遇到的困难。内容 Hilbert矩阵的定义是它是一个对称正定矩阵，而且随着n的增加迅速增加，其逆矩阵，这里1) 画出之间的曲线（可以用任何的一种范数）。你能猜出之间有何种关系吗？提出你的猜想并想法验证。用行范数for n=1:50for i=1:nfor j=1:nA(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);e

2、ndendresult1=0;for j=1:nresult1=result1+A(1,j);endresult1=log(result1);result2=0;for i=1:nfor j=1:nresult2=B(i,j)+result2;endresult(i)=log(result2);endm=max(result);x(n)=result1+m;endplot(1:50,x)对于更大的n值，由于Hilbert逆矩阵中的元素过大，溢出，故在此取50以内的n。图1 关系曲线图猜想之间存在线性关系验证：设在以上程序基础上，再添加; y=x; l=1:40; k=l; p=polyfit(

3、k,y,1) %一次多项式拟合p =3.5446 -3.0931% P=polyfit(k,y,2) %二次多项式拟合p = -0.0008 3.5778 -3.3253% P=polyfit(k,y,3) %三次多项式拟合 0.0000 -0.0033 3.6198 -3.4777% P=polyfit(k,y,4) %四次多项式拟合-0.0000 0.0002 -0.0082 3.6654 -3.5815% P=polyfit(k,y,5) %五次多项式拟合p = 0.0000 -0.0000 0.0007 -0.0156 3.7107 -3.6542从上式可以看出，高次项系数相对于一次项

4、和常数项系数要小很多，所以取2）设是的对角线元素开方构成的矩阵。，不难看出依然是对称矩阵，而且对角线元素都是1。把变成的技术称为预处理。画出之间的曲线（可以用任何一种范数）。你能对于预处理得出什么印象？本小题用2范数 clear n=500; c=; for k=2:nH=hilb(k);D=diag(sqrt(diag(H);D1=inv(D);H1=D1*H*D1;c=c,cond(H1)/cond(H);endC=log(c);k=2:n;plot(k,C,r-)图 2 随n的变化曲线图从图中给出了函数的变化曲线。我们观察到随着Hilbert矩阵阶数的增大，函数值在-6,4区间波动，主要

5、集中在-3,1区间。我们知道在时，有，在上图中，我们可以容易观察到，对于大部分，函数值都是小于或者等于零的，这说明经过预处理后的地条件数较小，由于条件数愈大，方程组的病态愈严重，也就愈难得到方程组比较准确地解，所以预处理在一定程度上改善了原Hilbert矩阵的特性。3）对于，给定不同的右端项。分别用以及，求解，比较计算结果。取n=4，b=1;2;3;4 b=1,2,3,4; H=hilb(4); H1=inv(H); x1=H1*bx1 = 1.0e+003 * -0.0640 0.9000 -2.5200 1.8200 D=diag(sqrt(diag(H); D1=inv(D); H2=D

6、1*H*D1; H3=inv(H2); x2=D1*H3*D1*bx2 = 1.0e+003 * -0.0640 0.9000 -2.5200 1.8200x3=Hbx3 = 1.0e+003 * -0.0640 0.9000 -2.5200 1.8200同理，当n=5时，b=1;2;3;4;5x1 = 1.0e+004 * 0.0125 -0.2880 1.4490 -2.4640 1.3230x2 = 1.0e+004 * 0.0125 -0.2880 1.4490 -2.4640 1.3230x3 = 1.0e+004 * 0.0125 -0.2880 1.4490 -2.4640 1.

7、3230当n=6时，b=1;2;3;4;5;6x1 = 1.0e+005 * -0.0022 0.0735 -0.5712 1.6632 -2.0160 0.8593x2 = 1.0e+005 * -0.0022 0.0735 -0.5712 1.6632 -2.0160 0.8593x3 = 1.0e+005 * -0.0022 0.0735 -0.5712 1.6632 -2.0160 0.8593当n=7时，b=1;2;3;4;5;6;7x1 = 1.0e+006 * 0.0003 -0.0161 0.1777 -0.7728 1.5593 -1.4636 0.5165x2 = 1.0e

8、+006 * 0.0003 -0.0161 0.1777 -0.7728 1.5593 -1.4636 0.5165x3 = 1.0e+006 * 0.0003 -0.0161 0.1777 -0.7728 1.5593 -1.4636 0.5165当n=8时，b=1;2;3;4;5;6;7;8x1 = 1.0e+007 * -0.0001 0.0032 -0.0469 0.2818 -0.8316 1.2757 -0.9754 0.2934x2 = 1.0e+007 * -0.0001 0.0032 -0.0469 0.2818 -0.8316 1.2757 -0.9754 0.2934x3

9、 = 1.0e+007 * -0.0001 0.0032 -0.0469 0.2818 -0.8316 1.2757 -0.9754 0.2934当n=9时，b=1;2;3;4;5;6;7;8;9x1 = 1.0e+007 * 0.0001 -0.0058 0.1095 -0.8649 3.4685 -7.6684 9.4594 -6.0952 1.5972x2 = 1.0e+007 * 0.0001 -0.0058 0.1095 -0.8649 3.4685 -7.6684 9.4594 -6.0952 1.5972x3 = 1.0e+007 * 0.0001 -0.0058 0.1095

10、-0.8649 3.4685 -7.6685 9.4595 -6.0952 1.5972当n=10时，b=1;2;3;4;5;6;7;8;9;10x1 = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -1.2107 3.5942 -6.3224 6.5105 -3.6228 0.8406x2 = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -1.2107 3.5943 -6.3227 6.5108 -3.6229 0.8406x3 = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -

11、1.2108 3.5947 -6.3233 6.5114 -3.6232 0.8407当n=11时，b=1;2;3;4;5;6;7;8;9;10;11x1 = 1.0e+009 * 0.0000 -0.0002 0.0046 -0.0564 0.3674 -1.3991 3.2758 -4.7725 4.2140 -2.0629 0.4294x2 = 1.0e+009 * 0.0000 -0.0002 0.0046 -0.0567 0.3687 -1.4038 3.2861 -4.7867 4.2259 -2.0685 0.4305x3 = 1.0e+009 * 0.0000 -0.0002

12、0.0046 -0.0567 0.3687 -1.4037 3.2858 -4.7862 4.2255 -2.0682 0.4305当n=12时，b=1;2;3;4;5;6;7;8;9;10;11;12 Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.692153e-017.x1 = 1.0e+010 * -0.0000 0.0000 -0.0007 0.0110 -0.0886 0.4225 -1.2686 2.4585 -3.0691 2.3821 -1.0

13、452 0.1980x2 = 1.0e+010 * -0.0000 0.0000 -0.0009 0.0133 -0.1055 0.4975 -1.4788 2.8408 -3.5190 2.7127 -1.1831 0.2229x3 = 1.0e+010 * -0.0000 0.0000 -0.0008 0.0123 -0.0978 0.4630 -1.3817 2.6636 -3.3099 2.5586 -1.1187 0.2113由此可见，当n较小时，三种方法得出的结果基本相同，随着n的增大，三种方法得出的结果的偏差也越来越大。4）取不同的并以的第一列为右端向量，用高斯-塞德尔迭代法求解

14、，观察其收敛性。最后你能对于有关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

15、,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)在最

16、大迭代次数内不收敛!最大迭代次数后的结果为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)在最大迭代次

17、数内不收敛!最大迭代次数后的结果为x = 0.6054 1.3942 0.2032 -0.4150（由于元素太多故不一一列出）通过推导，不难发现，起始的迭代向量为零时，对于不同的n值，均只需迭代两次就可以得出答案，且结果均是第一个元素为1，其余元素为0。这是因为右端向量b取的是Hilbert 矩阵的第一列。如果b取其它向量，则可以知道用迭代法求解时的求解误差比较大，不收敛。当起始的迭代向量不为零时，可以得到不同的答案，如上题中的取元素全为一的向量时，结果是发散。由此可以知道对于不同的初始条件，迭代结果也不同。所以用高斯-赛德尔迭代法解此方程组并不是对任意的初始向量和右端向量都收敛。Hilber

18、t 矩阵是病态的，对于对原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,计算矩阵DL=tril(-A)+D; %令A=D-L-U,计算矩阵LU=triu(-A)+D; %令A=