数值计算方法上机实验

PAGE1PAGE16重庆交通大学学生实验报告实验课程名称数值计算方法开课实验室理学院实验室学院土木建筑学院专业班级09级测绘工程专业2班学生姓名吴必坤学号09210205开课时间2010至2011学年第二学期评分细则评分考勤（10分）排版格式（10分）实验程序代码的正确性、完整性（40分）实验结果及分析（30分）实验创新（10分）总成绩教师签名张文忠目录实验一误差——定积分 2实验二非线性方程求根 4（一）牛顿迭代法 4（二）牛顿下山法 5实验三解线性代数方程组——列主元消元法 7实验四解线性代数方程组——追赶法 9实验五插值与拟合——三次样条差值 11实验一误差——定积分实验内容：用两种方法求解以下积分：,初值2.In=1-In-1*n2实验过程：1利用C语言编程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){doublex1;inti,n;x1=100;for(i=50;i>=0;i--){x1=(1-x1)/i;printf("i=%d,xi=%lf\n",i,x1);}}1.2.命令窗口中程序结果如下：2.1利用C语言编程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){doublex1=1-exp(-1);inti,n,x0;scanf("%d",&n);x0=1-exp(-1);printf("i=0,x0=0.6320205671\n");for(i=1;i<=n;i++){x1=1-i*x1;printf("i=%d,xi=%12.9lf\n",i,x1);}}1.2.命令窗口中程序结果如下：实验结果分析：从程序运行的结果可以看出，当输入n=14时，即运算到n=14时，结果与书上的结果并不一致，然而从n=1到n=10的运算结果是完全一致的，从n=11以后的运算结果可能是由于不同时期的计算机的处理数据的能力的不同或精度不同导致的。实验二非线性方程求根实验内容：用牛顿迭代法和牛顿下山法求解方程f（x）=，取近似值。1．牛顿迭代法：，直到。2．牛顿下山法：，（λ=1至合适），且直到|实验过程：1．用牛顿迭代法求解1.1利用C语言编程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){doublex=-0.99,y;intn=0;printf("%d,%lf\n",n,x);do{y=x;x=x-(x*x*x/3-x)/(x*x-1);n++;printf("n=%d,x=%lf\n",n,x);}while(fabs(y-x)>1e-5);}1.2.命令窗口结果截屏如下：2．用牛顿下山法求解2.1利用C语言编程求解，程序如下：#include<stdio.h>#include<math.h>doublef(doublex){ doublet;t=x*x*x/3-x;returnfabs(t);}voidmain(){doublex=-0.99,y;intk=0;printf("%d,%lf\n",k,x);do {y=x;intn=0;do{ x=y;x=x-pow(0.5,n)*(x*x*x/3-x)/(x*x-1); n++; }while(f(x)>=f(y));k++;printf("k=%d,n=%d,x=%lf\n",k,n,x);}while(fabs(x-y)>1e-5);}2.2.命令窗口结果截屏如下：实验结果分析：从运算截图中我们可以知道，一方面运用牛顿迭代法和牛顿下山迭代法的程序运算的结果是一致的，若出现小小的差异的话，可能是由于设置的变量类型不同导致，另一方面运用牛顿下山迭代法比运用一般牛顿迭代法计算的步骤要少一些，这证明对于同一个题，牛顿下山迭代法的优越性比一般的牛顿迭代要高。实验三解线性代数方程组——列主元消元法实验内容：用列主元法解线性方程组= 实验过程：利用C语言编程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){doublet,m1,m2,m3,x1,x2,x3,x4,a[4][5]={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};inti,j;for(i=0;i<4;i++) {for(j=0;j<5;j++)printf("%lf",a[i][j]);printf("\n");}for(i=0;i<4;i++) for(j=0;j<5;j++) if(a[i][0]<a[i+1][0]) t=a[i][0]; a[i][0]=a[i+1][0];a[i+1][0]=t;for(i=1;i<4;i++){m1=a[i][0]/a[0][0];for(j=0;j<5;j++) a[i][j]=a[i][j]-m1*a[0][j];}printf("第一次消元a[][]=\n");for(i=0;i<4;i++) {for(j=0;j<5;j++) printf("%lf",a[i][j]);printf("\n");}printf("\n");for(i=1;i<4;i++) for(j=0;j<5;j++) if(a[i][0]<a[i+1][0]) t=a[i][0]; a[i][0]=a[i+1][0];a[i+1][0]=t;for(i=2;i<4;i++){m2=a[i][1]/a[1][1];for(j=0;j<5;j++) a[i][j]=a[i][j]-m2*a[1][j];}printf("第二次消元a[][]=\n");for(i=0;i<4;i++) {for(j=0;j<5;j++) printf("%lf",a[i][j]);printf("\n");}printf("\n"); for(i=0;i<4;i++)for(j=0;j<5;j++) if(a[i][0]<a[i+1][0]) t=a[i][0]; a[i][0]=a[i+1][0];a[i+1][0]=t;for(i=3;i<4;i++){m3=a[i][2]/a[2][2];for(j=0;j<5;j++) a[i][j]=a[i][j]-m3*a[2][j];}printf("第三次消元a[][]=\n");for(i=0;i<4;i++) {for(j=0;j<5;j++) printf("%lf",a[i][j]); printf("\n");} printf("\n"); x4=a[3][4]/a[3][3]; x3=(a[2][4]-x4*a[2][3])/a[2][2]; x2=(a[1][4]-x4*a[1][3]-x3*a[1][2])/a[1][1]; x1=(a[0][4]-x4*a[0][3]-x3*a[0][2]-x2*a[0][1])/a[0][0]; printf("x1=%lf\nx2=%lf\nx3=%lf\nx4=%lf",x1,x2,x3,x4);}2.命令窗口结果截屏如下：实验结果分析：从程序中看，似乎程序很复杂，其实不然，中间有重复的语句，我这种方法虽然看起来比较过于冗长，但便于读者阅读，通俗易懂，我也曾想用循环体来代替中间的重复语句，但发现我编的循环体与程序中其他的循环体冲突，使程序无法运行出结果，怎么也理不清。既然计算机的运算能力很强，这点冗余度应该不是问题。实验四解线性代数方程组——追赶法实验内容：，用追赶法解三角方程组Ax=d，A=，d=实验过程：利用C语言编程求解，程序如下:#include<stdio.h>#include<math.h>voidmain(){inti;doublea[5]={0,-1,-1,-1,-1},b[5]={2,1,1,1,1},c[4]={2,2,2,2},f[5]={6,7,9,11,1};doubler[4],y[5],x[5];for(r[0]=c[0]/b[0],i=1;i<=3;i++)r[i]=c[i]/(b[i]-a[i]*r[i-1]);for(y[0]=f[0]/b[0],i=1;i<=4;i++)y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*r[i-1]);for(x[4]=y[4],i=3;i>=0;i--)x[i]=y[i]-r[i]*x[i+1];for(i=0;i<=4;i++)printf("%lf\n",x[i]);}2.命令窗口结果截屏如下:实验结果分析：运用追赶法与杜立特法有些类似，期间分解矩阵要运用到多个循环体，还有镶嵌式的循环体，这样很容易造成程序内冲突，于是运用另外的方法，定义的数组不代表一行，而是代表矩阵对角线上有元素的对角斜线。这样运算的过程就比较清晰了。实验五插值与拟合——三次样条差值实验内容：已知直升机旋转机翼外形曲线轮廓线上的某些点及端点处的一阶导数值，是S0=1.8648,S19=-0.046115试计算该曲线的横坐标为2,4,6,12,16,30,60,110,180,280,400,515该点处的纵坐标（要求该曲线具有二阶光滑度）。k01234560.523.18.017.9528.6539.6250.655.287949.413.8420.224.928.4431.1k78910111278104.6156.6208.6260.7312.53536.536.634.631.631.0k131415161718364.4416.346849450752020.914.87.83.71.50.2实验过程：1.利用C语言编程求解，程序如下：#include<stdio.h>#include<math.h>voidmain(){inti;doublex[19]={0.52,3.1,8.0,17.95,28.65,39.62,50.65,78.0,104.6,156.6,208.6,260.7,312.5,364.4,416.3,468.0,494.0,507.0,520.0},y[19]={5.28794,9.4,13.84,20.2,24.9,28.44,31.1,35.0,36.5,36.6,34.6,31.6,31.0,20.9,14.8,7.8,3.7,1.5,0.2};doubleh[18],u[18],v[18],g[18],a[19],b[19],c[19];for(i=0;i<19;i++){h[i]=x[i+1]-x[i];u[i]=h[i]/(h[i]+h[i+1]); v[i]=1-u[i]; g[i+1]=(6/(h[i+1]+h[i+2]))*((y[i+2]-y[i+1])/h[i+2]-(y[i+1]-y[i])/h[i+1]);}g[0]=6/(h[1])*((y[1]-y[0])/h[1]-1.8648);g[18]=6/(h[19])*(-0.046115-(y[18]-y[17])/h[17]);for(i=0;i<19;i++)printf("x[%d]:%lf",i,x[i]);printf("y[%d]:%lf",i,y[i]);for(i=0;i<19;i++)b[i]=2;a[0]=0,a[18]=1;for(i=1;i<=17;i++)a[i]=u[i-1];c[0]=1;for(i=1;i<=17;i++)c[i]=v[i-1];doubleB[18],Y[19],M[19];for(B[0]=c[0]/b[0],i=1;i<=17;i++)B[i]=c[i]/(b[i]-a[i]*B[i-1]);for(Y[0]=g[0]/b[0],i=1;