插值与拟合课件

1、 第3章 插值与拟合方法 随着社会的进步和收入水平的提高，汽车进入家庭已不再是奢望。但伴随而来的就是交通安全。“珍爱生命，安全出行”，并不仅仅是个口号，它关系到每个驾驶员的安全，也关系到每个驾驶员所在家庭的幸福和安定。驾驶时，车速过快、与前车距离过近，以致来不及刹车或制动距离不足，是造成绝大部分交通事故的主要原因。 统计上，刹车距离由反应距离和制动距离两部分组成，即 刹车距离为反应距离与制动距离之和。前者指从司机发现问题决定刹车到制动器开始起作用汽车行驶的距离，后者指从制动器开始起作用到汽车完全停止行驶的距离.第1页，共61页。为了了解刹车距离与车速的关系，美国交通部门进行了一系列刹车实验，实

2、验结果见表3-1所示。 问若车速分别为37、72英里/小时(分别约 60、115 Km/h)，问刹车距离是多少？保持多大车距才是安全的？第2页，共61页。 显然，实际观测没有针对这两个点的观测结果，这就需要我们根据已有的观测数据进行估算。进一步地，若要估算车速在区间20,80（英里/小时）内任意一点的反应距离、制动距离和刹车距离，应如何估算。 处理此类问题，插值方法与数据拟合方法是两类常见的建模方法第3页，共61页。3.1 插值法3.1.1 问题的提出插值问题的一般描述：若已知函数 （通常为未知）在给定的 个互不相同的观测点 上的函数值（通常为实验或观测值） ，希望寻求某一近似函数 , 使满足

3、 （3.1）则我们称此类问题为插值问题，近似函数 称为插值函数，观测点称为插值节点，式（3.1）称为插值条件，若令 , 则a,b称为插值区间。若 已找到，则在任一点 （ ） 上的函数值 就可以由其插值函数 近似估计。 第4页，共61页。那么应该如何构造插值函数呢？从中学的解析几何知识，我们知道：给定平面上两个互不相同的点可以确定一条直线，给定三个互不相同的点可以确定一条抛物线多项式，依此类推。这启示我们用多项式作为插值函数是一个很好的选择。事实上，多项式插值由于其易求导、求积分和足够的光滑性，在很多领域都有广泛的应用。设 是 个互不相同的观测点，要求一个次数不超过 的代数多项式 （3.2）使其

4、在插值节点上，满足 （3.3）则此类插值问题称为代数插值问题，称 为 次插值多项式。第5页，共61页。3.1.2 插值多项式的求法1 一般方法线性插值：给定两个互不相同的观测点 和 , 求一线性多项式使其通过这两个观测点，即 。显然 是平面上的一条直线，其表达式可采用两点式或点斜式直接给出，即 （3.4）当然，也可以利用代数方程组的方法求出待定参数 . 由插值条件， 通过这两个观测点，故有 解此线性方程组，可采用消元法，也可以采用矩阵方法直接求解. 详见3.1.3.第6页，共61页。二次插值：给定三个互不相同的观测点， , 和 ，求一个次数不超过2次的多项式 使其通过这三个观测点。求解方法与线

5、性插值完全类似，此处不再累述。二次插值又称抛物型插值。 次插值多项式：当 大于或等于2时，采用上述方法无法直接给出多项式的表达式，需要求解线性方程组。对 次插值多项式的确定，由于多项式中含有 +1 个待定系数，通常需要给定 +1个互不相同的观测点，由此可建立 +1元线性方程组，如下式： （3.5）直接解此线性方程组，通常比较麻烦，可通过数学软件(如Matlab)求解。第7页，共61页。3.1.3 Lagrange多项式插值方法线性插值：任给两个互不相同的观测点，求一个线性次多项式，使其满足插值条件。线性插值多项式可直接给出，如(3.4)式，但为了引出Lagrange插值多项式的构造思想，我们把

6、它重新组合 合并前两项，整理后得 （3.6）令则线性插值多项式可重写为 （3.7）注意到 都是线性多项式，二者的线性组合仍然至多是线性多项式。可以验证，由（3.7）定义的线性插值多项式一定满足插值条件，即 。且有 称 为分别对应于插值节点 的Lagrange线性插值基函数。第8页，共61页。抛物型插值：给定三个互不相同的观测点， , 和 ，求一个次数不超过2次的多项式 ，使其满足插值条件：受Lagrange 线性插值构造思想启发，我们类似地构造对应于插值节点 的二次插值基函数 ，使其满足首先确定 , 由于是二次多项式，且 ,则易知 是二次多项式 的根，因此其表达式一定可写为的形式，其中 为待定

7、系数。又由 , 代入上式得于是，可得第9页，共61页。类似地，可得进而，Lagrange二次插值（抛物型）多项式可表述为 （3.8）且也可以很容易地验证上式满足所要求的插值条件。 利用构造插值基函数的思想，可非常方便地给出 次Lagrange插值多项式的表达式，有兴趣的同学不妨试一下。理论上，只要给出足够多的观测点，就可以构造任意次插值多项式，但高次插值多项式存在着不可控制的数值震荡现象，在实际问题建模中一般不推荐使用。第10页，共61页。分段低次多项式插值方法：在实际问题观测中，一般会得到很多个观测点的观测结果，采用插值方法近似时，一般采取分段插值的方法。基本思想是：（1）把插值区间划分成若

8、干个小区间；（2）在每一小区间上用低次多项式进行插值；（3）在整个插值区间上就得到一个分段插值函数.第11页，共61页。假定给出 个互不相同的观测点 , 不妨设分段线性插值：把相邻两个插值节点作为一个插值子区间 , 则插值区间被划分为 个子区间，连接相邻两点 得 条线段, 这些线段组成一条折线, 这条折线就是我们构造的分段线性插值函数，记为 ，它具有如下特点。（1）在整个插值区间上， 连续 ,但在插值节点上不可导；（2）在第 个子区间 上， 的表达式为第12页，共61页。分段二次插值：若 , 把相邻三个插值节点组成一个插值子区间 ，则整个插值区间 被划分为 个子区间。在第 个子区间上 ，共有三

9、个插值节点 ， 为一二次插值多项式，表达式为例1 已知某函数的函数表如下： 用线性插值法估算 的近似值.第13页，共61页。解：由于 在插值节点 之间，故依此二点构造Lagrange线性插值多项式，并代入 得即 的近似估计值为2.4414.例2 已知观测数据如表3-3所示 试用二次插值方法求 处的插值.第14页，共61页。解： 取包含 的三个观测点 作为插值节点，作二次插值，并令 ,由(3.8)式，可得 =1.8903第15页，共61页。3.2 曲线拟合3.2.1 问题提出 利用插值方法求多项式函数作为未知函数的近似时，要求 1、所有插值节点互不相同，否则不可解； 2、近似函数曲线必须通过所有

10、观测点。 在实际观测或实验中，一般存在以下问题 1、为了得到更加准确、合理的观测结果，经常进行多次重复观测，插值节点互不相同的要求已不成立； 2、由于在观测过程中，常存在许多随机因素，如身高、体重的测量，受测量设备精度、发型、服装、站立方式等影响，测量结果不可避免地存在误差，甚至由于某些因素，误差很大。因此在考虑观测误差的因素下，要求近似函数曲线一定通过观测点已显得没有必要。 因此，只要要求近似函数在观测点上近似地满足插值条件，并使它们的整体误差最小就可以了。第16页，共61页。3.2.2 基本概念 给定函数 (未知)在观测点 上的观测值 ，寻求一近似函数(拟合曲线) ，使在所有观测点上，观测

11、值与近似函数的计算值之间的误差 总体上尽可能接近零，即要求 尽可能反映给定数据点的总体趋势，这就是函数逼近法，也称为曲线拟合法， 称为逼近函数或拟合函数，曲线 称为拟合曲线. 拟合函数的选择范围很广，如多项式，有理函数、指数函数、对数函数、三角多项式等，但具体选择何种函数，应综合多方面因素斟酌确定。比较简单和直观的方法是通过绘制观测点的散点图，进行观察、比较、猜测，然后根据观测结果和误差分析加以确定。第17页，共61页。3.2.3 最小二乘拟合方法 判断拟合曲线 尽可能逼近给定数据点 的标准有很多，如使最大偏差达到极小，所有偏差的绝对值之和取极小等，但因求解方法上的复杂性，实际使用起来并不方便

12、，实践中常用的一种曲线拟合方法就是最小二乘拟合方法. 对给定的数据点 ，选取拟合函数 ，使偏差 ， ，的平方和为最小，即： （3.9） 从几何意义上讲，就是求拟合曲线，使在给定的点 处，计算值与实际观测值的差的平方和最小。这种求近似函数的方法称为离散数据曲线拟合的最小二乘法，函数 称为这组数据的最小二乘拟合函数. 拟合程度的好坏，可以通过直接计算误差平方和的大小来反映。第18页，共61页。若拟合函数取为多项式，如取 次多项式 （3.10）则相应的最小二乘拟合问题就变为：求参数 , 使 （3.11）达到极小。由于拟合函数关于待定系数是线性的，故称该问题为线性最小二乘问题，又由于拟合函数是多项式故

13、也称多项式拟合问题。 最简单的多项式拟合是 和 两种情形，分别称为一次多项式拟合（直线拟合）和二次多项式拟合（抛物线拟合）。以下就这两种简单情形，进行详细讨论。第19页，共61页。3.2.4 拟合直线 1数学描述给定一组散点图呈线性变化的数据 设拟合函数 ，根据最小二乘准则（3.11）知，要求待定系数 ,使达到极小。由多元函数极值定理，一个必要条件是两个偏导数 等于零，整理得 （3.12）把给定数据 全部带入，用消去法解方程组得出这就是线性最小二乘拟合的数学表示，方程（3.12）称为正规方程.第20页，共61页。引入矩阵记号：则正规方程（3.12）可以表示为 （3.13）正规方程的矩阵表示可以

14、推广到 次多项式 的最小二乘拟合中，给定的数据 ，此时的正规方程用矩阵表示为式（3.13），对应的矩阵为： 第21页，共61页。利用正规方程求拟合多项式的系数，可分为如下几个步骤1、输入观测向量 2、生成超定矩阵G3、求4、计算 5、求解参数向量A。例3 观察下表数据，根据其分布趋势，进行曲线拟合，使拟合曲线和下列数据点间的偏差平方和达到最小.第22页，共61页。解：首先作出数据的散点图，如图3-2 所示根据上图，数据点近似直线，所以取 作为拟合函数根据正规方程求解步骤，分别计算如下(1)输入观测向量：第23页，共61页。 （2）生成超定矩阵G（3）求矩阵的转置（4）计算系数矩阵与右端向量（5

15、）把正规方程重新写成线性方程组的形式求解这个方程组，得 a=1.1, b=1.0. 即.第24页，共61页。3.2.5拟合幂曲线1数学描述给定一组散点图呈幂曲线变化的数据 ，设拟合函数 ， 为固定数，根据最小二乘准则（3.11）知，要极小化 最优化的必要条件为导数 等于零，即整理得第25页，共61页。于是，最小二乘拟合幂曲线的正规方程用矩阵表示为 （3.14）其中 ， ,解得 ，其中 是幂指数.例4 观察下表数据，根据其分布趋势，进行曲线拟合，使拟合曲线和下列数据点间的偏差平方和达到最小，并预测x=2.25时y的值.第26页，共61页。解：首先利用(x,y)作散点图，如图3-3 (a) 所示，

16、发现大概呈二次抛物线形式，为验证结果，再利用(x2,y)作图, 如图3-3 (b) 所示，二者呈良好的直线关系。故可以二次幂函数作为拟合函数。取拟合函数形式为，根据（3.14）式，给出求解过程如下： (a) (b)第27页，共61页。（1）输入并形成观测向量 (2) 生成超定矩阵G第28页，共61页。(3)计算 (4) 求参数a: 将x=2.25，代入 得y=16.1193第29页，共61页。3.3利用Matlab求解插值与拟合问题3.3.1 利用Matlab 求解一维插值问题 Matlab用于求解一维插值问题的基本函数为 interp1，给出了四种插值方法供选择，分别为: 最邻近点插值（零次

17、多项式插值）、线性插值、三次插值和三次样条插值，其调用的基本格式为: Interp1(x, y , Cx，Method)其中x, y 分别表示为数据点的横、纵坐标向量, x 必须单调.Cx 为需要插值的横坐标数据(或数组), Cx不能超出x的范围. Method 为可选参数, 对应于上述四种方法, 可从以下四个值中任选一个: nearest -最近邻点插值 linear -线性插值 spline-三次样条插值 cubic-三次插值 linear是缺省值. 当忽略选项时默认为线性插值.第30页，共61页。3.3.2 利用Matlab 求解曲线拟合问题 1 利用正规方程求解，可先把正规方程化为标准

18、矩阵代数方程的形式，然后利用 Matlab命令 x=Ab 2 利用多项式拟合函数polyfit( ) 求解： 调用格式： A=polyfit(x,y,m) 格式说明：x, y 分别对应观测点和观测值，是两个同维向量，m 是拟合多项式的次数；返回结果为所拟合的m次多项式的系数向量，排列次序为由幂次从高到低。 若要利用所求的拟合多项式估计某些点的函数值，可利用多项式求值函数 polyval( ) 调用格式： f=polyval(A,cx) 格式说明：A 为多项式函数的系数向量，按幂次由高到低排列，cx 为要估计的点，可以是一个值，也可以是一个向量。 第31页，共61页。3.3.3 应用示例例5 温

19、度预测问题在12小时内，每隔1小时测量一次温度，详细数据如下：（1）试分别用线性插值、三次样条插值估计在3.2、6.5、7.1、11.7h的温度值，用三次样条插值法每隔0.1h估计一次温度值并画出其图形.（2）用多项式拟合，估计在3.2、6.5、7.1、11.7h的温度值.第32页，共61页。问题分析： 已知一组数据，要估计插值点的值，既可以用插值法，也可以用数据拟合法，拟合出函数关系式，再求给定点的函数值（拟合值）.模型求解： 插值方法： 分别利用线性插值和三次样条方法估计, 编辑Matlab程序如下t=1 2 3 4 5 6 7 8 9 10 11 12;T=5 8 9 15 25 29

20、31 30 22 25 27 24;t0=3.2 6.5 7.1 11.7;T0=interp1(t,T,t0,linear)T1=interp1(t,T,t0,spline)t1=1:0.1:12;T2=interp1(t,T,t1,spline);第33页，共61页。plot(t,T,x);hold onplot(t1,T2,-);legend(观测值,拟合值);xlabel(时间);ylabel(温度);运行结果为：T0 =10.2000 30.0000 30.9000 24.9000T1 =9.6734 30.0427 31.1755 25.3820每隔0.1h估计一次温度值并作图结果

21、见图3-4第34页，共61页。（2）用三次多项式拟合，编程如下：t=1 2 3 4 5 6 7 8 9 10 11 12;T=5 8 9 15 25 29 31 30 22 25 27 24;plot(t,T,*);hold onp=polyfit(t,T,3); % 求三次拟合多项式的系数t0=3.2 6.5 7.1 11.7; % 插值点T0=polyval(p,t0) % 利用所求三次多项式，估计在插值点的值 plot(t0,T0,x); % 绘制观测点散点图hold ont1=1:0.1:12;T2=polyval(p,t1);第35页，共61页。plot(t1,T2,-) %T0为拟

22、合值title(三次多项式拟合)legend(观测值,拟合值,拟合曲线)xlabel(时间)ylabel(温度)运行结果得，利用三次多项式拟合在3.2、6.5、7.1、11.7h的估计值为： T0=14.8017 26.2500 27.3088 23.6551输出图形如图3-5所示。第36页，共61页。例 6 土豆产量与化肥的关系在农业生产试验研究中，对某地区土豆的产量与化肥的关系做了一实验，得到了氮肥、磷肥的施肥量与土豆产量的对应关系如下表：根据上表数据分别给出土豆产量与氮、磷肥的关系式.第37页，共61页。问题分析：首先画出土豆产量与氮施肥量的散点图，见图3-6.从图3-6可以看出，土豆产

23、量与氮肥量的关系是二次函数关系，因此可选取拟合函数为： 其中x和y分别为氮肥施肥量和土豆产量，a,b,c为待定系数.第38页，共61页。再画出磷肥量与土豆产量的散点图，见图3-7.从图3-7可以看出：当磷肥施肥量每公顷达到达到100公斤左右时，两侧曲线分别呈明显的线性关系.由此可选取分段的线性函数作为近似函数，即用 的观测点作一线性拟合函数，再用 的点作一线性拟合函数，最后用两个线性函数求出其分界点即可得分段线性函数.第39页，共61页。模型求解：（1）调用多项式拟合函数p= polyfit (x, y, n)及求值函数y= polyval(p, x)，编程如下：x=0 34 67 101 1

24、35 202 259 336 404 471;y=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75;p=polyfit(x,y,2)y0=polyval(p,x);plot(x,y,x,x,y0,-k)legend(观测值,拟合值)xlabel(氮肥施肥量(Kg/ha)ylabel(土豆产量(Kg)title(土豆产量与氮肥施肥量的二次多项式拟合)第40页，共61页。运行该程序得到拟合二次多项式的系数为p=-0.0003 0.1971 14.7416即所求拟合函数为土豆产量与氮肥量的拟合数据与实验数据的比较，如图3-8.第4

25、1页，共61页。（2）调用线性拟合函数A=regress(Y, G)，编程如下：x1=0 24 49 73 98;y1=33.46 32.47 36.06 37.96 41.04;L1=polyfit(x1,y1,1) %求前段线性拟合多项式x2=98 147 196 245 294 342;y2=41.04 40.09 41.26 42.17 40.36 42.73;L2=polyfit(x2,y2,1) %求后段线性拟合多项式z1=polyval(L1,x1); % 估计前段拟合值z2=polyval(L2,x2); % 估计后段拟合值plot(x1,y1,*,x1,z1,-b)hold

26、onplot(x2,y2,+,x2,z2,-b)xlabel(磷肥施肥量(Kg/ha)ylabel(土豆产量(Kg)hold off 第42页，共61页。运行得：L1 = 0.0844 32.0771L2 = 0.0059 39.9685即左侧拟合直线方程为： ，即左侧拟合直线方程为： ，联立求解这两个方程，得两条直线的交点， 100.53，则拟合函数为为便于观察，拟合多项式对观测数据的拟合情况，程序同时利用拟合结果作了对比分析，见图3-9所示。第43页，共61页。第44页，共61页。3.4 汽车刹车距离模型问题分析：反应距离由反应时间和车速决定，反应时间取决于司机个人状况(灵巧、机警、视野等

27、)和制动系统的灵敏性(从司机脚踏刹车板到制动器真正起作用的时间)，对于一般规则可以视反应时间为常数，且在这段时间内车速尚未改变.制动距离与制动器作用力(制动力)、车重、车速以及道路、气候等因素有关，制动器是一个能量耗散装置，制动力作的功被汽车动能的改变所抵消.设计制动器的一个合理原则是，最大制动力大体上与车的质量成正比，使汽车的减速度基本上是常数，这样，司机和乘客少受剧烈的冲击.至于道路、气候等因素，对于一般规则又可以看作是固定的.第45页，共61页。根据表3-1，取制动距离和刹车距离的平均值，分析一般情况下，反应距离、制动距离和刹车距离各自依赖于速率的情况，作图分析如图3-103-12所示

28、图3-10 反应距离与速率 图3-11制动距离与速率 图3-12刹车距离与速率第46页，共61页。由图3-10可以看出：反应距离跟速率成比例，图3-11-图3-12表明制动距离跟速率符合二次多项式，为了找到刹车距离随速率的变化函数，基于以上分析，作如下假设.基本假设（1）刹车距离 等于反应距离 和制动距离 之和；（2）反应距离 与车速 成正比，比例系数为反应时间 ；（3）刹车时使用最大制动力 ， 做的功等于汽车动能的改变，且 与车的质量 成正比.模型建立由假设（2）可得由假设（3），在作用力 下行驶距离 做的功，所做的功使车速从 变成0，动能的变化为 ，有第47页，共61页。又 ，按照牛顿第二

29、定理可知，刹车时的减速度 为常数，于是其中 ，由假设（1），刹车距离为即实际刹车距离的拟合多项式为（通常对人们刹车的反应时间经验估计为0.75秒，即拟合函数中 .）模型求解根据最小二乘拟合准则（3.11），要拟合参数 ，就要极小化第48页，共61页。最优的一个必要条件是两个偏导数 等于零，从而得到矩阵表示的正规方程为其中为了便于分析，我们对表3-1的数据做了适度整理，见表3-8所示。第49页，共61页。第50页，共61页。根据表3-8的数据，分别由制动距离、反应距离拟合参数 ，根据拟合幂曲线中式（3.14）拟合 ，编程如下：v=29.3 36.7 44 51.3 58.6 66 73.3 80

30、.7 88 95.3 102.7 110 117.3;d1=20 28 40.5 52.5 72 92.5 118 148.5 182 220.5 266 318 376;g=v.*v;k1=gd1;得到 k1 =0.0252d2=22 28 33 39 44 50 55 61 66 72 77 83 88;g=v;k2=gd2;得到 k2=0.7528所以 ，作图：第51页，共61页。K=k1 k2;d=v.*v v*K;plot(v,d,:);hold onplot(v,D,o);legend(d,D);title(plot of and D); 图3-13 拟合图第52页，共61页。由拟合函数计算刹车距离，与实际距离列表比较如下：第53页，共61页。模型检验为衡量上述拟合函数的优劣，用excel计算各点偏差如表3-10所示：第54页，共61页。根据式(3.11)，偏差平方和为 ，均方误差为 英尺小于一个车身的平均长度15英尺.这说明模型结果比较符合实际，曲线拟合的比较好.模型应用 在驾驶员培训班上经常给出