第十章第十章 插值与拟合方法建模插值与拟合方法建模在生产实际中，常常要处理由实验或测量所得到的一批离散数据，插值与拟合方法就是要通 过这些数据去确定某一类已经函数的参数，或寻求某个近似函数使之与已知数据有较高的拟合精 度。插值与拟合的方法很多，这里主要介绍线性插值方法、多项式插值方法和样条插值方法，以及 最小二乘拟合方法在实际问题中的应用。相应的理论和算法是数值分析的内容，这里不作详细介 绍，请参阅有关的书籍。1 数据插值方法及应用在生产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的一批离散样点， 要求由此建立变量之间的函数关系或得到样点之外的数据。与此有关的一类问题是当原始数据精度较高，要求确定一个初等函数（一般用多项式或分),( ,),(),(1100nnyxyxyx)(xPy 段多项式函数）通过已知各数据点（节点） ，即，或要求得函数在另外nixPyii, 1 , 0, )(一些点（插值点）处的数值，这便是插值问题。 1、分段线性插值这是最通俗的一种方法，直观上就是将各数据点用折线连接起来。如果bxxxan10那么分段线性插值公式为nixxxyxxxxyxxxxxPiii iii i iii, 2 , 1,)(1 11 1 1 可以证明，当分点足够细时，分段线性插值是收敛的。其缺点是不能形成一条光滑曲线。 例例 1、已知欧洲一个国家的地图，为了算出它的国土面积，对地图作了如下测量：以由西向东方 向为 x 轴，由南向北方向为 y 轴，选择方便的原点，并将从最西边界点到最东边界点在 x 轴上 的区间适当的分为若干段，在每个分点的 y 方向测出南边界点和北边界点的 y 坐标 y1 和 y2，这 样就得到下表的数据（单位：mm） 。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.0y1326555545250666668y2121122116838182868568根据地图的比例，18 mm 相当于 40 km。根据测量数据，利用 MATLAB 软件对上下边界进行线性多项式插值，分别求出上边界函数，下边界函数，利用求平面图形面积的数值积分方法将该面积近似分成若干个小)(2xf)(1xf长方形，分别求出这些长方形的面积后相加即为该面积的近似解。iiniinxffS )()(lim1 12式中，。,1iiixx这里线性插值和面积计算源程序如下： clear all x=7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0; y1=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; y2=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; newx=7:0.1:158; newy1=interp1(x,y1,newx,linear); newy2=interp1(x,y2,newx,linear); Area=sum(newy2- newy1)*0.1/182*1600 最后计算的面积约为 42414 平方公里。2、多项式插值设有次多项式mmmmmaxaxaxaxP 11 10)(通过所有个点，那么就有1n),( ,),(),(1100nnyxyxyxniyaxaxaxaimimm im i, 1 , 0,11 10可以证明当且时，这样的多项式存在且唯一。若要求得到函数表达式，nm nxxx10可直接解上面方程组。若只要求得函数在插值点处数值，可用下列 Lagrange 插值公式)()(, 00 nijjjijniinxxxxyxP多项式插值光滑但不具有收敛性，一般不宜采用高次多项式（如）插值。7m 例例 2、在万能拉拨机中有一个园柱形凸轮，其底园半径 R=300mm，凸轮的上端面不在同一平面 上，而要根据动杆位移变化的需要进行设计制造。按设计要求，将底园周 18 等分，旋转一周。第 个分点对应柱高，数据见下表。为了数控加工，需要计算出园周上任一i)18, 2 , 1 , 0(iyi点的柱高。 凸轮高度的数据（单位：mm）分点i0 和 1812345柱高502.8525.0514.3451.0326.5188.6分点i67891011柱高92.259.662.2102.7147.1191.6分点i121314151617柱高236.0280.5324.9369.4413.8458.3我们将园周展开，借助 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,5 02.8; plot(x,y,o);axis(0,2000,0,550);可见，可以用三次多项式插值，下面给出借助 MATLAB 软件画出的柱高插值曲线图(右上图)。 xi=0:2*pi*300; yi=interp1(x,y,xi,cubic); plot(xi,yi); 3、样条插值 这是最常用的插值方法。数学上所说的样条，实质上是指分段多项式的光滑连接。设有bxxxan10称分段函数为 k 次样条函数次样条函数,若它满足)(xS(1) 在每个小区间上是次数不超过次的多项式；)(xSk(2) 在上具有直到阶的连续导数。)(xS,ba1k用样条函数作出的插值称为样条插值。工程上广泛采用三次样条插值三次样条插值。 例例 3、某居民区的自来水是由一个园柱形的水塔提供。水塔高 12.2 米，直径 17.4 米。水塔由水 泵根据塔中水位高低自动加水，一般每天水泵工作两次。按照设计，当水塔内的水位降至约 8.2 米时，水泵自动启动加水；当水位升至约 10.8 米时，水泵停止工作。现在需要了解该居民区用 水规律，这可以通过用水率（单位时间的用水量）来反映。通过间隔一段时间测量水塔中的水 位来估算用水率。下表是某一天的测量记录数据，测量了 28 个时刻（单位：小时）的水位（单 位：米） ，但由于其中有 3 个时刻正遇到水泵在向水塔供水，而无水位记录（表中用符号/表示） 。时刻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.180先通过体积公式，利用上表中的水位高，得到不同时刻水塔中水的体积。为hdv2 4hitiv提高精度，采用二阶差商来估算时刻的水流速度，即。itiivtf2)(具体地，因为所有数据被水泵两次工作分割成三组数据，对每组数据的中间数据采用中心 差商，前后两个数据不能够采用中心差商，改用向前或向后差商。中心差商公式 )(1288121122iiiiii ittvvvvv向前差商公式 )(2341122iiiii ittvvvv向后差商公式 )(2431212 iiiii ittvvvv估算出水塔中水的流速（单位：立方米/小时）见下表。时刻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 画出水流速散点图。 t=0 0.921 1.843 2.949 3.871 4.978 5.9 7.006 7.982 8.967 10.954 12.032 12.954 13.875 14.982 15.903 16.826 17.931 19.037 19.959 20.839 22.958 23.88 24.986 25.908; r=54.516 42.320 38.085 41.679 33.297 37.814 30.748 38.455 32.122 41.718 73.686 76.434 71.686 60.19 68.333 59.217 52.011 56.626 63.023 54.859 55.439 57.602 57.766 51.891 36.464; plot(t,r,b+); % (t,r)表示时间和流速 title(流速散点图)；xlabel（时间(小时)）; ylabel(流速(立方米/小时)使用 MATLAB 软件中的三次样条插值命令得到用水率函数如下图所示。)(tfx0=t;y0=r; l,n=size (x0); dl=x0(n)-x0(1); x=x0(1):1/3600:x0(n); %被插值点 ys=interp1 (x0,y0,x,spline); %样条插值输出plot (x,ys); title(样条插值下的流速图)；xlabel（时间(小时)）; ylabel(流速(立方米/小时)