插值与数据拟合模型

1、5)6)5)6)第二讲插值与数据拟合模型函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟用插值还是拟合，有时容易确定，有时则并不明显。在数学建模过程中，常常需要确定一个变量依存于另一个或更多的变量的关系，即函数。但实际上确定函数的形式(线性形式、乘法形式、幂指形式或其它形式)时往往没有先验的依据。只能在收集的实际数据的基础上对若干合乎理论的形式进行试验，从中选择一个最能拟合有关数据，即最有可能反映实际问题的函数形式，这就是数据拟合问题。一、插值方法简介插值问题的提法是，已知n+1个节点(x,y),j二0,1,2

2、,n,其中x互不相同，不妨设a二xxx二b，求任一插值点x*(鼻x)处的插值y*。(x,y)可以看成是由某个函数01njjjy=g(x)产生的，g的解析表达式可能十分复杂，或不存在封闭形式。也可以未知。求解的基本思路是，构造一个相对简单的函数y二f(x)，使f通过全部节点，即f(x.)二y.(j二0，1，2,n)，再由f(x)计算插值，即y*二f(x*)。1拉格朗日多项式插值插值多项式从理论和计算的角度看，多项式是最简单的函数，设f(x)是n次多项式，记作1)L(x)=axn+axn-1hfax+annn-110对于节点(x,y)应有.L(x)=y,j=0,1,2,nnjj为了确定插值多项式L

3、(x)中的系数a,a，,a,a，将(1)代入(2)，有2)方程组(3)简写成nn-110axn+axnhfax+an0n-1010(axn+axnhfax+an1n-11110=y00=y13)axn+annn-1nxnhhaxha1n0 xnxn-100 xnxn-111jxnxn-1nn,A=(a,ann-1,a0)T,Y=(y0,y1，,yn)T注意detX是Vandermonde行列式，XA=Y利用行列式性质可得detX二H(xk0jkn4)因x互不相同，故detX丰0，于是方程(4)中A有唯一解，即根据n+1个节点可以确定唯一的n次插值多项式。拉格朗日插值多项式实际上比较方便的做法不

4、是解方程(4)求A,而是先构造一组基函数：、(x-x)(x-x)(x-x)(x-x)l(x)=01i+4n-il(x)是n次多项式，满足i:i-i+:n,i=0,1,2,n(x-x)(x-x)(x-x)(x-x)i0ii-1iih1inl(x)=ri=ji,j=0,1,2,nij|0,心j1/14 /14 /147)L(x)=y1(x)niii=o显然L(x)是满足(2)的n次多项式，由方程(4)解的唯一性，(7)式表示的L(x)的解与(1)式相nn同。(5)、(7)称拉格朗日插值多项式，用L(x)计算插值称拉格朗日多项式插值。n误差估计插值的误差通过插值多项式L(x)与产生节点(x,y)的g

5、(x)之差来估计，记作R(x)。虽然我们可njjn能不知道g(x)的解析表达式，但不妨设g(x)充分光滑，具有n+1阶导数。利用泰勒展开可以推出，对于任意xea,b。(n+1)!R(x)二g(x)-L(x)二-肝(x-x覧w(a,b)nn(n+1)!j若可以估计j=oIg(n+1)(g)|Mn+1实际上因为MIR(x)I减少；nn+1例将区间o,MnIR(x)I=+nIxxI,ge(a,b)n(n+1)!jj=o常难以确定，所以(10)式并不能给出精确的误差估计。但是可能看出，n增加,10)g越光滑，M越小，IR(x)I越小；x越接近x，IR(x)I越小。n+1njnn等分，用y=g(x)=c

6、osx产生n+1个节点，然后作拉格朗日插值多项式。用兀L(x)计算cosn6(取4位有效数字)。估计IR(x)I(取n=1,2)。n兀c),0。由(5)、(7)式丿兀x二+ox0=12x0-1t-0兀22则(%,yo)=(0,1),(x1，y1)=L1(x)=yo1o+”广1-2，则(x,y)=(0,1)，(x,y)=-,0.7071,(x,y)=-,00011v4丿22v2丿=0.6667若n二，由(5)、(7)式。1l6丿兀cos=L6L2(x)=yo1o+y111+y212=1(兀x4八兀Vo4、丫兀x2丿o八2丿(xo)x+o.7o71-兀、2丿兀兀42丿兀(x0)(x)4兀兀24丿兀

7、0丫12人(兀)(1)16(1)xX0.7071xxV4丿V2丿12V2丿1cos=L=0.8508。62V6丿=8兀2估计IRn(x)I:对于g(x)=C0Sx可设Mn+1=】，记节点间隔h=2。当xxj+1hh2IxxIIxxI-jj+14肝丨xx丨-2h-3hnhj4j=0于是(10)式给出1h2hn+1兀n+1丨R(x)|-2h-3h-nhn(n+1)!44(n+1)4(n+1)(2n)n+1可以算出n1234丨R(x)|n0.30.044.7x10-34.7x10-4,L16J216丿的误差在IR(x)l范围内。n兀cos的精确值是0.8660(4位有效数字)L61插值多项式的振荡用

8、拉格朗日插值多项式L(x)近似g(x)(axb)虽然随着节点个数的增加，L(x)的次数变大，nn多数情况下误差IR(x)I会变小，但n增加时，L(x)的光滑性变坏，有时会出现很大的振荡。理论上,nn当nTa时，在a,b内并不能保证L(x)处处收敛于g(x)。Runge给出了一个有名的例子：n1g(x)二,xg5,51+x2取x.=5+10j,j二0,1,2,n。对于n二2,4,6,作L(x)，会得到如下图所示的结果。可以看出,jnn对于较大的|x|，随着n的增加，L(x)的振荡越来越大，事实上可以证明，仅当|x|3.63时，才有nlimL(x)二g(x)，而在此区间外,L(x)是发散的。nnn

9、Ta高次插值多项式的这些缺陷，促使人们转而寻求简单的低次数多项式插值。2.分段线性插值简单地说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记作I(x)，它满足I(x)二y，且I(x)在每个小区间x,x上是线性函数(j二0,1,n)。nnjjnjj+1I(x)可以表示为n12)I(x)=Yyl(x)njjj=0 xxd,xxx(j=0时舍去)xxj1jxx亠，xxjj+10,xxx(j=n时舍去)jj+1其它13)jj1I(x)有良好的收敛性，即对于xga,b有，limI(x)=g(x)。nnnTa用I(x)计算x点的插值时，只用到x左右的两个节点，计算量与节点个数

10、n无关。但n越大，分段n越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。3.三次样条插值样条函数的由来分段线性插值虽然简单，n足够大时精度也相当高。但是折线在节点处显然不光滑，即I(x)在节点n处导数不连续。这影响了它在诸如机械加工等领域(希望插值曲线光滑)中的应用。所谓样条(Spline),来源于船舶、飞机等设计中描绘光滑外形曲线用的绘图工具。一根有弹性的细长木条用压铁固定在节点上，其它地方让它自然弯曲，如此画出的曲线称为样条曲线。因为这种曲线的曲率是处处连续的,所以要求样条函数的二阶导数连续。人们普遍使用的样条

11、函数是分段三次多项式。三次样条函数三次样条函数记作S(x),axb。要求它满足以下条件：在每个小区间x,x(i二1,.,n)上是3次多项式；i-1i在axb上二阶导数连续；TOC o 1-5 h zS(x)二y,i=0,1,n。(14)ii由条件a,不妨将S(x)记为S(x)二(x),xex,x,i二1,.,n)ii-1iS(x)二ax3+bx2+cx+d(15)iiiii其中a,b,c,d为待定系数，共4n个。由条件b,iiiiS(x)二S(x)iii+1iS(x)二S(x)i二1,2,n一1(16)iii+1iS(x)二S(x)iii+1i容易看出，(14)、(16)式共含有4n-2个方程

12、，为确定S(x)的4n个待定参数，尚需再给出2个条件。最常用的是所谓自然边界条件：S(x)=S(x)=0(17)0n可以证明，4n阶线性方程组(14)、(16)、(17)有唯一解，即S(x)被唯一确定。但是，这种解法的工作量太大，方程组又常呈病态，所以实际上要设计简便的解法。另外，像分段线性函数I(x)一样，三次样条函数S(x)也有良好的收敛性，即在相当一般的条件下，nlimS(x)=g(x)。4.用MATLAB作插值计算拉格朗日插值需先按照(5)、(7)式编写一个程序。设n个节点以数组x0,y0输入(注意：程序中用n个节点，而不是(5)、(7)式中的n+1个节点)，m个插值点以数组x输入。输

13、出数组y为m个插值。比如可以写一个名为lagr1.m的M文件。分段线性插值有现成的程序y=interp1(x0,y0,x)其中输入x0,y0,x和输出y的意义同上，数组长度自定义(x0,y0同长度，x,y同长度)。三次样条插值也有现成的程序y=interp1(x0,y0,x,spline)或y=spline(x0,y0,x)其中输入x0,y0,x和输出y的意义同上，数组长度自定义(x0,y0同长度，x,y同长度)。例对y二八I、,-5x5，用n(=11)个节点(等分)作上述三种插值，用m(=21)个插值点(等分)(1+x2)作图，比较结果。插值方法小结拉格朗日插值是高次多项式插值(n+1个节点

14、上用不超过n次的多项式)，插值曲线光滑，误差估计有表达式。但有振荡现象，收敛性不能保证。这种插值主要用于理论分析，实际意义不大。分段线性和三次样条插值是低次多项式插值，简单实用，收敛性有保证，但不光滑，三次样条插值的整体光滑性已大有提高，应用广泛，唯误差估计较困难。二、最小二乘法简介下面先看一个例子。例1“人口问题”是我国最大社会问题之一，估计人口数量和发展趋势是我们制定一系列相关政策的基础。有人口统计年鉴，可查的我国从1949年至1994年人口数据智料如下：年份1949195419591964196919741979198419891994人口数（百万）541.67602.66672.097

15、04.99806.71908.59975.421034.751106.761176.74分析：（1）在直角坐标系上作出人口数的图象。（2）估计出这图象近似地可看做一条直线。（3）用以下几种方法（之一）确定直线方程，并算出1999年人口数。方法一：先选择能反映直线变化的两个点，如（1949，541.67），（1984，1034.75）二点确定一条直线，方程为：N=14.088t-26915.842,代入t=1999,得N-12.46亿。方法二：可以多取几组点对，确定几条直线方程，将t=1999代入，分别求出人口数，再取其算数平值。方法三：可采用“最小二乘法”求出直线方程。最小二乘法简介设（x,y

16、）,（x,y），,（x,y）是直角平面坐标系下给出的一组数据，设xxx，我们可1122nn12n以把这组数据看作是一个离散的函数。根据观察，如果这组数据图象“很象”一条直线（不是直线），我们的问题是确定一条直线y=ax+b，使得它能最好的反映出这组数据的变化。对个别观察值来说，用直线y二ax+b的值来近似代替其观察值时，所产生的误差可能是正的，也可能是负的。为了不使它们相加彼此抵消，可用Iy-（ax+b）Iiii=1来表示用直线y=ax+b来近似代原来实验数据时所产生的误差。为了在数学上处理方便，又把上式改成d=y-（ax+b）2iii=1也就是说，我们选取常数a,b，使得总误差d达到最小。这

17、就是所谓的最小二乘法。用微分法不难求出上面最小值问题的驻点，这里不列出其结果。事实上，在MATLAB中已有现成的求最小二乘问题的函数polyfit，称为多项式拟合函数，并且这个函数允许多项式的次数可以是任意次的。除外，还可以用解线性方程组中的除法运算（矩阵除法）来求解。这两个方法的区别在于：用polyfit函数求拟合问题时，多项式的次数必须从0次到最高次数n之间每个次数都要出现。而如果需要选择一些次数进行拟合时，就可用矩阵除法运算来进行。矩阵除法还可以求一般的线性拟合问题，例如拟合函数不是多项式的线性拟合问题。上面例1中的问题就可以用polyfit来求解。例2用最小二乘法求一个形如y=a+bx

18、2的经验公式，数据如下：x1925313844y19.032.349.073.398.8解用求矩阵除法（因为要拟合的多项式缺了1次幕项，所以不能用polyfit函数）。代码如下。x=1925313844;y=19.032.349.073.398.8;x1=x.A2;x1=ones（5,1）,x1ab=x1yx0=19:0.2:44;y0=ab（1）+ab（2）*x0.A2;clfplot（x,y,o,x0,y0,-r）多项式拟合是线性拟合问题（注意：无论拟合的多项式次是多少，多项式拟合都是线性拟合!）。但在实际应用中，有时还需要作非线性拟合问题。所谓线性拟合问题是指：需要拟合的函数中的未知常数

19、都线性的。如函数y=a+bx2中,常数a,b是线性的。但y=ax+ebx、y=aebx中的常数a,b都是非线性。这种函数的拟合问题称为非线性拟合问题。有的非线性拟合问题可以化为线性拟合问题。例如在函数y=aebx中，两边取对数，得lny=lna+bx，再令a=lna，z=lny，则要拟合的函数就成z=a+bx，这样就变成线性拟合问题了。但也有不能化成线性拟合问题的情况，如函数y二ax+ebx就是这样。在MATLAB5.3中求非线性拟合问题的函数是lsqcurvefit。例3在区间-1,3内拟合函数y二ax+ebx。解用非线性拟合函数lsqcurvefit来拟合。先建立拟合函数。%建立拟合函数，

20、文件名是nxxyhhx.m，必须与函数名相同。%要拟合的函数中参数用x表示，即x(1)=ax(2)=b；%而拟合函数中x的值则用xdata表示。functionv=nxxyhhx(x,xdata)v=x(1)*xdata+exp(x(2)*xdata);以下指令在命令窗中进行。clf;x=linspace(-1,3,10);y1=2*x+exp(-0.1*x);%原型函数plot(x,y1,-k)holdony=y1+1.2*(rand(size(x)-0.5);%将原型函数加一些扰动plot(x,y,*g)x0=2.5,-0.5;a=lsqcurvefit(nxxyhhx,xO,x,y)%用

21、原始实验数据拟合函数nxxyhhx(x),vpa(a(1),a(2),8)%nxxyhhx(t)表达式中各项的系数。y2=nxxyhhx(a,x);plot(x,y2,-r)legend(原型函数，原始数据,用原始数据拟合的结果,4);三、血液流量问题小哺乳动物与小鸟的心跳速度比大哺乳动物与大鸟的快。如果动物的进化为每种动物确定了最佳心跳速度，为什么各种动物的最佳心跳速度不一样呢？由于热血动物的热量通过身体表面散失，所以它们要用大量的能量维持体温，而冷血动物在休息时只需要极少的能量，所以正在休息的热血动物似乎在维持体温。可以认为，热血动物可用的能量与通过肺部的血液流量成正比。试建立一个模型，将

22、体重与通过心脏的基础(即休息时的)血液流量联系起来，用下面的数据检验你的模型。有许多可得到脉搏数据但没有血液流量数据的动物，建立一个模型将体重与基础脉搏联系起来，用下面的数据检验你的模型。在检验你在(1)和(2)中的模型时会出现不一致，试进行分析。表一关于某些哺乳动物的数据哺乳动物名称兔山羊狗狗狗体重(千克)4.12416126.4基础血液流量(分升/分)5.331221211表二关于人类的数据年龄5101625334760体重(千克)18316668707270基础血液流量(分升/分)23335251434046脉搏(次/分)96906065687280表三关于小鸟类的数据表四关于大鸟类的数

23、据鸟类体重(克)脉搏(次/分)鸟类体重(克)脉搏(次/分)蜂鸟4615海鸥388401鹪鹩11450鸡1980312金丝雀16514秃鹰8310199麻雀28350火鸡875093鸽子130135驼鸟8000065表五关于哺乳动物的数据哺乳动物名称体重(千克)脉搏(次/分)哺乳动物名称体重(千克)脉搏(次/分)小蝙蝠0.006588海豹2025100小家鼠0.017500山羊3381仓鼠0.103347绵羊507080小猫0.117300猪1006080大家鼠0.252352马3804503455天竺鼠0.437269牛5004653兔1.34251象200030002550这里只对该问题作一

24、些拟合方面的练习。其它问题读者可自己进行讨论。符号用w表示动物的体重，单位：千克用v表示动物的基础血液流量，单位：公升/分用t表示动物的年龄，单位：岁用n表示动物的脉搏，单位：次/分假设动物的基础血液流量与动物的体重之间存在一定的函数关系V=f(w)，可以用表一中的数据来拟合这个函数。函数f(w)是一个什么样的函数呢？由于我们对“动物的基础血液流量与动物的体重”之间的关系并不清楚，所以只有根据表一中的数据得出函数f(w)一些性质。先将表一中的数据用MATLAB软件作出图形。从图上可以看出，这个函数关系V二f(w)应当是一个单调增加的函数。因此，拟合的函数如果不具有这一性质的话，就不能作为是好的

25、选择。一般地，可以假设函数f(w)是一个多项式，通常，这个多项式的次数不要超过3、4次，具体可根据拟合的效果来定。当然也可以用其它函数来拟合。为了提高拟合的效果，函数f(w)还可以用分段函数来拟合。以下是用分段函数拟合的结果：0.0033w4+0.1706w3-2.9727w2+21.3418w-43.0649,v=f(w)=we4.1,16-0.204756w2+9.31526w-74.6265,we16,24拟合函数图形是：问题1写出拟合函数v二f(w)和作出上面图形的MATLAB指令。同样可以拟合人的基础血液流量与体重之间的函数关系v=g(w)，可以用表二中的数据来拟合这个函数。这里用4

26、次多项式来拟合，拟合的结果是：v=0.000025w4+0.00351W30.1728w2+4.340w16.98,we18,72拟合函数图形是：问题2：写出拟合函数v=g(w)和作出上面图形的MATLAB指令。将上面拟合出来的函数v二f(w)和v二g(w)在它们的公共定义域18,24上的图形画出来，如下图所示。从图形上可以看出人类与动物之间的差异。人和动物在体重的公共区间上关于基础血液流量函数的图形sG=.-._tIk(xlr问题3：写出作出上面图形的MATLAB指令。下面考虑动物、人类的体重与基础脉搏的函数关系。假设人类的体重与基础脉搏之间的函数关系是w=H1(n)，利用表二中的数据来拟合

27、这个函数。这里用3次多项式拟合。拟合的结果是：w二H(n)二0.0002566n30.143347n2+16.2450n450.216,ne60,961其图形是：人的基础脉搏与体重之间的函数尸片恤;和原始数据图形1jUIIIIIII滞60-2040*riFC1*ifl/3-i-c2*w2+c3*w+c4+原始数据0IIIII=_11160657075808590951D0w问题4写出拟合函数w=H1(n)和作出上面图形的MATLAB指令。假设哺乳动物的体重与基础脉搏之间的函数关系是w=H(n)，利用表五中的数据来拟合这个函数。2这里用分段函数来拟合。由于当n100w=H2(n)彳7.56493

28、x1012“八2,n4000100.100.30300.150.45500.200.60建立模型的目的选择拖船的船型和船速，使冰山到达目的地后，可得到的每立方米淡水所花的费用最低，并与海水淡化的费用相比较。根据建模的目的和搜集到的有限的资料，需要对问题作如下的简化和假设：模型假设拖船航行过程中船速不变，航行不考虑天气等任何因素的影响，总航行距离为9600千米。山形状为球形，球面各点融化速率相同。冰山到达目的地后，每立方米冰可融化成0.85立方米水。模型构成首先需要知道冰山体积在运输过程中的变化情况;然后是计算航行过程中的燃料消耗；由此可以算出到达目的地后的冰山体积和运费。在计算过程中需要根据收

29、集到的数据拟合出经验公式。问题：1.理解建模过程中每一步的作用；用MATLAB软件求解这个模型，写出相应的MATLAB指令。五、水塔流量估计某居民区有一供居民用水的圆柱形水塔，一般可以通过测量其水位来估计水的流量。但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量。通常水泵每天供水一两次，每次约两个小时。水塔是一个高12.2米、直径17.4米的正圆柱。按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作。下表是某一天的水位测量记录、试估计任何时刻（包括水泵正供水时）从水塔

30、流出的水流量，及一天的总用水量。时刻（h）00.921.842.953.874.985.907.017.938.97水位（cm）968948931913898881869852839822时刻（h）9.9810.9210.9512.0312.9513.8814.9815.9016.8317.93水位（cm）/108210501021994965941918892时刻（h）19.0419.9620.8422.0122.9623.8824.9925.91水位（cm）866843822/105910351018六、加工工序的问题设有14件工件等待在一台机床上加工，某些工件的加工必须安排在另一些工件加

31、工完后才能进行第j件工件工件号j12345678910111213142028251642123210242040243616前期工件号i3，45，7，85，910，113，8，943，5，744，76，7，145，121，2，6（1）.若给出一个加工顺序，则确定了每个工件完工的时间（包括等待与加工两个阶段），试设计一个满足条件的加工顺序，使各个工件完工的时间之和最小。（2）.若第j号工件紧接着第i号工件完工后开工，机床需要准备时间是Ji+j9ij试设计一个满足条件的加工顺序，使机床花费的总时间最小。（3）.若工件完工时间超过一定时限u，则需支付一定的补尝费，其数额等于超过的时间与费用率卩.之i12/14 /14i /14积:工件号j1234567891011121314费用121015161011108541010812对于u=100,t