数据插值的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 on
plot(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 某山区一些点的高度表

1200
1600
2000
2400
2800
3200
3600
4000
4400
4800
1200
1130
1250
1280
1230
1040
900
500
700
780
750
1600
1320
1450
1420
1400
1300
700
900
850
840
380
2000
1390
1500
1500
1400
900
1100
1060
950
870
900
2400
1500
1200
1100
1350
1450
1200
1150
1010
880
1000
2800
1500
1200
1100
1550
1600
1550
1380
1070
900
1050
3200
1500
1550
1600
1550
1600
1600
1600
1550
1500
1500
3600
1480
1500
1550
1510
1430
1300
1200
980
850
750
4000
1450
1470
1320
1280
1200
1080
940
780
620
460
4400
1430
1440
1140
1110
1050
950
820
690
540
380
4800
1400
1410
960
940
880
800
690
570
430
290
解：编写M文件shanqu。m如下：
x=1200：400:4800;
y=1200：400：4800;
z=［1130 1250 1280 1230 1040 900 500 700 780 750；
1320 1450 1420 1400 1300 700 900 850 840 380；
1390 1500 1500 1400 900 1100 1060 950 870 900;
1500 1200 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。5

7。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.6
1． 问题的分析
解决零件加工问题的关键在于计算加密点的纵坐标,此坐标显然可以通过插值得到，而的最小值，则可以通过编程逐点搜索获得。
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，end
ymin,x（m)
运行得在［13,15]上的最小值为0.9828，相应的=13。8．
实验任务：
1．已知的部分函数表数据如下:

0
0.2
0.4
0.6
0.8
1.0
1。2
1。4
1.6

0
0。1987
0.3894
0.5646
0.7174
0.8415
0。9320
0.9854
0.9996
试根据表中的数据，分别利用分段线性插值和三次样条插值求.
2．在一丘陵地带测量高程，和方向每隔100m测一个点，得高度见表8—5，试作出该丘陵地带的地形图,并求出该地带的最高点和该点的高度.
表8—5 某丘陵测量数据表

100
200
300
400
100
636
697
624
478
200
698
712
630
478
300
680
674
598
412
400
662
626
552
334