第十章 插值与拟合方法建模在生产实际中，常常要处理由实验或测量所得到的一批离散数据，插值与拟合方法就是要通过这些数据去确定某一类已经函数的参数，或寻求某个近似函数使之与已知数据有较高的拟合精度。

插值与拟合的方法很多，这里主要介绍线性插值方法、多项式插值方法和样条插值方法，以及最小二乘拟合方法在实际问题中的应用。

相应的理论和算法是数值分析的内容，这里不作详细介绍，请参阅有关的书籍。

§1 数据插值方法及应用在生产实践和科学研究中，常常有这样的问题：由实验或测量得到变量间的一批离散样点，要求由此建立变量之间的函数关系或得到样点之外的数据。

与此有关的一类问题是当原始数据),(,),,(),,(1100n n y x y x y x 精度较高，要求确定一个初等函数)(x P y =（一般用多项式或分段多项式函数）通过已知各数据点（节点），即n i x P y i i ,,1,0,)( ==，或要求得函数在另外一些点（插值点）处的数值，这便是插值问题。

1、分段线性插值这是最通俗的一种方法，直观上就是将各数据点用折线连接起来。

如果b x x x a n =<<<= 10那么分段线性插值公式为n i x x x y x x x x y x x x x x P i i i i i i i i i i ,,2,1,,)(11111 =≤<--+--=-----可以证明，当分点足够细时，分段线性插值是收敛的。

其缺点是不能形成一条光滑曲线。

例1、已知欧洲一个国家的地图，为了算出它的国土面积，对地图作了如下测量：以由西向东方向为x 轴，由南向北方向为y 轴，选择方便的原点，并将从最西边界点到最东边界点在x 轴上的区间适当的分为若干段，在每个分点的y 方向测出南边界点和北边界点的y 坐标y1和y2，这样就得到下表的数据（单位：mm ）。

根据地图的比例，18 mm 相当于40 km 。

根据测量数据，利用MA TLAB 软件对上下边界进行线性多项式插值，分别求出上边界函数)(2x f ，下边界函数)(1x f ，利用求平面图形面积的数值积分方法—将该面积近似分成若干个小长方形，分别求出这些长方形的面积后相加即为该面积的近似解。

i i ni i n x f f S ∆-=∑=∞→)]()([lim 112ξξ式中，],[1i i i x x -∈ξ。

这里线性插值和面积计算源程序如下： 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 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 次多项式m m m m a x a x a x a x P ++++=--1110)(通过所有1+n 个点),(,),,(),,(1100n n y x y x y x ，那么就有n i y a x a x a x a i m i m m im i ,,1,0,1110 ==++++--可以证明当n m =且n x x x <<< 10时，这样的多项式存在且唯一。

若要求得到函数表达式，可直接解上面方程组。

若只要求得函数在插值点处数值，可用下列Lagrange 插值公式)()(,00∏∑≠==--=nij j ji j ni i n x x x x y x P多项式插值光滑但不具有收敛性，一般不宜采用高次多项式（如7>m ）插值。

例2、在万能拉拨机中有一个园柱形凸轮，其底园半径R=300mm ，凸轮的上端面不在同一平面上，而要根据动杆位移变化的需要进行设计制造。

按设计要求，将底园周18等分，旋转一周。

第i 个分点对应柱高)18,,2,1,0( =i y i ，数据见下表。

为了数控加工，需要计算出园周上任一点的柱高。

凸轮高度的数据（单位：mm ）分点i 0和18 1 2 3 4 5 柱高 502.8 525.0 514.3 451.0 326.5 188.6 分点i 6 7 8 9 10 11 柱高 92.2 59.6 62.2 102.7 147.1 191.6 分点i 12 13 14 15 16 17 柱高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,502.8];plot(x,y,’o ’);axis([0,2000,0,550]);可见，可以用三次多项式插值，下面给出借助MA TLAB 软件画出的柱高插值曲线图(右上图)。

xi=0:2*pi*300;yi=interp1(x,y,xi,’cubic ’); plot(xi,yi);3、样条插值这是最常用的插值方法。

数学上所说的样条，实质上是指分段多项式的光滑连接。

设有b x x x a n =<<<= 10称分段函数)(x S 为k 次样条函数,若它满足(1) )(x S 在每个小区间上是次数不超过k 次的多项式； (2) )(x S 在],[b a 上具有直到1-k 阶的连续导数。

用样条函数作出的插值称为样条插值。

工程上广泛采用三次样条插值。

例3、某居民区的自来水是由一个园柱形的水塔提供。

水塔高12.2米，直径17.4米。

水塔由水泵根据塔中水位高低自动加水，一般每天水泵工作两次。

按照设计，当水塔内的水位降至约8.2米时，水泵自动启动加水；当水位升至约10.8米时，水泵停止工作。

现在需要了解该居民区用水规律，这可以通过用水率（单位时间的用水量）来反映。

通过间隔一段时间测量水塔中的水位来估算用水率。

下表是某一天的测量记录数据，测量了28个时刻（单位：小时）的水位（单位：米），但由于其中有3个时刻正遇到水泵在向水塔供水，而无水位记录（表中用符号//表示）。

先通过体积公式h d v 24=，利用上表中的水位高h ，得到不同时刻i t 水塔中水的体积i v 。

为提高精度，采用二阶差商来估算i t 时刻的水流速度，即i i v t f 2)(-∇=。

具体地，因为所有数据被水泵两次工作分割成三组数据，对每组数据的中间数据采用中心差商，前后两个数据不能够采用中心差商，改用向前或向后差商。

中心差商公式 )(1288121122i i i i i i i t t v v v v v -+-+-=∇+--++向前差商公式 )(2341122i i ii i i t t v v v v --+-=∇+++向后差商公式 )(2431212----+-=∇i i i i i i t t v v v vt=[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 软件中的三次样条插值命令得到用水率函数)(t f 如下图所示。

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(‘流速(立方米/小时)’)。