数学实验3-插值与数值积分-安振华-2012011837.docx

清华大学数学实验报告实验3：插值与数值积分化学工程系 分2 安振华 2012011837【实验目的】1掌握用MATLAB 计算拉格朗日、分段线性、三次样条三种插值的方法，改变节点的数目，以及能对三种插值结果进行初步分析；2熟练掌握用MATLAB 及梯形公式、辛普森公式计算数值积分；3通过实例学习用插值和数值积分解决实际问题。【实验内容】预备：编制计算拉格朗日插值的M文件(lagr.m)function y=lagr(x0,y0,x)n=length(x0); m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end分段线性插值函数（interp.m）function y=interp1(x0,y0,x) n=length(x0);m=length(x); for i=1:m s=0; for j=1:n if x(i)=x0(j) s=y0(j); end if x(i)x0(j) s=s+(x(i)-x0(j-1)/(x(j)-x(j-1)*(y0(j)-y0(j-1); break end end y(i)=s; end【实验三：习题10】表3.10给出的x,y数据位于机翼断面的轮廓线上，y1与y2分别对应轮廓的上下线，假设需要得到x坐标每改变0.1时的y坐标。试完成加工所需数据，画出曲线，求机翼剖面的面积。表3.10 机翼剖面轮廓线数据x035791112131415y101.82.22.73.03.12.92.52.01.6y201.21.72.02.12.01.81.21.01.6【分析】由于给定的数据点是有限的，要想确定更多的有效数据，就要运用插值方法。而加工断面的面积，则应通过数值积分分别求出上下轮廓线对x 轴围成的面积，然后作差求得。【解答】运用学习的三种插值方法逐一计算，具体如下：1、拉格朗日多项式插值MATLAB程序:clf,shg,x0=0,3,5,7,9,11,12,13,14,15;y10=0,1.8,2.2,2.7,3,3.1,2.9,2.5,2,1.6;y20=0,1.2,1.7,2,2.1,2,1.8,1.2,1,1.6; %按照表格3.10输入原始数据x=0:0.1:15; %按题目要求以0.1的间隔产生插值点y1=lagr(x0,y10,x);y2=lagr(x0,y20,x); %在x方向计算拉格朗日插值subplot(1,2,1),plot(x0,y10,k,x0,y20,r)subplot(1,2,2),plot(x,y1,b,x,y2,b)S=trapz(x,y1)-trapz(x,y2) %计算机翼剖面面积运行结果：给出机翼剖面面积 S=40.30442、分段线性插值MATLAB程序:clf,shg,x0=0,3,5,7,9,11,12,13,14,15;y10=0,1.8,2.2,2.7,3,3.1,2.9,2.5,2,1.6;y20=0,1.2,1.7,2,2.1,2,1.8,1.2,1,1.6;x=0:0.1:15;y1=interp1(x0,y10,x)y2=interp1(x0,y20,x) %在x方向计算分段线性插值并输出结果subplot(1,2,1),plot(x0,y10,k,x0,y20,r)subplo(1,2,2),plot(x,y1,b,x,y2,b) %分段线性插值作图S=trapz(x,y1)-trapz(x,y2) %梯形公式计算面积运行结果：给出机翼剖面面积 S=10.75003、三次样条插值MATLAB程序:clf,shg,x0=0,3,5,7,9,11,12,13,14,15;y10=0,1.8,2.2,2.7,3,3.1,2.9,2.5,2,1.6;y20=0,1.2,1.7,2,2.1,2,1.8,1.2,1,1.6;x=0:0.1:15;y1=spline(x0,y10,x)y2=spline(x0,y20,x) %在x方向计算三次样条插值并输出结果subplot(1,2,1),plot(x0,y10,k,x0,y20,r)subplot(1,2,2),plot(x,y1,b,x,y2,b) %三次样条插值作图S=trapz(x,y1)-trapz(x,y2)运行结果：给出机翼剖面面积 S=11.3444整合程序，将三种差值方法集成在同一程序，可以直观地比较各自得出的结果。整合后的MATLAB程序:x0=0 3 5 7 9 11 12 13 14 15;y10=0 1.8 2.2 2.7 3.0 3.1 2.9 2.5 2.0 1.6;y20=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6;x=0:0.1:15;x0=0 3 5 7 9 11 12 13 14 15;a1=spline(x0,y10,x);a2=spline(x0,y20,x);b1=lagr(x0,y10,x);b2=lagr(x0,y20,x);c1=interp1(x0,y10,x);c2=interp1(x0,y20,x);subplot(1,3,1),plot(x,b1,r,x,b2,b)title(拉格朗日差值);subplot(1,3,2),plot(x,c1,r,x,c2,b)title(分段线性插值);subplot(1,3,3),plot(x,a1,r,x,a2,b)title(三次样条插值);subplot(1,3,1),plot(x,b1,r,x,b2,b)title(拉格朗日插值);subplot(1,3,2),plot(x,c1,r,x,c2,b)title(分段线性插值);subplot(1,3,3),plot(x,a1,r,x,a2,b)title(三次样条插值);运行结果：【结果分析】由图可知拉格朗日插值法得出图形与实际不符，不予考虑，可见拉格朗日插值方法的收敛性无法得到保证；分段线性插值的大致形状与原图相近，但它得到的图线不平滑，应用到机翼这样的部件上也是无法使用的。而三次样条插值比较光滑，误差最小。由梯形公式得到加工剖面面积为11.3444。根据梯形公式和辛普森公式的应用特点（已知被积函数的情况下应用辛普森函数更方便，而对于此处插值得到的数据则用梯形公式更方便），故此处应用梯形公式。又由于在插值之后，步长较短，因此推测梯形公式也较为准确。综上，此问题的最佳解决方案是采用三次样条插值方法，机翼剖面的面积为11.3444。 【实验三：习题11】图3.11是欧洲一个国家的地图，为算出它的国土面积，首先对该国地图做出如下测量：以由西向东方向为x轴，由南到北方向为y轴，选择方便的原点，并将从最西边界点到最东边界点在x轴上的区间适当的划分为若干段，在每个分点的y方向测出南边界点和北边界点的y坐标y1和y2，得到表3.11（单位：mm）。 根据地图的比例尺我们知道18mm相当于40km，试由测量数据计算该国国土的近似面积，与它的精确值41288km2作比较。图3.11表3.11 地图边界点数据x7.010.513.017.534.040.544.548.056.061.068.576.580.591.0y14445475050383030343634414546y24459707293100110110110117118116118118x96.0101.0104.0106.5111.5118.0123.5136.5142.0146.0150.0157.0158.0y143373328326555545250666668y2121124121121121122116838182868568【分析】由3.10解题思路知，可以先用三种插值计算方法分别进行插值，得到图像后分析选择与实际符合最好的方案进行计算。计算时可以选择辛普森积分或梯形积分方式进行积分运算。将结果与精确值进行比较和误差分析即可【解答】运用学习的三种插值方法分别计算，具体如下：clf,shg,x0=7,10.5,13,17.5,34,40.5,44.5,48,56,61,68.5,76.5,80.5,91,96,101,104,106.5,111.5,118,123.5,136.5,142,146,150,157,158;y10=44,45,47,50,50,38,30,30,34,36,34,41,45,46,43,37,33,28,32,65,55,54,52,50,66,66,68;y20=44,59,70,72,93,100,110,110,110,117,118,116,118,118,121,124,121,121,121,122,116,83,81,82,86,85,68;x=7:1:158;a1=lagr(x0,y10,x);a2=lagr(x0,y20,x);b1=interp1(x0,y10,x);b2=interp1(x0,y20,x);c1=spline(x0,y10,x);c2=spline(x0,y20,x);subplot(2,2,1),plot(x0,y10,r,x0,y20,r)title(A.实际地图 )subplot(2,2,4),plot(x,a1,m,x,a2,m)title(B.拉格朗日插值 );subplot(2,2,3),plot(x,b1,g,x,b2,g)title(C.分段线性插值);subplot(2,2,2),plot(x,c1,b,x,c2,b)title(D.三次样条插值);MATLAB程序:运行结果：由以上图像分析可得，拉格朗日插值所得结果与实际偏差非常大，顾舍弃该方法。选择三次样条插值或分线段性插值法分析计算。MATLAB程序:clf,shg,x0=7,10.5,13,17.5,34,40.5,44.5,48,56,61,68.5,76.5,80.5,91,96,101,104,106.5,111.5,118,123.5,136.5,142,146,150,157,158;y10=44,45,47,50,50,38,30,30,34,36,34,41,45,46,43,37,33,28,32,65,55,54,52,50,66,66,68;y20=44,59,70,72,93,100,110,110,110,117,118,116,118,118,121,124,121,121,121,122,116,83,81,82,86,85,68;x=7:1:158;b1=interp1(x0,y10,x);b2=interp1(x0,y20,x);c1=spline(x0,y10,x);c2=spline(x0,y20,x);S1=(trapz(x,b2)-trapz(x,b1)*(40/18)2 %分段线性插值梯形积分S2=(trapz(x,c2)-trapz(x,c1)*(40/18)2 %三次样条插值梯形积分运行结果：S1 = 42407S2 = 42458【结果分析】 这两种方法得到的最终结果都很相近，和精确值41288km2虽然也有一定误差，但并不是很大，分别为2.71%和2.83%，在误差可以接受的范围内。此题中，分段线性插值法的拟合效果最好，与实际轮廓最为相符。这与通常三次样条插值最适合的案例不符，原因可能是该题有关国家边界线，很多国界线沿经纬线分布，因而并不是严格的平滑曲线。 【实验三：习题12】在桥梁的一端每隔一段时间记录1min有几辆车过桥，得到下表的过桥车辆数量：时间车辆数/辆时间车辆数/辆时间车辆数/辆0:0029:001218:00222:00210:30519:00104:00011:301020:0095:00212:301221:00116:00514:00722:0087:00816:00923:0098:002517:002824:003试估计一天通过桥梁的车流量。【分析】表中所给21组数据（时间x0,车辆数y0）可视为节点，每间隔1min做一个插值点。根据前两题的经验，一般情况下三次样条插值法拟合较好，误差较小故此次直接可以通过三次样条插值法进行插值，得到曲线后求算面积即是一天通过桥梁的车流量。【解答】通过三次样条插值法进行插值，得到曲线后利用梯形积分法求算面积即是一天通过桥梁的车流量。MATLAB程序：x0=0 2 4 5 6 7 8 9 10.5 11.5 12.5 14 16 17 18 19 20 21 22 23 24;y0=2 2 0 2 5 8 25 12 5 10 12 7 9 28 22 10 9 11 8 9 3;x=0:1/60:24;q=spline(x0,y0,x);plot(x,q)title(一天内车流变化曲线)Q=60*trapz(x,q)运行结果：Q=12668【结果分析】作图如上图，可以看出图线光滑，拟合较好，综上所述，一天通过桥梁的车流量为：12668辆 【实验总结】