数据插值方法

2020/1/30差值方法第十章插值与拟合方法建模2020/1/30差值方法在生产实际中，常常要处理由实验或测量所得到的一批离散数据，插值与拟合方法就是要通过这些数据去确定某一类已经函数的参数，或寻求某个近似函数使之与已知数据有较高的拟合精度。插值与拟合的方法很多，这里主要介绍线性插值方法、多项式插值方法和样条插值方法，以及最小二乘拟合方法在实际问题中的应用。相应的理论和算法是数值分析的内容，这里不作详细介绍，请参阅有关的书籍。2020/1/30差值方法§1数据插值方法及应用在生产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的一批离散样点，要求由此建立变量之间的函数关系或得到样点之外的数据。与此有关的一类问题是当原始数据),(,),,(),,(1100nnyxyxyx精度较高，要求确定一个初等函数)(xPy（一般用多项式或分段多项式函数）通过已知各数据点（节点），即nixPyii,,1,0,)(，或要求得函数在另外一些点（插值点）处的数值，这便是插值问题。2020/1/30差值方法1、分段线性插值这是最通俗的一种方法，直观上就是将各数据点用折线连接起来。如果bxxxan10那么分段线性插值公式为nixxxyxxxxyxxxxxPiiiiiiiiii,,2,1,,)(11111可以证明，当分点足够细时，分段线性插值是收敛的。其缺点是不能形成一条光滑曲线。2020/1/30差值方法例1、已知欧洲一个国家的地图，为了算出它的国土面积，对地图作了如下测量：以由西向东方向为x轴，由南向北方向为y轴，选择方便的原点，并将从最西边界点到最东边界点在x轴上的区间适当的分为若干段，在每个分点的y方向测出南边界点和北边界点的y坐标y1和y2，这样就得到下表的数据（单位：mm）。2020/1/30差值方法x7.010.513.017.534.040.544.548.056.0y1444547505038303034y24459707293100110110110x61.068.576.580.591.096.0101.0104.0106.5y1363441454643373328y2117118116118118121124121121x111.5118.0123.5136.5142.0146.0150.0157.0158.0y1326555545250666668y21211221168381828685682020/1/30差值方法根据地图的比例，18mm相当于40km。根据测量数据，利用MATLAB软件对上下边界进行线性多项式插值，分别求出上边界函数)(2xf，下边界函数)(1xf，利用求平面图形面积的数值积分方法—将该面积近似分成若干个小长方形，分别求出这些长方形的面积后相加即为该面积的近似解。iiniinxffS)]()([lim112式中，],[1iiixx。这里线性插值和面积计算源程序如下：2020/1/30差值方法clearallx=[7.010.513.017.534.040.544.548.056.061.068.576.580.591.096.0101.0104.0106.5111.5118.0123.5136.5142.0146.0150.0157.0158.0];y1=[444547505038303034363441454643373328326555545250666668];y2=[4459707293100110110110117118116118118121124121121121122116838182868568];newx=7:0.1:158;newy1=interp1(x,y1,newx,’linear’);newy2=interp1(x,y2,newx,’linear’);Area=sum((newy2-newy1)*0.1/18^2*1600)最后计算的面积约为42414平方公里。2020/1/30差值方法2、多项式插值设有m次多项式mmmmaxaxaxaxP1110)(通过所有1n个点),(,),,(),,(1100nnyxyxyx，那么就有niyaxaxaxaimimmimi,,1,0,1110可以证明当nm且nxxx10时，这样的多项式存在且唯一。若要求得到函数表达式，可直接解上面方程组。2020/1/30差值方法若只要求得函数在插值点处数值，可用下列Lagrange插值公式)()(,00nijjjijniinxxxxyxP多项式插值光滑但不具有收敛性，一般不宜采用高次多项式（如m7）插值。2020/1/30差值方法例2、在万能拉拨机中有一个园柱形凸轮，其底园半径R=300mm，凸轮的上端面不在同一平面上，而要根据动杆位移变化的需要进行设计制造。按设计要求，将底园周18等分，旋转一周。第i个分点对应柱高)18,,2,1,0(iyi，数据见下表。为了数控加工，需要计算出园周上任一点的柱高。2020/1/30差值方法凸轮高度的数据（单位：mm）iiy0和18502.81525.02514.33451.04326.55188.6iiy692.2759.6862.29102.710147.111191.6iiy12236.013280.514324.915369.416413.817458.32020/1/30差值方法我们将园周展开，借助MATLAB软件画出对应的柱高曲线散点图(下图)。clear;close;x=linspace(0,2*pi*300,19);y=[502.8,525.0,514.3,451.0,326.5,188.6,92.2,59.6,62.2,102.7,147.1,191.6,236.0,280.5,324.9,369.4,413.8,458.3,502.8];plot(x,y,’o’);axis([0,2000,0,550]);2020/1/30差值方法2020/1/30差值方法可见，可以用三次多项式插值，下面给出借助MATLAB软件画出的柱高插值曲线图（下图）。xi=0:2*pi*300;yi=interp1(x,y,xi,’cubic’);plot(xi,yi);2020/1/30差值方法2020/1/30差值方法3、样条插值这是最常用的插值方法。数学上所说的样条，实质上是指分段多项式的光滑连接。设有bxxxan10称分段函数)(xS为k次样条函数,若它满足(1))(xS在每个小区间上是次数不超过k次的多项式;(2))(xS在],[ba上具有直到1k阶的连续导数。用样条函数作出的插值称为样条插值。工程上广泛采用三次样条插值。2020/1/30差值方法例3、某居民区的自来水是由一个园柱形的水塔提供。水塔高12.2米，直径17.4米。水塔由水泵根据塔中水位高低自动加水，一般每天水泵工作两次。按照设计，当水塔内的水位降至约8.2米时，水泵自动启动加水；当水位升至约10.8米时，水泵停止工作。现在需要了解该居民区用水规律，这可以通过用水率（单位时间的用水量）来反映。通过间隔一段时间测量水塔中的水位来估算用水率。2020/1/30差值方法下表是某一天的测量记录数据，测量了28个时刻（单位：小时）的水位（单位：米），但由于其中有3个时刻正遇到水泵在向水塔供水，而无水位记录（表中用符号//表示）。2020/1/30差值方法时刻00.9211.8432.9493.8714.9785.900水位9.6779.4799.3089.1258.9828.8148.686时刻7.0067.9828.9679.98110.92510.95412.032水位8.5258.3888.220////10.82010.500时刻12.95413.87514.98215.90316.82617.93119.037水位10.2109.9369.6539.4099.1808.9218.662时刻19.95920.83922.01522.95823.88024.98625.908水位8.4338.220//10.82010.59110.35410.1802020/1/30差值方法先通过体积公式hdv24，利用上表中的水位高h，得到不同时刻it水塔中水的体积iv。为提高精度，采用二阶差商来估算it时刻的水流速度，即iivtf2)(。具体地，因为所有数据被水泵两次工作分割成三组数据，对每组数据的中间数据采用中心差商，前后两个数据不能够采用中心差商，改用向前或向后差商。2020/1/30差值方法中心差商公式)(1288121122iiiiiiittvvvvv向前差商公式)(2341122iiiiiittvvvv向后差商公式)(2431212iiiiiittvvvv估算出水塔中水的流速（单位：立方米/小时）见下表。2020/1/30差值方法时刻00.9211.8432.9493.8714.9785.900流速54.51642.32038.08541.67933.29737.81430.748时刻7.0067.9828.9679.98110.92510.95412.032流速38.45532.12241.718////73.68676.434时刻12.95413.87514.98215.90316.82617.93119.037流速71.68660.19068.33359.21752.01156.62663.023时刻19.95920.83922.01522.95823.88024.98625.908流速54.85955.439//57.60257.76651.89136.464先用MATLAB画出水流速散点图。2020/1/30差值方法t=[00.9211.8432.9493.8714.9785.97.0067.9828.96710.95412.03212.95413.87514.98215.90316.82617.93119.03719.95920.83922.95823.8824.98625.908];r=[54.51642.32038.08541.67933.29737.81430.74838.45532.12241.71873.68676.43471.68660.1968.33359.21752.01156.62663.02354.85955.43957.60257.76651.89136.464];2020/1/30差值方法plot(t,r,’b+’);%(t,r)表示时间和流速title(‘流速散点图’);xlabel(’时间(小时)’);ylabel(‘流速(立方米/小时)’)2020/1/30差值方法2020/1/30差值方法使用MATLAB软件中的三次样条插值命令得到用水率函数)(tf如下图所示。x0=t;y0=r;[l,n]=size(x0);x=x0(1):1/3600:x0(n);%被插值点ys=interp1(x0,y0,x,’spline’);%样条插值输出plot(x,ys);title(‘样条插值下的流速图’)；xlabel（’时间(小时)’）;ylabel(‘流速(立方米/小时)’)2020/1/30差值方法