试验报告.doc

计算方法实 验 报 告 册姓 名： 学 号： 班 级： 实验名称： 实验所属课程： 实验室(中心)： 指 导 教 师 ： 实验完成时间： 目 录实验一 数值稳定性实验.实验二 非线性方程求根的数值解法.实验三 线性方程组的数值解法（一）.实验四 线性方程组的数值解法（二）.实验五 分段插值和拉格朗日插值实验.实验六 复合辛普森求积分实验实验名称： 实验一 数值稳定性实验实验内容：设 （1）从尽可能精确的近似值出发，利用递推式计算的近似值。（2）从较粗糙的估计值出发，利用递推式计算的近似值。（3）分析所得结果的可靠性以及出现这种现象的原因。实验过程：（1）先求的：程序代码：#includemain()double I0=0.182322,In=0,It,n=1.0;It=I0;printf(nttItn);for(n=1;n=20;)In=-5*It+1.0/n;It=In;printf(%ft%.15fn,n,It);n=n+1.0;（2）以（1）中用算的的近似最为较粗糙的估计值。程序代码：#includemain()double I20=42267422.243506,In=0,It,n=1.0;It=I20;printf(nttItn);for(n=20;n0;)In=-0.2*It+1.0/(5*n);It=In;printf(%ft%.15fn,n-1,It);n=n-1.0;实验结果分析：第一小题程序运行结果如下：但是当初始的时运行结果为 显然，此时由于初始值的微小差异使得在经过20次的计算后结果已有了很大的差异。对同一数学问题，为了求得其解，往往可以设计出多种不同的算法。而不同的算法，在执行过程中引入的舍入误差以及舍入误差的累积情况，也往往是不同的。假定初始值的误差为E，中间不断产生新误差，考察由E引起的误差累积是否增长。在第二小题中，如果用作为近似值计算，计算得到的结果为0.182322，说明算法的稳定性不是问题，之所以出现上述问题，个人认为与小数点后保留的位数和初始值的小数点位数有关。实验名称： 实验二 非线性方程求根的数值解法实验内容： 求方程在区间上的近似值，要求误差不超过。（1）利用二分法；（2）取初值，利用简单迭代法；（3）取初值，利用牛顿迭代法。实验过程： （1）二分法程序代码：#includestdio.h#includemath.hfloat function(float x)float y;y=pow(x,3)-x-1;return(y);main()float a,b,E1=0.000001,E2=0.000001,x,y,y1,y2;int k;printf(请输入a和b.n);printf(a=),scanf(%f,&a);printf(b=),scanf(%f,&b);y1=function(a);y2=function(b);printf(kttxttyn);if(y1*y20)printf(注意：输入a和b的值无法计算出结果。n);elsek=1;dox=(a+b)/2,y=function(x);if(y-E1&yE1) printf(%ft%ft%ftn,k,x,y);break;elseif(y1*yE2|b-a=E2)k=k+1;while(b-aE2|b-a=E2);（2）简单迭代法；迭代函数为x=x+113程序代码：#includestdio.h#includemath.hmain()float x0,x1,E=0.000001;int k=0,N=17;printf(请输入x0.n);printf(x0=),scanf(%f,&x0);printf(kttx1n);dok=k+1;x1=pow(x0+1,1.0/3);if(x1-x0-E&x1-x0E)printf(%dtt%fn,k,x1);break;elsex0=x1; if(k=N)printf(求算失败！n);while(k-E&x1-x0E)printf(%dtt%fn,k,x1);break;elsex0=x1;if(k=N)printf(求算失败！n);break;while(kN);实验结果分析：程序一运行结果：程序二运行结果如下：程序三运行结果如下：对三个程序的运行结果进行比较，发现牛顿迭代法有明显的优越性，在得到相同结果的前提下，用牛顿迭代法编的程序循环部分运行的次数最少。当运行一个大型的程序时，程序运行所花的时间是必须考虑的。时间太长，严重影响程序的整体性能，所以这时候要寻找一种尽量减少时间的算法。实验名称： 实验三 线性方程组的数值解法（一）实验内容：求解线性方程组 并比较计算结果精度（方程组精确解为）（1） 顺序消去法；（2） 列主元消去法。实验过程：(1) 顺序消元法解该方程组的代码如下：#includestdio.hmain()float A45=1.1348,3.8326,1.1651,3.4017,9.5342,0.5301,1.7875,2.5330,1.5435,6.3941,3.4129,4.9317,8.7643,1.3142,18.4231,1.2371,4.9998,10.6721,0.0147,16.9237;float T45;float X5;int i,j,k;printf(按顺序消元法，将方程组转化为上三角矩阵形式如下：n);for(j=0;j3;j+)for(i=1;i4;i+)Tij=Ai+jj/Ajj;for(k=j;k5;k+)Ai+jk=Ai+jk-Ajk*Tij;for(i=0;i4;i+)for(j=0;j5;j+)printf(%.4ft,Aij);printf(n);printf(方程组的解如下：n);X4=A34/A33;X3=(A24-A23*X4)/A22;X2=(A14-A13*X4-A12*X3)/A11;X1=(A04-A03*X4-A02*X3-A01*X2)/A00;printf(x1=%.4ftnx2=%.4ftnx3=%.4ftnx4=%.4ft,X1,X2,X3,X4);（2）列主元消去法解该方程的程序代码如下：#includestdio.hmain()float A45=1.1348,3.8326,1.1651,3.4017,9.5342,0.5301,1.7875,2.5330,1.5435,6.3941,3.4129,4.9317,8.7643,1.3142,18.4231,1.2371,4.9998,10.6721,0.0147,16.9237;float T45;float X5,MAX5,ZD5;int i,j,k,l,m;printf(按列组消元法，将方程组转化为上三角矩阵形式如下：n);for(j=0;j3;j+)MAXj=Ajj;for(l=1;lMAXj)MAXj=Al+jj;for(l=0;l4-j;l+)if(MAXj=Al+jj)for(m=j;m5;m+)ZDm=Ajm;Ajm=Al+jm;Al+jm=ZDm;for(i=1;i4;i+)Tij=Ai+jj/Ajj;for(k=j;k5;k+)Ai+jk=Ai+jk-Ajk*Tij;for(i=0;i4;i+)for(j=0;j5;j+)printf(%.4ft,Aij);printf(n);printf(方程组的解如下：n);X4=A34/A33;X3=(A24-A23*X4)/A22;X2=(A14-A13*X4-A12*X3)/A11;X1=(A04-A03*X4-A02*X3-A01*X2)/A00;printf(x1=%.4ftnx2=%.4ftnx3=%.4ftnx4=%.4ft,X1,X2,X3,X4);实验结果分析：(1)顺序消元法程序运行结果如下所示：(2)列组消元法程序运行结果如下所示： 在此题中，用顺序消元法解方程的时候，当把矩阵化成上三角矩阵时，第三个方程的X3的系数与其他系数相比，显得格外的大。第二个方程的X2的系数与其他系数相比，又显得格外的小。而当用列组消元法，这种情况得到了很大程度的解决，对结果的准确性有了一定的提高。 在程序运行的时候弹出下图所示的对话框：但是程序能正常运行处结果。 在改实验的程序代码中，大量的运用了FOR语句，很多运算都是通过循环来解决的。使得程序的代码行数得到了大大的降低。实验名称： 实验四 线性方程组的数值解法（二）实验内容：利用杜利特尔分解方法求解线性方程组 实验过程：程序代码如下：#includestdio.hmain()float A66=1,2,4,7,11,16,2,3,5,8,12,17,4,5,6,9,13,18,7,8,9,10,14,19,11,12,13,14,15,20,16,17,18,19,20,21,L66=1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,U66=1,2,4,7,11,16,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;float B6=10,11,12,13,14,15;float Y6,X6;int i,j;for(i=0;i6;i+) Li0=Ai0/L00;U11=(A11-L10*U01)/L11;for(i=2;i6;i+) Li1=(Ai1-Li0*U01)/U11;U1i=(A1i-L10*U0i)/L11;U2i=(A2i-L20*U0i-L21*U1i)/L22;for(i=3;i6;i+) Li2=(Ai2-Li0*U02-Li1*U12)/U22;U3i=(A3i-L30*U0i-L31*U1i-L32*U2i)/L33;for(i=4;i6;i+) Li3=(Ai3-Li0*U03-Li1*U13-Li2*U23)/U33;U4i=(A4i-L40*U0i-L41*U1i-L42*U2i-L43*U3i)/L44;L54=(A54-L50*U04-L51*U14-L52*U24-L53*U34)/U44;U55=(A55-L50*U05-L51*U15-L52*U25-L53*U35-L54*U45)/L55;printf(矩阵A：n);for(i=0;i6;i+)for(j=0;j6;j+)printf(%.0ft,Aij);printf(n);printf(矩阵L：n);for(i=0;i6;i+)for(j=0;j6;j+)printf(%.0ft,Lij);printf(n);printf(矩阵U：n);for(i=0;i6;i+)for(j=0;j6;j+)printf(%.0ft,Uij);printf(n);Y0=B0;Y1=B1-L10*Y0;Y2=B2-L20*Y0-L21*Y1;Y3=B3-L30*Y0-L31*Y1-L32*Y2;Y4=B4-L40*Y0-L41*Y1-L42*Y2-L43*Y3;Y5=B5-L50*Y0-L51*Y1-L52*Y2-L53*Y3-L54*Y4;printf(Y的值：n);printf(Y1=%.4fnY2=%.4fnY3=%.4fnY4=%.4fnY5=%.4fnY6=%.4fn,Y0,Y1,Y2,Y3,Y4,Y5);X5=Y5/U55;X4=(Y4-U45*X5)/U44;X3=(Y3-U35*X5-U34*X4)/U33;X2=(Y2-U25*X5-U24*X4-U23*X3)/U22;X1=(Y1-U15*X5-U14*X4-U13*X3-U12*X2)/U11;X0=(Y0-U05*X5-U04*X4-U03*X3-U02*X2-U01*X1)/U00;printf(方程的解X的值：n);printf(X1=%.4fnX2=%.4fnX3=%.4fnX4=%.4fnX5=%.4fnX6=%.4fn,X0,X1,X2,X3,X4,X5);实验结果分析：程序运行结果如下：实验名称： 实验五 分段插值和拉格朗日插值实验实验内容：从函数表x0.00.10.1950.30.4010.5f(x)0.398940.396950.391420.381380.368120.35206出发，用下列方法计算f(0.15), f(0.31), f(0.47)的近似值；（1） 分段线性插值（2） 分段二次插值（3） 全区间上拉格朗日插值实验过程：（1）分段线性插值的代码如下：#includestdio.hmain()float X6=0.0,0.1,0.195,0.3,0.401,0.5,Y6=0.39894,0.39695,0.39142,0.38138,0.36812,0.35206;int i,n,t,YorN;float x,y;doprintf(请输入x的值：nx=);scanf(%f,&x);for(i=0;i6;i+)if(xXi)t=i;break; y=Yi-1*(x-Xi)/(Xi-1-Xi)+Yi*(x-Xi-1)/(Xi-Xi-1);printf(近似值Y为:%fn,y);printf(是否继续进行插值计算？n请输入1 or 0(1代表继续，0代表不继续)。n);scanf(%d,&YorN);while(YorN=1);（2）分段二次插值的代码如下：#includestdio.hmain()float X6=0.0,0.1,0.195,0.3,0.401,0.5,Y6=0.39894,0.39695,0.39142,0.38138,0.36812,0.35206;int i,n,t,YorN;float x,y;float XCHZ5;doprintf(请输入x的值：nx=);scanf(%f,&x);for(i=0;i6;i+)if(xXi)t=i;break; y=Yi-2*(x-Xi-1)*(x-Xi)/(Xi-2-Xi-1)*(Xi-2-Xi)+Yi-1*(x-Xi-2)*(x-Xi)/(Xi-1-Xi-2)*(Xi-1-Xi)+Yi*(x-Xi-2)*(x-Xi-1)/(Xi-Xi-2)*(Xi-Xi-1);printf(近似值Y为:%fn,y);printf(是否继续进行插值计算？n请输入1 or 0(1代表继续，0代表不继续)。n);scanf(%d,&YorN);while(YorN=1);（3）全区间上拉格朗日插值程序代码如下：#includestdio.hmain()float X6=0.0,0.1,0.195,0.3,0.401,0.5,Y6=0.39894,0.39695,0.39142,0.38138,0.36812,0.35206;int i,n,t,YorN;float x,y;float XCHZ5,F0,F15,F24,F33,F42,F51;F0=Y0;for(i=0;i5;i+)F1i=(Yi+1-Yi)/(Xi+1-Xi);for(i=0;i4;i+)F2i=(F1i+1-F1i)/(Xi+2-Xi);for(i=0;i3;i+)F3i=(F2i+1-F2i)/(Xi+3-Xi);for(i=0;i2;i+)F4i=(F3i+1-F3i)/(Xi+4-Xi);for(i=0;i1;i+)F5i=(F4i+1-F4i)/(Xi+5-Xi);doprintf(请输入x的值：nx=);scanf(%f,&x);y=F0+F10*(x-X0)+F20*(x-X0)*(x-X1)+F30*(x-X0)*(x-X1)*(x-X2)+F40*(x-X0)*(x-X1)*(x-X2)*(x-X3)+F50*(x-X0)*(x-X1)*(x-X2)*(x-X3)*(x-X4);printf(近似值Y为:%fn,y