第六章插值与拟合

1、n9nXnnan4XnnXonX1n 4Xon 4X1=(an,an = ,，ao)T，丫 二(yo,yn)T(1)(3)第六章插值与拟合在实际中，常常需要确定一个变量依存于另一个或更多的变量的关系，即函数。但实际上确定函数的形式(线性形式、乘法形式、幕指形式或其它形式)时往往没有先验的依据。 只能在实验或测量得到的离散数据的基础上进行试验，这时候我们常常会用到插值与拟合的方法。插值与拟合就是要通过已知的数据去确定某一类已知函数的参数或寻找某个近似函数， 使所得到的近似函数对已知数据有较高的拟合进度。如果要求这个近似函数经过所有已知的数据点，则称此类问题为插值问题。当所给的数据较多时， 用插值

2、方法所得到的插值函数会很复杂，所以，通常插值方法用于数据较少的情况。其实，通常情况下数据都是由观测或试 验得到的，往往会带有一定的随机误差，因而，要求近似函数通过所有的数据点也是没有必要的。如果不要求近似函数通过所有的数据点，而是要求他能较好地反映数据的整体变化趋势，则解决这类问题的方法称为数据拟合。虽然插值与拟合都是要构造已有数据的近似函数，但因对近似要求的准则不同，因此二者在数学方法上有很大的差异。第一节一般插值方法1. 问题的提出插值问题的一般提法：已知n -1个节点(xj, yj) (j二0,1/ n)，其中xj (j = 1,2 - , n)互不相同，要求构造一个函数y = f(x)

3、使得y f (Xj) (j =0,1,n)。我们通常称这样一类问题为插值问题，并称构造的函数y = f (x)为插值函数，xj (j = 1 , 2 , n为插值节点，yj = f (xj) (j = o , 1 n 为插值条件。2. 多项式插值从理论和计算的角度看，多项式是最简单的函数，设f (X)是n次多项式，记作Ln(x) =anXn an"2 d x a°对于节点(Xj,yj)应有Ln (Xj ) = yj , j - 0,1,2，,n为了确定插值多项式Ln(x)中的系数a.,a.亠，c ,a°，将(1 )代入(2), 有an Xo ' a“ 4X

4、0 ' a1Xo a - yoanX； - an4X；a a。= %1X11X1n 4xn方程组(3)简写成1X1(4)XA=Y注意：detX是Vandermonde行列式，利用行列式性质可得detX ：| 丨区-Xj)0 勺:k m因xj互不相同，故detX = 0，于是方程(4 )中A有唯一解，即根据n 1个节点可以 确定唯一的n次插值多项式。3.拉格朗日(Lagrange)插值多项式实际上比较方便的做法不是解方程(4)求A,而是先构造一组基函数：(5)li(X)=(X-X0)(X-XX-Xii)(X-Xn).。佔(Xi X。)(Xi Xi_L)(X x 卅)(Xi -Xn)|j(

5、x)是n次多项式，满足li(Xj)o,1, j i, j =0,12 ,ni = j(6)(4)(4)n(7)Ln(x) yli(x)显然Ln(X)是满足(2)的i =0n次多项式，由方程(4)解的唯一性，(7)式表示的Ln(x)的解与(1)式相同。(5)、( 7)称拉格朗日插值多项式，用Ln(x)计算插值称拉格朗日多项式插值。4.牛顿(Newton)插值构造门次多项式 Nn(x) = f (X0) f(X°,X1)(X -X0)f(X°,X1,X2)(X-X0)(X-xj+f (X。，,，Xn)(X x°)(x Xj(X X称该多项式为牛顿插值多项式，其中上f

6、(x0) f (xjf(Xo,XJ01(二个节点，一阶差商)X。一X1f(X0,X1)- f(X1,X2)f(X°,X1,X2)0 -(二个节点，二阶差商)X0 X2f (Xo,%，,Xn4)- f (%，焉)f(X0,X1,.,Xn)0 1 口1 已 (n+1 个节点，n 阶差商)X0 Xn实际上，牛顿插值公式是拉格朗日插值公式的一种变形，二者是等价的。另外还有著名的埃艾米特(Hermite )插值等。5.分段线性插值简单地说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插 值函数，记作ln(X)，它满足ln(Xj)=yj，且ln(X)在每个小区间Xj ,Xj 1

7、上是线性函数 (j =01, n)。I n (x)可以表示为(12)ln(x)yjlj (x)j=0s»XjXj _ Xj 二(j 0时舍去)x x 丄(13)l j (x) = , Xj 兰 x 兰 Xj 卅 (j = n时舍去)Xj _Xj 卅0,其它In(x)有良好的收敛性，即对于 Xa,b有，lim ln(x)=g(x)。n >:用ln(x)计算X点的插值时，只用到 X左右的两个节点，计算量与节点个数n无关。但n越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。Matlab中分段线性

8、插值有现成的程序：y=interp1(x0， y0，x)其中输入x0、y0、x分别表示节点数据(xj,yj) (j = 0,1/ n)和插值点以数组，输出y为插值点对应的插值。数组长度自定义(x0和y0同长度，x和y同长度)。第二节样条函数插值方法1. 样条函数的由来分段线性插值虽然简单， n足够大时精度也相当高。 但是折线在节点处显然不光滑， 即 In(x)在节点处导数不连续。这影响了它在诸如机械加工等领域(希望插值曲线光滑)中的应用。所谓样条(Spline)，来源于船舶、飞机等设计中描绘光滑外形曲线用的绘图工具。一根 有弹性的细长木条用压铁固定在节点上，其它地方让它自然弯曲，如此画出的曲线

9、称为样条曲线。因为这种曲线的曲率是处处连续的，所以要求样条函数的二阶导数连续。人们普遍使用的样条函数是分段三次多项式。2. 三次样条函数三次样条函数 记作S(x),a空x乞b。要求它满足以下条件：a) 在每个小区间Xjx/i =1,，n)上是3次多项式；b) 在a乞x乞b上二阶导数连续；c) S(n)二 yj =0,1, ,n。(14)由条件a,不妨将S(x)记为S(x) =<S (x),x Xi：,Xi,i =1,，nS(x) =ajX3 +bjX2 +CjX + di( 15)其中ai ,bi,Ci,di为待定系数，共4n个。由条件b，S(Xi) =Si41(Xi)*S：(Xi)=S

10、i；1(Xi)i =1,2，n-1( 16)0(xi) = SE(Xi)容易看出，(14)、(16)式共含有4n-2个方程，为确定 S(x)的4n个待定参数，尚需再 给出2个条件。最常用的是所谓自然边界条件：S”(X°)=S“(Xn)=O( 17 )可以证明，4 n阶线性方程组(14)、(16)、(17)有唯一解，即S(x)被唯一确定。但是，这 种解法的工作量太大，方程组又常呈病态，所以实际上要设计简便的解法。Matlab中三次样条插值有现成的程序：y=interp1 ( xO， yO, x， 'spline') 或 y=spli ne(x0,y0,x)其中输入x0、

11、y0、x分别表示节点数据(xj, yj) (j =0,1，n)和插值点以数组，输出 y为插值点对应的插值。数组长度自定义(x0和y0同长度，x和y同长度)。另外，像分段线性函数ln(X)样，三次样条函数S(x)也有良好的收敛性，即在相当一般的条件下，lim S(x)二g(x)。n_jpc3. 插值方法小结拉格朗日插值是高次多项式插值(n+1个节点上用不超过n次的多项式)，插值曲线光滑，误差估计有表达式。但有振荡现象，收敛性不能保证。这种插值主要用于理论分析，实 际意义不大。分段线性和三次样条插值是低次多项式插值，简单实用，收敛性有保证，但不光滑，三次样条插值的整体光滑性已大有提高，应用广泛，唯

12、误差估计较困难。第三节Matlab在插值问题中的应用1拉格朗日插值 拉格朗日插值没有有现成的程序，需要编译实现。1.1 Matlab 实现函数Lagrange.mfunction y=lagra nge(x0,y0,x) n=le ngth(x0);m=le ngth(x); for i=1:mz=x(i);s=0.0;for k=1: np=1.0;for j=1: n if j=k p=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;end y(i)=s; end1.2应用：x=1 2 3 4; y=0 -5 -6 3;lagra nge(x,y,2.5

13、)ans =-6.37502. Runge现象及分段线性插值Runge 现象Runge在本世纪初发现：在-1,1上用n 1个等距结点作插值多项式Pn(x)，使其在各结点的值与函数 y(x) =1/(1 25x2)在结点的值相等。但在n 时，插值多项式R(x)在区间中部趋于y(x)。但对于0.726 <1 x I < 1的x , Pn(x)严重发散。2通过下面的例子，以图形的方式体会Runge现象(令y(x) = 1/(1 X )x=-5:1:5; y=1./(1+x.A2); x0=-5:0.1:5;y0=lagra nge(x,y,x0); y1=1./(1+x0.A2);%绘制

14、图形plot(x0,y0,'-r')hold on22.2 Matlab实现分段插值维插值interp1yi=i nterp1(x,y,xi)对(x,y)进行插值，计算插值点xi的函数值yi=in terp1(y,xi)默认x=1:n , n是向量y的元素个数plot(x0,y1,'-b')2.3例题用一维线性插值解决Runge现象x=0:0.1:10;yi=interp1(x,y,xi, 'method ' 指定特定算法插值，method可以是如下字符串lin ear线性插值spli ne三次样条插值cubic三次插值要求：x是单调，但不要求连

15、续等距。如果x连续等距，可以选用快速插值法。调用函数时只需在 method前加” ”，女口”spliney2=i nterp1(x,y,x0);plot(x0,y2,'*m')T0.80.60.400.2-0.2-0.4-0.6-0.8正弦曲线的插值示例y=si n( x);xi=0:0.25:10;yi=i nterp1(x,y,xi);plot(x,y, ' ,xi,yi)3. 三次样条插值y3=i nterp1(x,y,x0,'*spli ne');y3=spli ne(x,y,xO);plot(x0,y3,'-g')第四节数据拟合

16、1.方法介绍在实际生活中，往往需要从一组实验数据(xj,yj) (j =0,1,n)中寻找出变量x , y之间的函数关系。由于观测数据不可避免出现误差，因此并不需要y二f (x) 一定要经过所有 的点，而只要求在给定点 Xj上误差 d = f(Xj)-yj按某种标准达到最小。 通常用欧式范数2I -l作为误差量度的标准。这就是所谓的最小二乘法。注意：数据拟合与插值的最大区别在于拟合需要给出一个曲线方程的具体解析形式, 插值只需求出该点的内插数值。2 线性拟合线性拟合以最简单的一次线性方程f(x)=a 1X+a°拟合数据。按最小二乘法，a1 ,a 0需满足n2R八，(yi -(a。 a

17、1x)最小，因此可以通过i 4dR« ca0kca1-0=0求出此时的a1,a0 。事实上，在 MATLAB中已有现成的求最小二乘问题的函数polyfit，称为多项式拟合函数，并且这个函数允许多项式的次数可以是任意次的。除外，还可以用解线性方程组中的除法运算(矩阵除法)来求解。这两个方法的区别在于：用polyfit函数求拟合问题时，多项式的次数必须从0次到最高次数n之间每个次数都要出现。而如果需要选择一些次数进行拟 合时，就可用矩阵除法运算来进行。矩阵除法还可以求一般的线性拟合问题，例如拟合函数不是多项式的线性拟合问题。多项式拟合是线性拟合问题(注意：无论拟合的多项式次是多少，多项式

18、拟合都是线性 拟合。)但在实际应用中，有时还需要作非线性拟合问题。所谓线性拟合问题是指：需要拟合的函数中的未知常数都线性的。如函数y = a bx2中，常数a,b是线性的。但y = ax ebx、y = aebx中的常数a,b都是非线性。这种函数的 拟合问题称为非线性拟合问题。有的非线性拟合问题可以化为线性拟合问题。3 可线性化的非线性模型模型形式变换后形式变量和参数的变化YXa1a2ax1 1b +111by 一1 +bxy axayxaaaV 1 xb1x1by .x by aayaaaxxb22 xy21b2y 一 22b - xy aaxxaaby = axIn y = b 1 n x

19、 + 1 n aln yln xln abbxy = aeIn y = bx +ln aln yxln ab"b2y = ae2 x ln y =b-+ln aln y2 xln a1b22 2a2b2y2 b2b2 x22 y2 xb2b2y _ b _2 x a2 a第五节Matlab在拟合问题中的应用1线性拟合及多项式拟合120以最高次为i的多项式拟合数据点(x,y)1008060402000.511.522.533.544.550Jployfit(x,y,i)例 1 试验 ployfit(x,y,i)x=0 1 2 3 4 5;y=0 21 62 70 77 110; coe

20、f=polyfit(x,y,1);a仁coef(1),a0=coef(2); ybest=a1*x+aO;s=sum(y-ybest).A2); axis(-1,6,-20,120);plot(x,y, '*')hold onplot(x,ybest)例2如下给出从二阶到十阶多项式拟合曲线的比较程序，并给出拟合曲线x=0 1 2 3 4 5;2例3用最小二乘法求一个形如 y =a bx的经验公式，数据如下:x1925313844y19.032.349.073.398.8y=0 21 62 70 77 110;xi=0:0.2:5;for n=2:10bb=polyfit(x,

21、y,n);yi=polyval(bb,xi); plot(xi,yi,x,y, '* ')title(int2str(n),'次多项式拟合曲线') grid on pauseend120解用求矩阵除法(因为要拟合的多项式缺了1次幕项，所以不能用polyfit函数)。x=19 25 31 38 44;y=19.0 32.3 49.0 73.3 98.8;x1=x.A2;x1=ones(5,1),x1'yx0=19:0.2:44;y0=ab(1)+ab (2)*x0.A2;plot(x,y, 'o')hold onplot(x0,y0, &#

22、39;-r')例4在区间-1,3内拟合函数y = ax ebx。解用非线性拟合函数curvefit来拟合。先建立拟合函数。%建立拟合函数，文件名是nxxyhhx.m，必须与函数名相同。%要拟合的函数中参数用x表示，即x(1)=a x(2)=b ；%而拟合函数中x的值则用xdata表示。function v=n xxyhhx(x,xdata)v=x(1)*xdata+exp(x(2)*xdata);以下指令在命令窗中进行。clf;x=li nspace(-1,3,10);y1=2*x+exp(-0.1*x); % 原型函数plot(x,y1,'-k')hold ony=y

23、1+1.2*(ra nd(size(x)-0.5);% 将原型函数加一些扰动plot(x,y,'*g')x0=2.5,-0.5;a=curvefit(' nxxyhhx',x0,x,y)% 用原始实验数据拟合函数n xxyhhx (x),vpa(a(1),a(2),8)% nxxyhhx (t)表达式中各项的系数。y2=n xxyhhx(a,x);plot（x,y2,'）lege nd（'原型函数，原始数据','用原始数据拟合的结果,4）;第六节建模实例1.血液流量问题小哺乳动物与小鸟的心跳速度比大哺乳动物与大鸟的快。如果动物的进

24、化为每种动物确定了最佳心跳速度，为什么各种动物的最佳心跳速度不一样呢？由于热血动物的热量通过身 体表面散失，所以它们要用大量的能量维持体温，而冷血动物在休息时只需要极少的能量， 所以正在休息的热血动物似乎在维持体温。可以认为，热血动物可用的能量与通过肺部的血液流量成正比。（1）试建立一个模型，将体重与通过心脏的基础（即休息时的）血液流量联系起来， 用下面的数据检验你的模型。（2）有许多可得到脉搏数据但没有血液流量数据的动物，建立一个模型将体重与基础 脉搏联系起来，用下面的数据检验你的模型。（3）在检验你在（1 ）和（2）中的模型时会出现不一致，试进行分析。鸟类体重（克）脉搏 （次 /分）蜂鸟4

25、:615鹪鹩11450金丝雀16514麻雀28350鸽子130135表五关于哺乳动物的数据鸟类体重（克）脉搏 （次 /分）海鸥388401鸡1980312秃鹰8310199火鸡875093驼鸟8000065表一关于某些哺乳动物的数据哺乳动物名称兔山羊狗狗狗体重（千克）4.12416126.4基础血液流量（分升/分）5.331十221211表二关于人类的数据年龄5101625334760体重（千克）18316668707270基础血液流量（分升/分）23335251434046脉搏（次/分）96906065687280表三 关于小鸟类的数据表四 关于大鸟类的数据哺乳动物名 称体重（千克）脉搏（次

26、/分）哺乳动物名 称体重（千克）脉搏（次/分）小蝙蝠0.006588海豹2025100小家鼠0.017500山羊3381仓鼠r 0.103347绵羊r 50r 7080小猫0.117300猪1006080大家鼠0.252352马P 3804503455天竺鼠0.437269牛5004653兔1.34251象200030002550这里只对该问题作一些拟合方面的练习。其它问题读者可自己进行讨论。 符号用W表示动物的体重，单位：千克用V表示动物的基础血液流量，单位：公升/分用t表示动物的年龄，单位：岁用n表示动物的脉搏，单位：次 /分假设动物的基础血液流量与动物的体重之间存在一定的函数关系V =

27、f(W)，可以用表一中的数据来拟合这个函数。函数f (w)是一个什么样的函数呢？由于我们对“动物的基础血液流量与动物的体重” 之间的关系并不清楚，所以只有根据表一中的数据得出函数f (w) 一些性质。先将表一中的数据用 MATLAB 软件作出图形。从图上可以看出，这个函数关系 V = f(w)应当是一个单调增加的函数。因此，拟合的函数如果不具有这一性质的话，就不 能作为是好的选择。一般地，可以假设函数f(w)是一个多项式，通常，这个多项式的次数不要超过3、4次，具体可根据拟合的效果来定。当然也可以用其它函数来拟合。为了提高拟合的效果，函数f(w)还可以用分段函数来拟合。以下是用分段函数拟合的结

28、果：'-0.0033W4 +0.1706W3 2.9727W2 +21.3418W-43.0649,v = f (w) = «w壬4.1,16-0.204756W2 +9.31526W-74.6265,w 16,24拟合函数图形是：基础血液流量与动物的体重之间的函 如曲 原始数据图形问题1写出拟合函数 V=f(w)和作出上面图形的 MATLAB指令。同样可以拟合人的基础血液流量与体重之间的函数关系V = g(w)，可以用表二中的数据来拟合这个函数。这里用4次多项式来拟合，拟合的结果是：V = -0.000025W40.00351W3 -0.1728W24.340w-16.98

29、,w 18,72拟合函数图形是：人的基砒血液流量与体重之间的函数丫利帥和原始数据图形问题2:写出拟合函数 V二g(w)和作出上面图形的 MATLAB指令。将上面拟合出来的函数 v = f(w)和v =g(w)在它们的公共定义域18,24上的图形画 出来，如下图所示。从图形上可以看出人类与动物之间的差异。人和动物在体重的公共区间上黄于基础血液流量函数的图形问题3:写出作出上面图形的 MATLAB指令。 下面考虑动物、人类的体重与基础脉搏的函数关系。假设人类的体重与基础脉搏之间的函数关系是w = Hi(n)，利用表二中的数据来拟合这个函数。这里用 3次多项式拟合。拟合的结果是：wnHjn) =0.

30、0002566? -0.143347n216.2450n-450.216, n 60,96其图形是：人的基础脉撐与体重之间的函数庐片伸闲I原始数据图形Q I1IIII60 -40+ n=c+c+原始数据° 1160657075 冊 669095100问题4:写出拟合函数 wH"n)和作出上面图形的 MATLAB指令。假设哺乳动物的体重与基础脉搏之间的函数关系是w二H2(n)，利用表五中的数据来拟合这个函数。这里用分段函数来拟合。由于当 n : 100时，w变化激烈，所以用多项式已不能描述其变化的规律，可用其它函数来拟合。拟合的结果是:-0.0345386n 14.7835,

31、 n _ 100w = H2 (n) = < 7.56493汉 1012 “6,n <100l n 图形如下。哺乳动物的基础脉搏与体重的函数w=H2(nn原始数据图形3000200010Q00100200300400500600n问题5:写出拟合函数 w = H2(n)和作出上面图形的 MATLAB 指令。 假设小鸟类、大鸟类的体重与基础脉搏之间的函数关系分别是w = H 31(n)和-1000w = H32(n)，利用表三、四中的数据来拟合这两个函数。拟合的结果是：w 二 H31(n) - -0.0000021n3 0.0031454n2 -1.61275n 295.685 n

32、135,615102.04631 10n315050°0n +00 +c,m+c4+原始数据20040060080011X,n 匸65,401 其图形如下。鸟类基础脉搏与体堇N间的函数护巴”仆rH豊和原始数图形问题6:写出拟合函数 w二H 31 (n)和w = H 32 (n)和作出上面图形的 MATLAB指令。 下面考虑人类的基础血液流量与基础脉搏的函数关系。假设人类的基础血液流量与基础脉搏之间的函数关系是v =U (n),利用表二中的数据来拟合这个函数。拟合结果是：v =U(n) - -0.00217132n30.482394n2 -37.3316n 989.322,n 60,9

33、6人类的基础血液济量与基础脉搏之间的函数v=U(n丙口原始数据图形60 111111rII30Fc+n+cn+c+原始数据1002Q LiI丄1l_6065707580859095n考虑复合函数V = gH1( n) =-.2489e-4*(.256557e-3* 门人3-.143347* 门人2+16.2450* n-450.216)M+.3506e-2 *(.256557e-3* 门人3-.143347* 门人2+16.2450* n-450.216)A3-.1728*(.256557e-3* 门人3-.143347* n2+1 6.2450* n-450.216)A2+.1113457380e-2* 门人3-.622125980* 门人2+70.5033000* n-1970.917440;用MATLAB软件画出上面两个函数的图形：人类的基础血液流量与基础脉搏之间的函数戶U