如何利用MATLAB进行数据插值?

文章目录

  • 前言
  • 1 引例-零件加工问题
  • 2 数据插值的计算机制
  • 3 数据插值的实现方法
  • 3 应用案例1-粮储仓的通风控制问题
  • 4 应用案例2-机动车刹车距离问题
  • 5 应用案例3-沙盘制作问题
  • 总结

前言

本文是科学计算与MATLAB语言课程的第5章第3、4小结的学习笔记,通过查阅本文,可以轻松掌握利用MATLAB进行数据插值。Enjoyyourreading!Enjoy\ your\ reading!Enjoy your reading!
欢迎大家👍,收藏⭐,转发🚀,
如有问题、建议请您在评论区留言💬。

1 引例-零件加工问题

在飞机制造中,机翼的加工是一项关键技术。由于机翼尺寸很大,通常在图纸中只能标出一些关键点的数据。下表给出了某型飞机机翼的下缘轮廓线数据,求x每改变0.1时y的值。

x 0 3 5 7 9 11 12 13 14 15
y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6
x=[0,3,5,7,9,11,12,13,14,15];
y=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
x1=0:0.1:15;
y1=interp1(x,y,x1,'spline');
plot(x1,y1)

fig1

2 数据插值的计算机制

从数学上来说,数据插值是一种函数逼近的方法。

xxx x1x2…xk…xnx_1 \ x_2 \ …x_k…x_nx1 x2 ...xk...xn
yyy y1y2…yk…yny_1\ y_2\ …y_k…y_ny1 y2 ...yk...yn

y=f(x)y=f(x)y=f(x)
它的实质就是用一个近似函数ϕ(x)\phi(x)ϕ(x)来逼近未知函数f(x)f(x)f(x),然后利用这个近似函数ϕ(x)\phi(x)ϕ(x)进行插值。

3 数据插值的实现方法

数据插值的实现方法
在MATLAB中,一维插值函数为interp1(),其调用格式为:
Y1=interp1(X,Y,X1,method)
该语句将根据X、Y的值,计算函数在X1处的值。其中,X、Y是两个等长的已知向量,分别表示采样点和采样值。X1是一个向量或标量,表示要插值的点。
method参数用于指定插值方法,常用的取值有以下四种:
(1)linear:线性插值,默认方法。将与插值点靠近的两个数据点用直线连接,然后在直线上选取对应插值点的数据。
(2)nearest:最近点插值。选择最近样本点的值作为插值数据。
(3)pchip:分段3次埃尔米特插值。采用分段三次多项式,除满足插值条件,还需满足在若干节点处相邻段插值函数的一阶导数相等,使得曲线光滑的同时,还具有保形性。
(4)spline:3次样条插值。每个分段内构造一个三次多项式,使其插值函数除满足插值条件外,还要求在各节点处具有连续的一阶和二阶导数。

思考:为什么这两种插值方法都用3次多项式而不用更高次的?
这里就要提一下龙格现象了,龙格(Runge)发现多项式插值并非次数越高越精确!
四种方法的比较:

x=[0,3,5,7,9,11,12,13,14,15]; 
y=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6]; 
x1=0:0.1:15; 
y1=interp1(x,y,x1,'spline'); %3次样条插值
subplot(2,2,1)
plot(x1,y1)
legend('3次样条插值')
hold on
y2=interp1(x,y,x1,'linear'); %线性插值
subplot(2,2,2)
plot(x1,y2,'r')
legend('线性插值')
y3=interp1(x,y,x1,'pchip'); %分段3次埃尔米特插值
subplot(2,2,3)
plot(x1,y3,'g')
legend('分段3次埃尔米特插值')
y4=interp1(x,y,x1,'nearest'); %最近点插值
subplot(2,2,4)
plot(x1,y4,'b')
legend('最近点插值')

fig2
线性插值和最近点插值方法比较简单。其中线性插值方法的计算量与样本点n无关。n越大,误差越小。
3次埃尔米特插值和3次样条插值都能保证曲线的光滑性。相比较而言,3次埃尔米特插值具有保形性;而3次样条插值要求其二阶导数也连续,所以插值函数的性态更好。
二维插值函数:
MATLAB中的二维插值函数为interp2(),其调用格式为:
Z1=interp2(X,Y,Z,X1,Y1,method)
其中,X、Y是两个向量,表示两个参数的采样点,Z是采样点对应的函数值。X1、Y1是两个标量或向量,表示要插值的点。

3 应用案例1-粮储仓的通风控制问题

在某粮情自动测控系统中,根据粮温、粮湿计算平衡点湿度,与大气湿度进行比较,再根据通风模拟情况决定是否自动进行通风。已测得平衡点湿度与粮温、粮湿关系的部分数据如下表,请推算相应范围内温度每变化1度、湿度每变化1个点的平衡点湿度。

平衡点湿度与粮温、粮湿度关系表

(第一列:粮温,第一行:粮湿,其余:平衡点湿度)

粮温\粮湿 20 30 40 50 60 70 80 90
0 8.9 10.32 11.3 12.5 13.9 15.3 17.8 21.3
5 8.7 10.8 11 12.1 13.2 14.8 16.55 20.8
10 8.3 9.65 10.88 12 13.2 14.6 16.4 20.5
15 8.1 9.4 10.7 11.9 13.1 14.5 16.2 20.3
20 8.1 9.2 10.8 12 13.2 14.8 16.9 20.9
x=20:10:90;
y=(0:5:20)';
z=[8.9,10.32,11.3,12.5,13.9,15.3,17.8,21.3;8.7,10.8,11,12.1,13.2,14.8,16.55,20.8;8.3,9.65,10.88,12,13.2,14.6,16.4,20.5;8.1,9.4,10.7,11.9,13.1,14.5,16.2,20.3;8.1,9.2,10.8,12,13.2,14.8,16.9,20.9];
xi=20:90;
yi=(0:20)';
zi=interp2(x,y,z,xi,yi,'spline');
surf(xi,yi,zi)

fig3

4 应用案例2-机动车刹车距离问题

在车辆行驶中,从驾驶员看到障碍物开始,到作出判断而采取制动措施停车所需的最短距离叫停车视距。停车视距由三部分组成:一是驾驶员反应时间内行驶的距离(即反应距离);二是开始制动到车辆完全停止所行驶的距离(即制动距离);三是车辆停止时与障碍物应该保持的安全距离。其中,制动距离主要与行驶速度和路面类型有关。根据测试,某型车辆在潮湿天气于沥青路面行驶时,其行车速度(单位:km/h)与制动距离(单位:m)的关系如下表所示。

速度 20 30 40 50 60 70 80 90 100 110 120 130 140 150
制动距离 3.15 7.08 12.59 19.68 28.34 38.57 50.4 63.75 78.71 95.22 113.29 132.93 154.12 176.87

假设驾驶员的反应时间为10s,安全距离为10m。请问:
①根据某驾驶员的实际视力和视觉习惯,其驾驶时的有效视距为120m,则其在该路面行车时,时速最高不能超过多少(结果取整)?
②若以表中数据为参考,设计一条最高时速为125km/h的高速公路,则设计人员应该保证驾驶者在公路上任一点的可视距离为多少米?

设速度为vvv,停车视距为ddd,反应距离为d1d_1d1,制动距离为d2d_2d2,安全距离为d3d_3d3,反应时间为asa_sas,则
d=d1+d2+d3d=d_1+d_2+d_3d=d1+d2+d3其中,d1=asvd_1=a_svd1=asvd2d_2d2为v的函数,d3d_3d3已知。
第一问:根据某驾驶员的实际视力和视觉习惯,其驾驶时的有效视距为120m,则其在该路面行车时,时速最高不能超过多少(结果取整)?
已知反应时间为10s,安全距离为10m,可采用解方程方法:
10v+d2+10=12010v+d_2+10=12010v+d2+10=120存在的问题是,d2d_2d2vvv的函数,但是函数关系未知,方程不可解。
下面考虑数据插值方法,以表格中的数据为样本,进行数据插值,计算出与120m的停车视距所对应的速度指标。
编程思路:
第一步:建立速度和停车视距向量。
第二步:以1为单位,对采样区间内所有速度进行插值,计算出相应的停车视距。
第三步:求出停车视距120所对应的速度。
第四步:绘图展示。
如何根据停车视距120找到对应的速度?
第一步:令代表停车视距的向量di减去120,再取绝对值,得到一个新的向量x。
第二步:将x按升序排列,并记录最小元素的序号,该序号即为停车视距120所对应的速度数据在向量vi中的序号。
第三步:根据序号取得速度数据。

v=20:10:150;
vs=v.*(1000/3600);
d1=10.*vs;
d2=[3.15,7.08,12.59,19.68,28.34,38.57,50.4,63.75,78.71,95.22,113.29,132.93,154.12,176.87];
d3=10;
d=d1+d2+d3;
vi=20:1:150;
di=interp1(v,d,vi,'spline');
x=abs(di-120);
[y,i]=sort(x);
vi(i(1))
plot(vi,di,vi(i(1)),di(i(1)),'rp')

fig4
停车视距的增长随着车速增加呈非线性增长。速度越快,要求视线越远。
第二问:设计一条最高时速为125km/h的高速公路,则设计人员应该保证驾驶者在公路上任一点的可视距离为多少米?

>>j=find(vi==125);
>>di(j)
>>ans=480.14
>>plot(vi,di,125,480.14,'rp')

5 应用案例3-沙盘制作问题

某地面部队分成红蓝两方在指定的陌生区域(平面区域[0,2000]*[0,2000]内,单位:m)进行作战演习。在演习过程中,红方侦查单位已经测得一些地点的高程如下表所示。

y/x 0 200 400 600 800 1000 1200 1400 1600 1800
0 2000 2000 2001 1992 1954 1938 1972 1995 1999 1999
200 2000 2002 2006 1908 1533 1381 1728 1959 1998 2000
400 2000 2005 2043 1921 977 897 1310 1930 2003 2000
600 1997 1978 2009 2463 2374 1445 1931 2209 2050 2003
800 1992 1892 1566 1971 2768 2111 2653 2610 2121 2007
1000 1991 1875 1511 1556 2221 1986 2660 2601 2119 2007
1200 1996 1950 1797 2057 2849 2798 2608 2303 2052 2003
1400 1999 1999 2079 2685 3390 3384 2781 2165 2016 2000
1600 2000 2002 2043 2271 2668 2668 2277 2049 2003 2000
1800 2000 2000 2004 2027 2067 2067 2027 2004 2000 2000

①根据表中数据,制作军事沙盘。
②在演习范围内,占领最大高地的一方将获得居高临下的优势。请问红方应第一时间抢占哪块区域。

解题思路:
第一问:用二维插值估算数据,以方便制作军事沙盘。
第二问:在插值的基础上,绘制等高线图,找到最大高地。

x=0:200:1800;
y=x';
z=[2000,2000,2001,1992,1954,1938,1972,1995,1999,1999;
2000,2002,2006,1908,1533,1381,1728,1959,1998,2000;
2000,2005,2043,1921,977,897,1310,1930,2003,2000;
1997,1978,2009,2463,2374,1445,1931,2209,2050,2003;
1992,1892,1566,1971,2768,2111,2653,2610,2121,2007;
1991,1875,1511,1556,2221,1986,2660,2601,2119,2007;
1996,1950,1797,2057,2849,2798,2608,2303,2052,2003;
1999,1999,2079,2685,3390,3384,2781,2165,2016,2000;
2000,2002,2043,2271,2668,2668,2277,2049,2003,2000;
2000,2000,2004,2027,2067,2067,2027,2004,2000,2000];
surf(x,y,z);
x1=0:100:1800;
y1=x1';
z1=interp2(x,y,z,x1,y1,'spline');
surf(x1,y1,z1);
x2=0:50:1800;
y2=x2';
z2=interp2(x1,y1,z1,x2,y2,'spline');
surf(x2,y2,z2);
contour(x2,y2,z2,12)

fig5
从图中可以看出最大高地的位置。

总结

这篇文章是否对您有帮助呢?最后
欢迎大家👍 收藏⭐ 转发🚀,
如有问题、建议请您在评论区留言💬。

Published by

风君子

独自遨游何稽首 揭天掀地慰生平