1、第十章 插值与拟合方法建模在生产实际中，常常要处理由 实验或测量所得到的一批离散数据，插值与拟合方法就是要通过这些数据去确定某一类已经函数的参数，或 寻求某个近似函数使之与已知数据有 较高的拟合精度。插值与拟合的方法很多，这里主要介绍线性插值方法、多项式插值方法和样条插值方法，以及最小二乘拟合方法在实际问题中的应用。相 应的理论和算法是数 值分析的内容， 这里不作详细介绍， 请参阅有关的书籍。1 数据插值方法及应用在生产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的一批离散样点，要求由此建立变量之间的函数关系或得到样点之外的数据。与此有关的一类问题是当原始数据精度较高，要求确定一
2、个初等函数 （一般用多项式或分),(),(,10nyxyx )(xPy段多项式函数）通过已知各数据点（节点） ，即 ，或要求得函数在另外nixyi ,10,)(一些点（插值点）处的数值，这便是插值问题。1、分段线性插值这是最通俗的一种方法，直观上就是将各数据点用折线连接起来。如果 bxxan10那么分段线性插值公式为 niyxyxPiiiii ,21,)( 1111 可以证明，当分点足够细时，分段线性插值是收敛的。其缺点是不能形成一条光滑曲线。例 1、已知欧洲一个国家的地图，为了算出它的国土面积，对地图作了如下测量：以由西向东方向为 x 轴，由南向北方向为 y 轴，选择方便的原点，并将从最西边
3、界点到最东边界点在 x 轴上的区间适当的分为若干段，在每个分点的 y 方向测出南边界点和北边界点的 y 坐标 y1 和 y2，这样就得到下表的数据（单位：mm） 。x 7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0y1 44 45 47 50 50 38 30 30 34y2 44 59 70 72 93 100 110 110 110x 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5y1 36 34 41 45 46 43 37 33 28y2 117 118 116 118 118 121 124 12
4、1 121x 111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0y1 32 65 55 54 52 50 66 66 68y2 121 122 116 83 81 82 86 85 68根据地图的比例，18 mm 相当于 40 km。根据测量数据，利用 MATLAB 软件对上下边界进行线性多项式插值，分别求出上边界函数，下边界函数 ，利用求平面图形面积的数值积分方法将该面积近似分成若干个小)(2xf )(1xf长方形，分别求出这些长方形的面积后相加即为该面积的近似解。 iiniixffS)(lm112式中， 。,1iiix这里线性插值和面
5、积计算源程序如下：clear allx=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
6、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、多项式插值设有 次多项式mmmaxxaxP110)(通过所有 个点 ，那么就有1n),(,(10nyy niyxxamimii ,0,1 可以证明当 且 时，这样的多项式存在且唯一。若要求得到函数表达式，mn10可直接解上面方
7、程组。若只要求得函数在插值点处数值，可用下列 Lagrange 插值公式)()(,0nijjinxyxP多项式插值光滑但不具有收敛性，一般不宜采用高次多项式（如 ）插值。7m例 2、在万能拉拨机中有一个园柱形凸轮，其底园半径 R=300mm，凸轮的上端面不在同一平面上，而要根据动杆位移变化的需要进行设计制造。按设计要求，将底园周 18 等分，旋转一周。第 个分点对应柱高 ，数据见下表。为了数控加工，需要计算出园周上任一i )18,20(iy点的柱高。凸轮高度的数据（单位： mm）分点 i0 和 18 1 2 3 4 5柱高 502.8 525.0 514.3 451.0 326.5 188.6
8、分点 i6 7 8 9 10 11柱高 92.2 59.6 62.2 102.7 147.1 191.6分点 12 13 14 15 16 17柱高 236.0 280.5 324.9 369.4 413.8 458.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,502.8
9、;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、样条插值这是最常用的插值方法。数学上所说的样条，实质上是指分段多项式的光滑连接。设有 bxxan10称分段函数 为 k 次样条函数,若它满足)(xS(1) 在每个小区间上是次数不超过 次的多项式；k(2) 在 上具有直到 阶的连续导数。)(x,ba1用样条函数作出的插值称为样条插值。工程上广泛采用三次样条插值。例 3、某居民区的自来水
10、是由一个园柱形的水塔提供。水塔高 12.2 米，直径 17.4 米。水塔由水泵根据塔中水位高低自动加水，一般每天水泵工作两次。按照设计，当水塔内的水位降至约 8.2米时，水泵自动启动加水；当水位升至约 10.8 米时，水泵停止工作。现在需要了解该居民区用水规律，这可以通过用水率（单位时间的用水量）来反映。通过间隔一段时间测量水塔中的水位来估算用水率。下表是某一天的测量记录数据，测量了 28 个时刻（单位：小时）的水位（单位：米） ，但由于其中有 3 个时刻正遇到水泵在向水塔供水，而无水位记录（表中用符号/表示） 。时刻 0 0.921 1.843 2.949 3.871 4.978 5.900
11、水位 9.677 9.479 9.308 9.125 8.982 8.814 8.686时刻 7.006 7.982 8.967 9.981 10.925 10.954 12.032水位 8.525 8.388 8.220 / / 10.820 10.500时刻 12.954 13.875 14.982 15.903 16.826 17.931 19.037水位 10.210 9.936 9.653 9.409 9.180 8.921 8.662时刻 19.959 20.839 22.015 22.958 23.880 24.986 25.908水位 8.433 8.220 / 10.820
12、10.591 10.354 10.180先通过体积公式 ，利用上表中的水位高 ，得到不同时刻 水塔中水的体积 。为hdv24hitiv提高精度，采用二阶差商来估算 时刻的水流速度，即 。it iivtf2)(具体地，因为所有数据被水泵两次工作分割成三组数据，对每组数据的中间数据采用中心差商，前后两个数据不能够采用中心差商，改用向前或向后差商。中心差商公式 )(128212iiiii tvvv向前差商公式 )(3412iiii t向后差商公式 )(3122iiitvv估算出水塔中水的流速（单位：立方米/小时）见下表。时刻 0 0.921 1.843 2.949 3.871 4.978 5.900
13、流速 54.516 42.320 38.085 41.679 33.297 37.814 30.748时刻 7.006 7.982 8.967 9.981 10.925 10.954 12.032流速 38.455 32.122 41.718 / / 73.686 76.434时刻 12.954 13.875 14.982 15.903 16.826 17.931 19.037流速 71.686 60.190 68.333 59.217 52.011 56.626 63.023时刻 19.959 20.839 22.015 22.958 23.880 24.986 25.908流速 54.85
14、9 55.439 / 57.602 57.766 51.891 36.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.
15、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(流速( 立方米/ 小时)
