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