计算方法大作业第一次.docx

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)*factori

2、al(n-j);endendresult1=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

3、=polyfit(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

5、,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*bx1 = 1.0e+003 * -0.0640 0.9000 -2.5200 1.8200 D=diag(sqrt(diag(H); D1=inv

6、(D); H2=D1*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 -

7、2.4640 1.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.5165

8、x2 = 1.0e+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

9、 0.2934x3 = 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.005

10、8 0.1095 -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

11、 0.2330 -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

12、 -0.0002 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.

13、3821 -1.0452 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; gaus

15、sseidel(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

16、,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

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

18、收敛。Hilbert 矩阵是病态的，对于对原Hilbert 矩阵进行预处理后的新Hilbert 矩阵的条件数在一定范围内呈波动的变化规律。从整体上观察，对于大多数的n，进行预处理后的Hilbert 矩阵地条件数有明显的降低，这就说明预处理改善了大多数Hilbert 矩阵的性质。用迭代法求解方程时，如果系数矩阵为Hilbert 矩阵，则求解结果的误差较大，根据迭代法基本定理可知用高斯-赛德尔迭代法解系数矩阵为Hilbert 矩阵时是不收敛的。当n越大时，Hn的病态越严重。附录function n,x=gaussseidel(A,b,X,nm,w)%用高斯-赛德尔迭代法求解方程组Ax=b%输入：A

19、为方程组的系数矩阵，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=D-L-U,计算矩阵UM=inv(D-L)*U; %计算迭代矩阵g=inv(I-inv(D)*L)*(inv(D)*b); %计算迭代格式中的常数项%下面是迭代过程while n=max x=M*X+g; %用迭代格式进行迭代 if norm(x-X,2)w disp(迭代次数为);n disp(方程组的解为);x return; %上面：达到精度要求就结束程序，输出迭代次数和方程组的解 end X=x;n=n+1;end%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）disp(在最大迭代次数内不收敛!);disp(最大迭代次数后的结果为);x