北航数值分析B第二题报告完整版

1、数学分析B 计算实习题目第二题 1 算法设计方案1. 总体思路：由于需要四种方法进行细化，所以大体需要四个模块，实际上采用在main函数里面进行声明和调用几个子函数，f1(),f2(),f3(),分别是线性插值，二次插值和三次插值，prepare()函数用于求关于Mi的线性方程的参数，chase()函数用于采用追赶法求出Mi，three_moment()函数用于求出三弯矩法的y表达式，求出各种方法下的y的表达式后，进而x就容易求了，利用final_data()函数将数据输出，在matlab下画出细化后的轮廓线。2. 线性插值：插值法均采用牛顿法，d用于存放一阶差商值，这个值不仅在线性插值中要用

2、，在二次插值与三次插值中都需要它进行计算，同理，二次插值中的二阶差商也要用作三次插值中的计算。线性插值算出一阶差商后基本就可以代入公式得到y的表达式了。3. 二次插值：每三个点为一个插值区域，共六个插值区域，利用差商表计算得出二阶差商，展开牛顿二次插值多项式，得出结果。4. 三次插值：每四个点为一个插值区域，共四个插值区域，这时计算量就已经有点大了，为简化运算，采用不单独计算，而是在二次插值的基础上进行计算，利用二次插值的已有数据，但是有两个二次插值在二次插值函数中并没有得到计算，即f(x3,x4,x5)和f(x7,x8,x9)，所以在二次插值函数中要添加相关代码。5. 样条插值：计算最为复杂

3、，但其实并不比三次插值容易错，因为只是公式很多，所以化为三个函数来完成，每个函数负责一部分计算。6. 输出数据：前面的计算只是算出公式，对于不同的插值法，y的分段表达式不一样，采用if，else if语句以及for循环并用到了pow()函数算幂，完成对表中数据x从1到190的细化。2 全部源程序#include<stdio.h>#include<math.h>double x13=0.00,4.74,9.50,19.00,38.00,57.00,76.00,95.00,114.00,133.00,152.00,171.00,190.00;double y13=0.00,

4、5.32,8.10,11.97,16.15,17.10,16.34,14.63,12.16,9.69,7.03,3.99,0.00;double d12=0;/存放一阶差商double da12=0;/存放线性插值常数项double dt8=0;/存放二阶差商double dta8=0;/存放二阶插值x项系数double dtb8=0;/存放二阶插值常数项double ddt4=0;/存放二阶差商用作计算三阶差商double dth4=0;/存放三阶差商double dtha4=0;/存放三阶插值附加x2项系数double dthb4=0;/存放三阶插值附加x项系数double dthc4=0

5、;/存放三阶插值附加常数项double h13=0;/存放样条插值的一些参数double a13=0;double b13=0;double c13=0;double e13=0;double p13=0;double q13=0;double z13=0;double m13=0;double ta12=0;double tb12=0;double tc12=0;double td12=0;double ans1190=0;double ans2190=0;double ans3190=0;double ans4190=0;void main()/调用几个子程序void f1(); /线性插

6、值void f2();/二次插值void f3();/三次插值void prepare();/样条插值计算关于Mn的线性方程组的有关参数void chase();/追赶法解线性方程组void three_moment();/三弯矩法void final_data();f1();f2();f3();prepare();chase();three_moment(); final_data(); void f1()/线性插值函数int i=0;printf("分段线性插值法：nn");for(i=0;i<12;i+)di=(yi+1-yi)/(xi+1-xi);dai=yi

7、-di*xi;printf("插值区间：(%6.2f,%6.2f)t ",xi,xi+1);printf("插值函数：P1(x)=%6.2fx+%6.2fn ",di,dai);for(i=0;i<=190;i+)/计算具体数据if(i<5)ans1i=d0*i+da0;else if(i<10)ans1i=d1*i+da1;else if(i<20)ans1i=d2*i+da2;else if(i<39)ans1i=d3*i+da3;else if(i<58)ans1i=d4*i+da4;else if(i<7

8、7)ans1i=d5*i+da5;else if(i<96)ans1i=d6*i+da6;else if(i<115)ans1i=d7*i+da7;else if(i<134)ans1i=d8*i+da8;else if(i<153)ans1i=d9*i+da9;else if(i<172)ans1i=d10*i+da10;else if(i<191)ans1i=d11*i+da11;void f2()/二次插值int i=0;printf("nn");printf("二次插值法;nn");for(i=0;i<6

9、;i+)dti=(d2*i+1-d2*i)/(x2*i+2-x2*i);dtai=d2*i-dti*(x2*i+x2*i+1);dtbi=y2*i-d2*i*x2*i+dti*x2*i*x2*i+1; printf("插值区间：(%6.2f,%6.2f) ",x2*i,x2*i+2);printf("插值函数：P1(x)=(%6.5fx2)+(%6.3fx)+(%.2f)n ",dti,dtai,dtbi); dt6=(d4-d3)/(x5-x3);dta6=d3-dt6*(x3+x4);dtb6=y3-d3*x3+dt6*x3*x4;/计算以x3x4x

10、5插值的函数，为三次插值计算做准备dt7=(d10-d9)/(x11-x9);dta7=d9-dt7*(x9+x10);dtb7=y9-d9*x9+dt7*x9*x10;/计算以x7x8x9插值的函数，为三次插值计算做准备for(i=0;i<=190;i+)/计算具体数据if(i<10)ans2i=dt0*pow(i,2)+dta0*i+dtb0;else if(i<39)ans2i=dt1*pow(i,2)+dta1*i+dtb1;else if(i<77)ans2i=dt2*pow(i,2)+dta2*i+dtb2;else if(i<115)ans2i=dt

11、3*pow(i,2)+dta3*i+dtb3;else if(i<153)ans2i=dt4*pow(i,2)+dta4*i+dtb4;else if(i<191)ans2i=dt5*pow(i,2)+dta5*i+dtb5;void f3()/三次插值int i=0;printf("nn"); printf("三次插值法;nn");ddt0=(d2-d1)/(x3-x1);ddt1=(d4-d3)/(x5-x3);ddt2=(d8-d7)/(x9-x7);ddt3=(d10-d9)/(x11-x9);dth0=(ddt0-dt0)/(x3-

12、x0);dth1=(dt2-ddt1)/(x6-x3);dth2=(ddt2-dt3)/(x9-x6);dth3=(dt5-ddt3)/(x12-x9);for(i=0;i<4;i+)dthai=(x3*i+x3*i+1+x3*i+2)*dthi;dthbi=(x3*i*x3*i+1+x3*i*x3*i+2+x3*i+1*x3*i+2)*dthi;dthci=(x3*i*x3*i+1*x3*i+2)*dthi; printf("插值区间：(%6.2f,%6.2f)n ",x0,x3);printf("插值函数：P1(x)=(%.6fx3)+(%.2fx2)+

13、(%.2fx)+(%.2f)n ",dth0,dt0-dtha0,dta0+dthb0,dtb0-dthc0); printf("插值区间：(%6.2f,%6.2f)n ",x3,x6); printf("插值函数：P1(x)=(%.6fx3)+(%.2fx2)+(%.2fx)+(%.2f)n ",dth1,dt6-dtha1,dta6+dthb1,dtb6-dthc1); printf("插值区间：(%6.2f,%6.2f)n ",x6,x9); printf("插值函数：P1(x)=(%.6fx3)+(%.2f

14、x2)+(%.2fx)+(%.2f)n ",dth2,dt3-dtha2,dta3+dthb2,dtb3-dthc2);printf("插值区间：(%6.2f,%6.2f)n ",x9,x12); printf("插值函数：P1(x)=(%.6fx3)+(%.2fx2)+(%.2fx)+(%.2f)n ",dth3,dt7-dtha3,dta7+dthb3,dtb7-dthc3);for(i=0;i<=190;i+)/计算具体数据if(i<20)ans3i=dth0*pow(i,3)+(dt0-dtha0)*pow(i,2)+(dt

15、a0+dthb0)*i+dtb0-dthc0;else if(i<77)ans3i=dth1*pow(i,3)+(dt6-dtha1)*pow(i,2)+(dta6+dthb1)*i+dtb6-dthc1;else if(i<134)ans3i=dth2*pow(i,3)+(dt3-dtha2)*pow(i,2)+(dta3+dthb2)*i+dtb3-dthc2;else if(i<191)ans3i=dth3*pow(i,3)+(dt7-dtha3)*pow(i,2)+(dta7+dthb3)*i+dtb7-dthc3;void prepare()/样条插值计算关于Mn的

16、线性方程组的有关参数int i=0;int j=0;for(i=1;i<13;i+) hi=xi-xi-1;for(j=1;j<=11;j+)aj=hj+1/(hj+hj+1);cj=1-aj;bj=(yj+1-yj)/hj+1-(yj-yj-1)/hj)*6/(hj+hj+1);/printf("a=%.3f,b=%.3f,c=%.3fn",aj,bj,cj);c12=0;b0=0;b12=0;a0=0;void chase()/追赶法解线性方程组(注意此模块中a对应追赶法中的c；c对应追赶法中的d)int i=0;int j=0;p0=2;z0=b0/p0;

17、for(i=0;i<=11;i+)qi=ai/pi;pi+1=2-ci+1*qi;zi+1=(bi+1-ci+1*zi)/pi+1; m12=z12;for(j=11;j>=0;j-) mj=zj-qj*mj+1; /printf("m=%.5f",mj);void three_moment()/三弯矩法int i=0;int j=0;printf("nn");printf("三次样条插值结果：");for(j=1;j<=12;j+)taj-1=yj-1/hj-mj-1*hj/6;tbj=yj/hj-mj*hj/6;

18、tcj-1=mj-1/(6*hj);tdj=mj/(6*hj);for(i=1;i<=12;i+)printf("插值区间：(%6.2f,%6.2f)n ",xi-1,xi);printf("插值函数s(x)=%.3f(%.2f-x)3+%.3f(x-%.2f)3+%.3f(%.2f-x)+%.3f(x-%.2f)n",tcj-1,xi,tdj,xi-1,tai-1,xi,tbi,xi-1);for(i=0;i<=190;i+)/计算具体数据if(i<5) ans4i=tc0*pow(x1-i,3)+td1*pow(i-x1,3)+ta

19、0*(x1-i)+tb1*(i-x0); else if(i<10)ans4i=tc1*pow(x2-i,3)+td2*pow(i-x1,3)+ta1*(x2-i)+tb2*(i-x1);else if(i<20)ans4i=tc2*pow(x3-i,3)+td3*pow(i-x2,3)+ta2*(x3-i)+tb3*(i-x2);else if(i<39)ans4i=tc3*pow(x4-i,3)+td4*pow(i-x3,3)+ta3*(x4-i)+tb4*(i-x3);else if(i<58)ans4i=tc4*pow(x5-i,3)+td5*pow(i-x4,

20、3)+ta4*(x5-i)+tb5*(i-x4);else if(i<77) ans4i=tc5*pow(x6-i,3)+td6*pow(i-x5,3)+ta5*(x6-i)+tb6*(i-x5);else if(i<96)ans4i=tc6*pow(x7-i,3)+td7*pow(i-x6,3)+ta6*(x7-i)+tb7*(i-x6);else if(i<115)ans4i=tc7*pow(x8-i,3)+td8*pow(i-x7,3)+ta7*(x8-i)+tb8*(i-x7);else if(i<134)ans4i=tc8*pow(x9-i,3)+td9*po

21、w(i-x8,3)+ta8*(x9-i)+tb9*(i-x8);else if(i<153)ans4i=tc9*pow(x10-i,3)+td10*pow(i-x9,3)+ta9*(x10-i)+tb10*(i-x9);else if(i<172)ans4i=tc10*pow(x11-i,3)+td11*pow(i-x10,3)+ta10*(x11-i)+tb11*(i-x10);else if(i<191)ans4i=tc11*pow(x12-i,3)+td12*pow(i-x11,3)+ta11*(x12-i)+tb12*(i-x11);void final_data()

22、/输出最终数据int j=0;printf("y值线性插值二次插值三次插值样条插值n");printf("x值n");for(j=0;j<=190;j+)printf("%d%.2f %.2f %.2f %.2f n",j,ans1j,ans2j,ans3j,ans4j);三全部数据分段线性插值法：插值区间：( 0.00, 4.74) 插值函数：P1(x)= 1.12x+ 0.00 插值区间：( 4.74, 9.50) 插值函数：P1(x)= 0.58x+ 2.55 插值区间：( 9.50, 19.00) 插值函数：P1(x)=

23、 0.41x+ 4.23 插值区间：( 19.00, 38.00) 插值函数：P1(x)= 0.22x+ 7.79 插值区间：( 38.00, 57.00) 插值函数：P1(x)= 0.05x+ 14.25 插值区间：( 57.00, 76.00) 插值函数：P1(x)= -0.04x+ 19.38 插值区间：( 76.00, 95.00) 插值函数：P1(x)= -0.09x+ 23.18 插值区间：( 95.00,114.00) 插值函数：P1(x)= -0.13x+ 26.98 插值区间：(114.00,133.00) 插值函数：P1(x)= -0.13x+ 26.98 插值区间：(13

24、3.00,152.00) 插值函数：P1(x)= -0.14x+ 28.31 插值区间：(152.00,171.00) 插值函数：P1(x)= -0.16x+ 31.35 插值区间：(171.00,190.00) 插值函数：P1(x)= -0.21x+ 39.90二次插值法;插值区间：( 0.00, 9.50) 插值函数：P1(x)=(-0.05667x2)+( 1.391x)+(0.00) 插值区间：( 9.50, 38.00) 插值函数：P1(x)=(-0.00657x2)+( 0.595x)+(3.04) 插值区间：( 38.00, 76.00) 插值函数：P1(x)=(-0.00237

25、x2)+( 0.275x)+(9.12) 插值区间：( 76.00,114.00) 插值函数：P1(x)=(-0.00105x2)+( 0.090x)+(15.58) 插值区间：(114.00,152.00) 插值函数：P1(x)=(-0.00026x2)+(-0.065x)+(22.99) 插值区间：(152.00,190.00) 插值函数：P1(x)=(-0.00132x2)+( 0.265x)+(-2.85)三次插值法;插值区间：( 0.00, 19.00) 插值函数：P1(x)=(0.002330x3)+(-0.09x2)+(1.50x)+(0.00) 插值区间：( 19.00, 76

26、.00) 插值函数：P1(x)=(0.000037x3)+(-0.01x2)+(0.62x)+(3.04) 插值区间：( 76.00,133.00) 插值函数：P1(x)=(0.000018x3)+(-0.01x2)+(0.58x)+(0.38) 插值区间：(133.00,190.00) 插值函数：P1(x)=(-0.000014x3)+(0.01x2)+(-0.94x)+(65.55)三次样条插值结果：插值区间：( 0.00, 4.74) 插值函数s(x)=0.000(4.74-x)3+1.122(x-0.00)3+0.000(4.74-x)+1.256(x-0.00)插值区间：( 4.74

27、, 9.50) 插值函数s(x)=0.000(9.50-x)3+1.122(x-4.74)3+1.252(9.50-x)+1.704(x-4.74)插值区间：( 9.50, 19.00) 插值函数s(x)=0.000(19.00-x)3+1.122(x-9.50)3+0.858(19.00-x)+1.286(x-9.50)插值区间：( 19.00, 38.00) 插值函数s(x)=0.000(38.00-x)3+1.122(x-19.00)3+0.682(38.00-x)+0.876(x-19.00)插值区间：( 38.00, 57.00) 插值函数s(x)=0.000(57.00-x)3+1

28、.122(x-38.00)3+0.876(57.00-x)+0.914(x-38.00)插值区间：( 57.00, 76.00) 插值函数s(x)=0.000(76.00-x)3+1.122(x-57.00)3+0.914(76.00-x)+0.867(x-57.00)插值区间：( 76.00, 95.00) 插值函数s(x)=0.000(95.00-x)3+1.122(x-76.00)3+0.867(95.00-x)+0.779(x-76.00)插值区间：( 95.00,114.00) 插值函数s(x)=0.000(114.00-x)3+1.122(x-95.00)3+0.779(114.0

29、0-x)+0.637(x-95.00)插值区间：(114.00,133.00) 插值函数s(x)=0.000(133.00-x)3+1.122(x-114.00)3+0.637(133.00-x)+0.513(x-114.00)插值区间：(133.00,152.00) 插值函数s(x)=0.000(152.00-x)3+1.122(x-133.00)3+0.513(152.00-x)+0.371(x-133.00)插值区间：(152.00,171.00) 插值函数s(x)=0.000(171.00-x)3+1.122(x-152.00)3+0.371(171.00-x)+0.222(x-152

30、.00)插值区间：(171.00,190.00) 插值函数s(x)=0.000(190.00-x)3+1.122(x-171.00)3+0.222(190.00-x)+0.000(x-171.00)y值 线性插值 二次插值 三次插值 样条插值x值0 0.00 0.00 0.00 0.001 1.12 1.33 1.41 1.572 2.24 2.56 2.65 2.633 3.37 3.66 3.74 3.804 4.49 4.66 4.70 5.035 5.47 5.54 5.52 5.546 6.06 6.31 6.24 6.277 6.64 6.96 6.87 6.898 7.22 7.

31、50 7.41 7.419 7.81 7.93 7.88 7.8810 8.30 8.33 8.30 8.3211 8.71 8.79 8.68 8.7612 9.12 9.23 9.04 9.2013 9.53 9.66 9.38 9.6214 9.93 10.08 9.73 10.0415 10.34 10.49 10.09 10.4516 10.75 10.88 10.48 10.8517 11.16 11.25 10.91 11.2418 11.56 11.62 11.41 11.6119 11.97 11.97 11.97 11.9720 12.19 12.31 12.30 12.3

32、121 12.41 12.63 12.61 12.6422 12.63 12.95 12.91 12.9523 12.85 13.24 13.19 13.2424 13.07 13.53 13.47 13.5225 13.29 13.80 13.73 13.7926 13.51 14.06 13.98 14.0427 13.73 14.31 14.22 14.2828 13.95 14.54 14.45 14.5129 14.17 14.76 14.67 14.7230 14.39 14.97 14.87 14.9231 14.61 15.16 15.07 15.1132 14.83 15.3

33、4 15.25 15.2933 15.05 15.51 15.43 15.4634 15.27 15.66 15.59 15.6235 15.49 15.81 15.74 15.7636 15.71 15.93 15.89 15.9037 15.93 16.05 16.02 16.0338 16.15 16.15 16.15 16.1539 16.20 16.24 16.27 16.2640 16.25 16.33 16.38 16.3741 16.30 16.41 16.48 16.4642 16.35 16.49 16.57 16.5543 16.40 16.57 16.65 16.634

34、4 16.45 16.63 16.73 16.7145 16.50 16.70 16.80 16.7746 16.55 16.76 16.86 16.8347 16.60 16.81 16.91 16.8948 16.65 16.86 16.96 16.9349 16.70 16.91 17.00 16.9850 16.75 16.95 17.03 17.0151 16.80 16.98 17.06 17.0452 16.85 17.02 17.08 17.0653 16.90 17.04 17.09 17.0854 16.95 17.06 17.10 17.0955 17.00 17.08

35、17.11 17.1056 17.05 17.09 17.11 17.1057 17.10 17.10 17.10 17.1058 17.06 17.10 17.09 17.0959 17.02 17.10 17.07 17.0860 16.98 17.09 17.05 17.0761 16.94 17.08 17.03 17.0562 16.90 17.07 17.00 17.0263 16.86 17.04 16.97 17.0064 16.82 17.02 16.94 16.9665 16.78 16.99 16.90 16.9366 16.74 16.95 16.86 16.8967

36、16.70 16.91 16.82 16.8568 16.66 16.87 16.77 16.8069 16.62 16.82 16.72 16.7570 16.58 16.76 16.67 16.7071 16.54 16.71 16.62 16.6572 16.50 16.64 16.57 16.5973 16.46 16.57 16.51 16.5374 16.42 16.50 16.46 16.4775 16.38 16.42 16.40 16.4176 16.34 16.34 16.34 16.3477 16.25 16.27 16.28 16.2778 16.16 16.20 16

37、.22 16.2079 16.07 16.12 16.15 16.1380 15.98 16.04 16.08 16.0581 15.89 15.96 16.01 15.9782 15.80 15.88 15.93 15.9083 15.71 15.80 15.85 15.8184 15.62 15.71 15.76 15.7385 15.53 15.62 15.67 15.6486 15.44 15.53 15.58 15.5587 15.35 15.44 15.49 15.4688 15.26 15.35 15.39 15.3789 15.17 15.25 15.29 15.2790 15

38、.08 15.15 15.18 15.1791 14.99 15.05 15.08 15.0792 14.90 14.95 14.97 14.9693 14.81 14.85 14.86 14.8594 14.72 14.74 14.75 14.7495 14.63 14.63 14.63 14.6396 14.50 14.52 14.51 14.5197 14.37 14.41 14.39 14.3998 14.24 14.29 14.27 14.2799 14.11 14.17 14.15 14.15100 13.98 14.05 14.02 14.02101 13.85 13.93 13.90 13.90102 13.72 13.81 13.77 13.77103 13.59 13.68 13.64 13.64104 13.46 13.55 13.51 13.50105 13.33 13.42 13.38 13.37106 13.20 13.29 13.24