数据插值的matlab实现.doc

实验8 插值实验目的：1了解插值的基本原理； 2熟练掌握插值的Matlab实现；3会用插值方法解决一些实际问题。实验内容：在数学中，有时需要查表，如对数表。在具体查表时，需要的数据表中可能没有，这时一般可以先找出它相邻的数，再从表中查出其相应结果，然后按一定的关系把这些相邻的数以及它相应的结果加以修正，就可求出要查数的数据结果的近似值.这个修正关系就是一种插值.在实践中，常常需要测量某些数据，但由于客观条件的限制，所测得的数据可能不够细密，满足不了实践的需要，这时便可以通过插值方法对数据进行加密处理。此外，对于给定的离散数据对，如果要找一个函数来近似描述其对应关系，常常也需要插值。8.1 插值简介微积分中，用函数描述一条平面曲线，但在实际问题中，此函数关系往往通过实验观测得到的一组数据给出，如表8-1。表8-1 怎样根据表8-1中的数据，找出自变量与因变量之间的关系的表达式？可以采用插值。所谓插值，就是根据给定的数据表，寻找一个解析形式的函数，来近似地代替实际关系，并且满足。对于解析函数，可选择不同的类型，最常用的是代数多项式。不同的插值问题，构造插值多项式的方法有时也不相同，较为常用的有拉格朗日（Lagrange）插值、分段线性插值和样条插值。8.2 一维插值的MATLAB实现 MATLAB中常用的一维插值指令为interp1(x,y,cx,method) 其中，x，y分别表示数据点的横坐标和纵坐标向量；cx为待插值的横坐标数据向量，它不能超出x的范围；method为插值的使用方法，包括：nearest（最近邻点插值）、linear（线性插值）、spline（三次样条插值）和cubic（三次插值）。 例1 试用分段线性插值和三次样条插值分别对函数（）进行插值，并与实际曲线图比较。解：假设已知给定曲线的初始节点为-5，且步长为1，则用分段线性插值和三次样条插值分别对其插值作图，相应程序如下：x0=-5:5;y0=1./(1+x0.2); %产生节点（x0,y0)x=-5:0.1:5;y=1./(1+x.2); %产生插值点x y2=interp1(x0,y0,x);figure(2),plot(x,y,b,x,y2,k:),grid %使用分段线性插值并作图y3=interp1(x0,y0,x,spline);figure(3),plot(x,y,b,x,y3,k:),grid %使用三次样条插值并作图运行得图形见图8.1和图8.2。 图8.1 分段线性插值与实际曲线比较图 图8.2 三次样条插值与实际曲线比较图由图可见，分段线性插值与三次样条插值效果较好，而后者更好，其插值曲线与所给曲线几乎重合。此外需要注意的是，分段线性插值函数在节点处不可导，所以在实际应用时如果光滑程度要求较高，则应采用三次样条插值。例2 在一天24小时内,从零点开始每间隔2小时测得的环境温度为(摄氏度)12,9,9,10,18,24,28,27,25,20,18,15,13,推测在每一秒时的温度.并利用不同的插值方法描绘温度曲线键入：x=0:2:24;y=12 9 9 10 18 24 28 27 25 20 18 15 13;xi=0:1/3600:24;yi=interp1(x,y,xi,nearest);hold onplot(xi,yi,r); yi=interp1(x,y,xi,linear);plot(xi,yi,g);yi=interp1(x,y,xi,spline); plot(xi,yi,b);yi=interp1(x,y,xi,cubic); plot(xi,yi,y);8.3 二维插值的MATLAB实现对二维插值，MATLAB分别给出了针对插值基点为网格点的插值函数interp2和针对插值基点为散乱节点的插值函数griddata，其格式分别为（1）N维插值函数interpN（x,y,z,cx,cy, method）其中N可以为2,3,如N=2为二维插值,调用格式为 zi=interp2(x,y,z,Xi,Yi,method)其中x,y为横纵坐标上的坐标点,(x,y)=mashgrid(x,y)生成平面网格点,z为观测到的在网格点上的二元函数值.(x,y,z)构成空间插值节点.引入两个向量xi,yi.xi为横坐标上的插值点,yi为纵坐标上的插值点.便可给出Xi,Yi=meshgrid(xi,yi)zi为新的或者是加细了的网格点上产生的插值结果(函数值). method 表示采用的插值方法.nearest 最邻近插值,linear线性插值,cubic双三次插值.缺省时表示线性插值.所有的插值方法都要求x和y是单调的网格,x和y可以是等距的也可以是不等距的。（2）griddata（x,y,z,cx,cy, method）其中，x，y，z分别为数据点的横坐标、纵坐标和竖坐标向量，要求x和y的分量值必须单调；cx，cy分别为待插值点的横坐标和纵坐标的数据向量，它不能超出x的范围；method为插值的使用方法，包括：nearest（最近邻点插值）、linear（线性插值）、v4（MATLAB中所提供的插值方法）和cubic（三次插值）。例3 在某山区（平面区域内，单位：m）测得一些点的高度（单位：m）见表8-2，试作出该山区的地貌图。表8-2 某山区一些点的高度表 120016002000240028003200360040004400480012001130125012801230104090050070078075016001320145014201400130070090085084038020001390150015001400900110010609508709002400150012001100135014501200115010108801000280015001200110015501600155013801070900105032001500155016001550160016001600155015001500360014801500155015101430130012009808507504000145014701320128012001080940780620460440014301440114011101050950820690540380480014001410960940880800690570430290解：编写M文件shanqu.m如下：x=1200:400:4800;y=1200:400:4800;z=1130 1250 1280 1230 1040 900 500 700 780 750;13201450 1420 1400 1300 700 900 850 840 380;1390 1500 1500 1400 900 1100 1060 950 870 900;15001200 1100 1350 1450 1200 1150 1010 880 1000;1500 1200 1100 1550 1600 1550 1380 1070 900 1050;1500 1550 1600 1550 1600 1600 1600 1550 1500 1500;1480 1500 1550 1510 1430 1300 1200 980 850 750;1450 1470 1320 1280 1200 1080 940 780 620 460;1430 1440 1140 1110 1050 950 820 690 540 380;1400 1410 960 940 880 800 690 570 430 290;figure(1),mesh(x,y,z), xlabel(X),ylabel(Y),zlabel(Z) %用原始数据作山区表面图x1,y1=meshgrid(1200:50:4800,1200:50:4800); z1=interp2(x,y,z,x1,y1,spline);figure(2),mesh(x1,y1,z1), xlabel(X),ylabel(Y),zlabel(Z) %将数据加密，作山区表面图运行得图形见图8.3和图8.4。 图8.3原始数据山区地貌图 图8.4插值后山区地貌图例4 在某水道（平面区域75200,-90150单位：m）测得一些点的深度，数据见表8-3，已知某船只的吃水线为5米，试画出该水道的海底地貌图及船的禁入区。表8-3 某水道一些点的水深数据129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.57.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5 4 8 6 8 6 8 8 9 9 8 8 9 4 9解： 在MATLAB中编写M文件haiyum如下：x=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5;y=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5;z=-4 8 6 8 6 8 8 9 9 8 8 9 4 9;cx,cy=meshgrid(75:5:200,-90:5:150);cz=griddata(x,y,z,cx,cy,cubic);figure(1),mesh(cx,cy,cz);view(-60,30);figure(2), contour(cx,cy,cz,-5,-5,k) %绘制等高线运行程序得图形见图8.5和图8.6，其中图8.6中实线圈定区域为船只禁入区。 图8.5 某水道的海底地貌图 图8.6 船只禁入区域图例5 轮船的甲板成近似半椭圆面形，为了得到甲板的面积，首先测得横向最大相间为8.534米，然后等间距地测得纵向高度，自左向右分别为：0.914，5.060，7.772，8.717，9.083，9.144，9.083，8.992，8.687，7.376，2.073。计算甲板的面积。解：横向最大相间为8.534米，然后等间距地测得纵向高度,共有11个值，所以应该是把8.534米分成12分，对应的值为纵向高度；以左边零点位坐标原点，建立坐标系。线性插值得到图形(如图8.7)，再用数值积分可求面积x=linspace(0,8.534,13); y=0 0.914 5.060 7.772 8.717 9.083 9.144 9.083 8.992 8.687 7.376 2.073 0;x0=0:0.001:8.534;y1=interp1(x,y,x0); %线性插值得到x0处对应纵坐标的值figure,plot(x,y,k*,x0,y1,-r);S=trapz(y1)*0.001 %进行数值积分求甲板的面积运行程序得到S =54.6894图8.7 甲板横向坐标对应的高度值8.4 差值的建模示例：零件加工问题问题：待加工零件的外形根据手艺要求由一组数据给出（在平面情况下），用数控机床加工时刀具必须沿这些数据点前进，由于刀具每次只能沿方向或方向走非常小的一步，所以需要将已知数据加密，得到加工所要求的步长很小的坐标已知某机翼断面的下轮廓线上的部分数据，见表8-4。现在需要得到坐标每改变0.1时的坐标试完成加工所需数据，画出曲线，并求1315范围内的最小值表8-4 机翼断面下轮廓线上的部分数据 0 3 5 7 9 11 12 13 14 15 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.61 问题的分析解决零件加工问题的关键在于计算加密点的纵坐标，此坐标显然可以通过插值得到，而的最小值，则可以通过编程逐点搜索获得。2 模型的建立及其求解根据零件的加工要求，以所给数据为插值节点，以0. 1为步长，在=0到=15范围内作三次样条插值，相应的MATLAB指令为：x0=0,3,5,7,9,11,12,13,14,15;y0=0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6;x=0:0.1:15;y=interp1(x0,y0,x,spline);plot(x0,y0, k+,x,y, r),grid运行可得加密点的纵坐标，从而得下轮廓插值曲线见图8.8。 图8.8 某机翼断面的下轮廓插值曲线下面指令可在1315范围内搜索的最小值。for i=1:length(x) if x(i)=13 j=i; end,end %搜索数组x中13所对应的元素的下标 ymin=y(j);for k=j+1:length(y) %逐点搜索数组y的最小值 if y(k)ymin ymin=y(k);m=k; end,endymin,x(m)运行得在13,15上的最小值为0.9828