第3章插值与拟合

第3章插值与拟合方法

随着社会的进步和收入水平的提高，汽车进入家庭已不再是奢望。但伴随而来的就是交通安全。“珍爱生命，安全出行”，并不仅仅是个口号，它关系到每个驾驶员的安全，也关系到每个驾驶员所在家庭的幸福和安定。驾驶时，车速过快、与前车距离过近，以致来不及刹车或制动距离不足，是造成绝大部分交通事故的主要原因。

统计上，刹车距离由反应距离和制动距离两部分组成，即刹车距离为反应距离与制动距离之和。前者指从司机发现问题决定刹车到制动器开始起作用汽车行驶的距离，后者指从制动器开始起作用到汽车完全停止行驶的距离.

为了了解刹车距离与车速的关系，美国交通部门进行了一系列刹车实验，实验结果见表3-1所示。

问若车速分别为37、72英里/小时(分别约60、115Km/h)，问刹车距离是多少？保持多大车距才是安全的？

显然，实际观测没有针对这两个点的观测结果，这就需要我们根据已有的观测数据进行估算。进一步地，若要估算车速在区间[20,80]（英里/小时）内任意一点的反应距离、制动距离和刹车距离，应如何估算。

处理此类问题，插值方法与数据拟合方法是两类常见的建模方法

3.1插值法

3.1.1问题的提出

插值问题的一般描述：若已知函数（通常为未知）在给定的个互不相同的观测点上的函数值（通常为实验或观测值），希望寻求某一近似函数,使满足

（3.1）

则我们称此类问题为插值问题，近似函数称为插值函数，观测点称为插值节点，式（3.1）称为插值条件，若令,则[a,b]称为插值区间。

若已找到，则在任一点（）上的函数值就可以由其插值函数近似估计。

那么应该如何构造插值函数呢？从中学的解析几何知识，我们知道：给定平面上两个互不相同的点可以确定一条直线，给定三个互不相同的点可以确定一条抛物线多项式，依此类推。这启示我们用多项式作为插值函数是一个很好的选择。事实上，多项式插值由于其易求导、求积分和足够的光滑性，在很多领域都有广泛的应用。

设是个互不相同的观测点，要求一个次数不超过的代数多项式

（3.2）

使其在插值节点上，满足

（3.3）

则此类插值问题称为代数插值问题，称为次插值多项式。

3.1.3Lagrange多项式插值方法

线性插值：任给两个互不相同的观测点，求一个线性次多项式，使其满足插值条件。线性插值多项式可直接给出，如(3.4)式，但为了引出Lagrange插值多项式的构造思想，我们把它重新组合

合并前两项，整理后得

（3.6）

令

则线性插值多项式可重写为（3.7）

注意到都是线性多项式，二者的线性组合仍然至多是线性多项式。可以验证，由（3.7）定义的线性插值多项式一定满足插值条件，即。且有

称为分别对应于插值节点的Lagrange线性插值基函数。

抛物型插值：给定三个互不相同的观测点，,和，求一个次数不超过2次的多项式，使其满足插值条件：

受Lagrange线性插值构造思想启发，我们类似地构造对应于插值节点的二次插值基函数，使其满足

首先确定,由于是二次多项式，且,则易知是二次多项式的根，因此其表达式一定可写为

的形式，其中为待定系数。又由,代入上式得

于是，可得

类似地，可得

进而，Lagrange二次插值（抛物型）多项式可表述为

（3.8）

且也可以很容易地验证上式满足所要求的插值条件。

利用构造插值基函数的思想，可非常方便地给出次Lagrange插值多项式的表达式，有兴趣的同学不妨试一下。理论上，只要给出足够多的观测点，就可以构造任意次插值多项式，但高次插值多项式存在着不可控制的数值震荡现象，在实际问题建模中一般不推荐使用。

分段低次多项式插值方法：

在实际问题观测中，一般会得到很多个观测点的观测结果，采用插值方法近似时，一般采取分段插值的方法。基本思想是：

（1）把插值区间划分成若干个小区间；

（2）在每一小区间上用低次多项式进行插值；

（3）在整个插值区间上就得到一个分段插值函数.

假定给出个互不相同的观测点,不妨设

分段线性插值：把相邻两个插值节点作为一个插值子区间,则插值区间被划分为个子区间，连接相邻两点得条线段,这些线段组成一条折线,这条折线就是我们构造的分段线性插值函数，记为，它具有如下特点。

（1）在整个插值区间上，连续,但在插值节点上不可导；

（2）在第个子区间上，的表达式为

分段二次插值：若,把相邻三个插值节点组成一个插值子区间，则整个插值区间被划分为个子区间。在第个子区间上，共有三个插值节点，为一二次插值多项式，表达式为

例1已知某函数的函数表如下：

用线性插值法估算的近似值.

解：由于在插值节点之间，故依此二点构造Lagrange线性插值多项式，并代入得

即的近似估计值为2.4414.

例2已知观测数据如表3-3所示

试用二次插值方法求处的插值.

解：取包含的三个观测点作为插值节点，作二次插值，并令,由(3.8)式，可得

=1.8903

3.2曲线拟合

3.2.1问题提出

利用插值方法求多项式函数作为未知函数的近似时，要求

1、所有插值节点互不相同，否则不可解；

2、近似函数曲线必须通过所有观测点。

在实际观测或实验中，一般存在以下问题

1、为了得到更加准确、合理的观测结果，经常进行多次重复观测，插值节点互不相同的要求已不成立；

2、由于在观测过程中，常存在许多随机因素，如身高、体重的测量，受测量设备精度、发型、服装、站立方式等影响，测量结果不可避免地存在误差，甚至由于某些因素，误差很大。因此在考虑观测误差的因素下，要求近似函数曲线一定通过观测点已显得没有必要。

因此，只要要求近似函数在观测点上近似地满足插值条件，并使它们的整体误差最小就可以了。

3.2.2基本概念

给定函数(未知)在观测点上的观测值，寻求一近似函数(拟合曲线)，使在所有观测点上，观测值与近似函数的计算值之间的误差

总体上尽可能接近零，即要求尽可能反映给定数据点的总体趋势，这就是函数逼近法，也称为曲线拟合法，

称为逼近函数或拟合函数，曲线称为拟合曲线.

拟合函数的选择范围很广，如多项式，有理函数、指数函数、对数函数、三角多项式等，但具体选择何种函数，应综合多方面因素斟酌确定。比较简单和直观的方法是通过绘制观测点的散点图，进行观察、比较、猜测，然后根据观测结果和误差分析加以确定。

3.2.3最小二乘拟合方法

判断拟合曲线尽可能逼近给定数据点

的标准有很多，如使最大偏差达到极小，所有偏差的绝对值之和取极小等，但因求解方法上的复杂性，实际使用起来并不方便，实践中常用的一种曲线拟合方法就是最小二乘拟合方法.

对给定的数据点，选取拟合函数，使偏差，，的平方和为最小，即：

（3.9）

从几何意义上讲，就是求拟合曲线，使在给定的点

处，计算值与实际观测值的差的平方和最小。这种求近似函数的方法称为离散数据曲线拟合的最小二乘法，函数称为这组数据的最小二乘拟合函数.拟合程度的好坏，可以通过直接计算误差平方和的大小来反映。

若拟合函数取为多项式，如取次多项式

（3.10）

则相应的最小二乘拟合问题就变为：求参数,使

（3.11）

达到极小。由于拟合函数关于待定系数是线性的，故称该问题为线性最小二乘问题，又由于拟合函数是多项式故也称多项式拟合问题。

最简单的多项式拟合是和两种情形，分别称为一次多项式拟合（直线拟合）和二次多项式拟合（抛物线拟合）。以下就这两种简单情形，进行详细讨论。

3.2.4拟合直线

1．数学描述

给定一组散点图呈线性变化的数据设拟合函数，根据最小二乘准则（3.11）知，要求待定系数,使

达到极小。由多元函数极值定理，一个必要条件是两个偏导数等于零，整理得

（3.12）

把给定数据全部带入，用消去法解方程组得出

这就是线性最小二乘拟合的数学表示，方程（3.12）称为正规方程.

引入矩阵记号：

则正规方程（3.12）可以表示为

（3.13）

正规方程的矩阵表示可以推广到次多项式的最小二乘拟合中，给定的数据，此时的正规方程用矩阵表示为式（3.13），对应的矩阵为：

利用正规方程求拟合多项式的系数，可分为如下几个步骤

1、输入观测向量

2、生成超定矩阵G

3、求

4、计算

5、求解参数向量A。

例3观察下表数据，根据其分布趋势，进行曲线拟合，使拟合曲线和下列数据点间的偏差平方和达到最小.

解：首先作出数据的散点图，如图3-2所示

根据上图，数据点近似直线，所以取作为拟合函数

根据正规方程求解步骤，分别计算如下

(1)输入观测向量：

（2）生成超定矩阵G

（3）求矩阵的转置

（4）计算系数矩阵与右端向量

（5）把正规方程重新写成线性方程组的形式

求解这个方程组，得a=1.1,b=1.0.即.

3.2.5拟合幂曲线

1．数学描述

给定一组散点图呈幂曲线变化的数据，设拟合函数，为固定数，根据最小二乘准则（3.11）知，要极小化

最优化的必要条件为导数等于零，即

整理得

于是，最小二乘拟合幂曲线的正规方程用矩阵表示为

（3.14）

其中

，,解得，其中是幂指数.

例4观察下表数据，根据其分布趋势，进行曲线拟合，使拟合曲线和下列数据点间的偏差平方和达到最小，并预测x=2.25时y的值.

解：首先利用(x,y)作散点图，如图3-3(a)所示，发现大概呈二次抛物线形式，为验证结果，再利用(x2,y)作图,如图3-3(b)所示，二者呈良好的直线关系。故可以二次幂函数作为拟合函数。取拟合函数形式为，根据（3.14）式，给出求解过程如下：

(a)(b)

（1）输入并形成观测向量

(2)生成超定矩阵G

(3)计算

(4)求参数a:

将x=2.25，代入得y=16.1193

3.3利用Matlab求解插值与拟合问题

3.3.1利用Matlab求解一维插值问题

Matlab用于求解一维插值问题的基本函数为interp1，给出了四种插值方法供选择，分别为:最邻近点插值（零次多项式插值）、线性插值、三次插值和三次样条插值，其调用的基本格式为:

Interp1(x,y,Cx，’Method’)

其中x,y分别表示为数据点的横、纵坐标向量,x必须单调.

Cx为需要插值的横坐标数据(或数组),Cx不能超出x的范围.

Method为可选参数,对应于上述四种方法,可从以下四个值中任选一个:

‘nearest’---最近邻点插值

‘linear’---线性插值

‘spline’-----三次样条插值

‘cubic’------三次插值

‘linear’是缺省值.当忽略选项时默认为线性插值.

3.3.2利用Matlab求解曲线拟合问题

1利用正规方程求解，可先把正规方程化为标准矩阵代数方程的形式，然后利用Matlab命令

x=A\b

2利用多项式拟合函数polyfit()求解：

调用格式：A=polyfit(x,y,m)

格式说明：x,y分别对应观测点和观测值，是两个同维向量，m是拟合多项式的次数；返回结果为所拟合的m次多项式的系数向量，排列次序为由幂次从高到低。

若要利用所求的拟合多项式估计某些点的函数值，可利用多项式求值函数polyval()

调用格式：f=polyval(A,cx)

格式说明：A为多项式函数的系数向量，按幂次由高到低排列，cx为要估计的点，可以是一个值，也可以是一个向量。

3.3.3应用示例

例5温度预测问题

在12小时内，每隔1小时测量一次温度，详细数据如下：

（1）试分别用线性插值、三次样条插值估计在3.2、6.5、7.1、11.7h的温度值，用三次样条插值法每隔0.1h估计一次温度值并画出其图形.

（2）用多项式拟合，估计在3.2、6.5、7.1、11.7h的温度值.

问题分析：

已知一组数据，要估计插值点的值，既可以用插值法，也可以用数据拟合法，拟合出函数关系式，再求给定点的函数值（拟合值）.

模型求解：

插值方法：分别利用线性插值和三次样条方法估计,编辑Matlab程序如下

t=[123456789101112];

T=[589152529313022252724];

t0=[3.26.57.111.7];

T0=interp1(t,T,t0,‘linear’)

T1=interp1(t,T,t0,‘spline’)

t1=1:0.1:12;

T2=interp1(t,T,t1,‘spline’);

plot(t,T,‘x’);

holdon

plot(t1,T2,‘-’);

legend(‘观测值’,‘拟合值’);

xlabel(‘时间’);

ylabel(‘温度’);

运行结果为：

T0=10.200030.000030.900024.9000

T1=9.673430.042731.175525.3820

每隔0.1h估计一次温度值并作图结果见图3-4

（2）用三次多项式拟合，编程如下：

t=[123456789101112];

T=[589152529313022252724];

plot(t,T,'*');

holdon

p=polyfit(t,T,3);%求三次拟合多项式的系数

t0=[3.26.57.111.7];%插值点

T0=polyval(p,t0)%利用所求三次多项式，估计在插值点的值

plot(t0,T0,'x');%绘制观测点散点图

holdon

t1=1:0.1:12;

T2=polyval(p,t1);

plot(t1,T2,‘-’)%T0为拟合值

title(‘三次多项式拟合’)

legend(‘观测值’,‘拟合值’,‘拟合曲线’)

xlabel(‘时间’)

ylabel(‘温度’)

运行结果得，利用三次多项式拟合在3.2、6.5、7.1、11.7h的估计值为：

T0=

14.801726.250027.308823.6551

输出图形如图3-5所示。

例6土豆产量与化肥的关系

在农业生产试验研究中，对某地区土豆的产量与化肥的关系做了一实验，得到了氮肥、磷肥的施肥量与土豆产量的对应关系如下表：

根据上表数据分别给出土豆产量与氮、磷肥的关系式.

问题分析：

首先画出土豆产量与氮施肥量的散点图，见图3-6.

从图3-6可以看出，土豆产量与氮肥量的关系是二次函数关系，因此可选取拟合函数为：

其中x和y分别为氮肥施肥量和土豆产量，a,b,c为待定系数.

再画出磷肥量与土豆产量的散点图，见图3-7.

从图3-7可以看出：当磷肥施肥量每公顷达到达到100公斤左右时，两侧曲线分别呈明显的线性关系.由此可选取分段的线性函数作为近似函数，即用的观测点作一线性拟合函数，再用的点作一线性拟合函数，最后用两个线性函数求出其分界点即可得分段线性函数.

模型求解：

（1）调用多项式拟合函数p=polyfit(x,y,n)及求值函数y=polyval(p,x)，编程如下：

x=[03467101135202259336404471];

y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];

p=polyfit(x,y,2)

y0=polyval(p,x);

plot(x,y,'x',x,y0,'-k')

legend('观测值','拟合值')

xlabel('氮肥施肥量(Kg/ha)')

ylabel('土豆产量(Kg)')

title('土豆产量与氮肥施肥量的二次多项式拟合')

运行该程序得到拟合二次多项式的系数为

p=

-0.00030.197114.7416

即所求拟合函数为

土豆产量与氮肥量的拟合数据与实验数据的比较，如图3-8.

（2）调用线性拟合函数A=regress(Y,G)，编程如下：

x1=[024497398];

y1=[33.4632.4736.0637.9641.04];

L1=polyfit(x1,y1,1)%求前段线性拟合多项式

x2=[98147196245294342];

y2=[41.0440.0941.2642.1740.3642.73];

L2=polyfit(x2,y2,1)%求后段线性拟合多项式

z1=polyval(L1,x1);%估计前段拟合值

z2=polyval(L2,x2);%估计后段拟合值

plot(x1,y1,'*',x1,z1,'--b')

holdon

plot(x2,y2,'+',x2,z2,'--b')

xlabel('磷肥施肥量(Kg/ha)')

ylabel('土豆产量(Kg)')

holdoff

运行得：

L1=

0.084432.0771

L2=

0.005939.9685

即左侧拟合直线方程为：，

即左侧拟合直线方程为：，

联立求解这两个方程，得两条直线的交点，100.53，则拟合函数为

为便于观察，拟合多项式对观测数据的拟合情况，程序同时利用拟合结果作了对比分析，见图3-9所示。

3.4汽车刹车距离模型

问题分析：反应距离由反应时间和车速决定，反应时间取决于司机个人状况(灵巧、机警、视野等)和制动系统的灵敏性(从司机脚踏刹车板到制动器真正起作用的时间)，对于一般规则可以视反应时间为常数，且在这段时间内车速尚未改变.制动距离与制动器作用力(制动力)、车重、车速以及道路、气候等因素有关，制动器是一个能量耗散装置，制动力作的功被汽车动能的改变所抵消.设计制动器的一个合理原则是，最大制动力大体上与车的质量成正比，使汽车的减速度基本上是常数，这样，司机和乘客少受剧烈的冲击.至于道路、气候等因素，对于一般规则又可以看作是固定的.

根据表3-1，取制动距离和刹车距离的平均值，分析一般情况下，反应距离、制动距离和刹车距离各自依赖于速率的情况，作图分析如图3-10—3-12所示

图3-10反应距离与速率图3-11制动距离与速率

图3-12刹车距离与速率

由图3-10可以看出：反应距离跟速率成比例，图3-11-图3-12表明制动距离跟速率符合二次多项式，为了找到刹车距离随速率的变化函数，基于以上分析，作如下假设.

基本假设

（1）刹车距离等于反应距离和制动距离之和；

（2）反应距离与车速成正比，比例系数为反应时间；

（3）刹车时使用最大制动力，做的功等于汽车动能的改变，且与车的质量成正比.

模型建立

由假设（2）可得

由假设（3），在作用力下行驶距离做的功，所做的功使车速从变成0，动能的变化为，有

又，按照牛顿第二定理可知，刹车时的减速度为常数，于是

其中，由假设（1），刹车距离为

即实际刹车距离的拟合多项式为

（通常对人们刹车的反应时间经验估计为0.75秒，即拟合函数中.）

模型求解

根据最小二乘拟合准则（3.11），要拟合参数，就要极小化

最优的一个必要条件是两个偏导数等于零，从而得到矩阵表示的正规方程为

其中

为了便于分析，我们对表3-1的数据做了适度整理，见表3-8所示。

根据表3-8的数据，分别由制动距离、反应距离拟合参数，根据拟合幂曲线中式（3.14）拟合，编程如下：

v=[29.336.74451.358.66673.380.78895.3102.7110117.3];

d1=[202840.552.57292.5118148.5182220.5266318376];

g=[v‘.*v’];

k1=g\d1‘;

得到k1=0.0252

d2=[22283339445055616672778388];

g=[v’];

k2=g\d2‘;

得到k2=0.7528

所以，作图：

K=[k1k2];

d=[v’.*v’v’]*K’;

plot(v,d,’