六插值与拟合

(1)(3)第六章插值与拟合在实际中，常常需要确定一个变量依存于另一个或更多的变量的关系，即函数。但实际上确定函数的形式(线性形式、乘法形式、幕指形式或其它形式)时往往没有先验的依据。 只能在实验或测量得到的离散数据的基础上进行试验，这时候我们常常会用到插值与拟合的方法。插值与拟合就是要通过已知的数据去确定某一类已知函数的参数或寻找某个近似函数， 使所得到的近似函数对已知数据有较高的拟合进度。如果要求这个近似函数经过所有已知的数据点，则称此类问题为插值问题。当所给的数据较多时， 用插值方法所得到的插值函数会很复杂，所以，通常插值方法用于数据较少的情况。其实，通常情况下数据都是由观测或试 验得到的，往往会带有一定的随机误差，因而，要求近似函数通过所有的数据点也是没有必要的。如果不要求近似函数通过所有的数据点，而是要求他能较好地反映数据的整体变化趋势，则解决这类问题的方法称为数据拟合。虽然插值与拟合都是要构造已有数据的近似函数，但因对近似要求的准则不同，因此二者在数学方法上有很大的差异。第一节一般插值方法1. 问题的提出插值问题的一般提法：已知n -1个节点(xj, yj) (j二0,1/ n)，其中xj (j = 1,2 - , n)互不相同，要求构造一个函数y = f(x)使得yj = f (召)(j =0,1,n)。我们通常称这样一类问题为插值问题，并称构造的函数y = f (x)为插值函数，xj (j = 1 , 12, n为插值节点，X =f(Xj) (j = o ,IHn 为插值条件。2. 多项式插值从理论和计算的角度看，多项式是最简单的函数，设f (X)是n次多项式，记作Ln(x) =anXn an2 d x a对于节点(Xj,yj应有Ln (Xj ) = yj , j - 0,1,2，,n为了确定插值多项式Ln(x)中的系数a.,a.亠，c ,a，将(1 )代入(2), 有an Xo a“ 4X0 a1Xo a - yon9nXnnan4XnnXonX1n 4Xon 4X1=(an,an = ,ao)T,Y 二(yo, y1,yn)Tn 4Xn1 / 14anX； - an4X；a a。= %XA=Y注意：detX是Vandermonde行列式，利用行列式性质可得detX ：| 丨区-Xj)o 勺:k m因xj互不相同，故detX = 0，于是方程(4 )中A有唯一解，即根据n 1个节点可以 确定唯一的n次插值多项式。3.拉格朗日(Lagrange)插值多项式实际上比较方便的做法不是解方程(4)求A,而是先构造一组基函数:li(x)二(X Xo)(X Xid)(XXi 1)(X Xn)(Xi Xo) (Xi Xy)(X X 1) (Xi -Xn),i =0,1,2, ,n(5)li(x)是n次多项式，满足h(Xj) J： = i,j =0,12，n(6)0, i和nLn(x) yli(X)(7)i =0显然Ln(x)是满足(2)的n次多项式，由方程(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 (x0) f (xjf(X,X1)01(二个节点，一阶差商)X0 Xjf(X0,XJ - f(XX2)f(X,X1,X2)-(二个节点，二阶差商)X0 X2f (X0,Xt，,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上是线性函数 (j 71, , n)。I n (x)可以表示为(12)ln(X)yjlj(X)j =0gXjxj _ Xj 二(j 0时舍去)x x -(13)I 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中分段线性插值有现成的程序：y=interp1(x0， y0，x)其中输入x0、y0、x分别表示节点数据(xj,yj) (j = 0,1/ n)和插值点以数组，输出y为插值点对应的插值。数组长度自定义(x0和y0同长度，x和y同长度)。第二节样条函数插值方法1. 样条函数的由来分段线性插值虽然简单， n足够大时精度也相当高。 但是折线在节点处显然不光滑， 即 In(x)在节点处导数不连续。这影响了它在诸如机械加工等领域(希望插值曲线光滑)中的应用。所谓样条(Spline)，来源于船舶、飞机等设计中描绘光滑外形曲线用的绘图工具。一根 有弹性的细长木条用压铁固定在节点上，其它地方让它自然弯曲，如此画出的曲线称为样条曲线。因为这种曲线的曲率是处处连续的，所以要求样条函数的二阶导数连续。人们普遍使用的样条函数是分段三次多项式。2. 三次样条函数三次样条函数 记作S(x),a空x乞b。要求它满足以下条件：a) 在每个小区间Xi*。=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)=Si；1(Xi)i =1,2，n-1( 16)0(xi) = SE(Xi)容易看出，(14)、(16)式共含有4n-2个方程，为确定 S(x)的4n个待定参数，尚需再 给出2个条件。最常用的是所谓自然边界条件：3 / 14S”(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、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次的多项式)，插值曲线光滑，误差估计有表达式。但有振荡现象，收敛性不能保证。这种插值主要用于理论分析，实 际意义不大。分段线性和三次样条插值是低次多项式插值，简单实用，收敛性有保证，但不光滑，三次样条插值的整体光滑性已大有提高，应用广泛，唯误差估计较困难。第三节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)ans =-6.37502. Runge现象及分段线性插值2.1 Runge 现象Runge在本世纪初发现：在-1,1上用n 1个等距结点作插值多项式Pn(x)，使其在各结点的值与函数 y(x) =1/(1 25x2)在结点的值相等。但在n 时，插值多项式R(x)在区间中部趋于y(x)。但对于0.726 1 x I iiiiiqy=0 21 62 70 77 110;100coef=polyfit(x,y,1);80*a1=coef(1),a0=coef(2);t-4-ybest=a1*x+aO;60- ” s=sum(y-ybest).A2);40k-axis(-1,6,-20,120);20l+-pIot(x,y, *)hoId on0iiIiiijjjC0.511.522.533.544.55pIot(x,ybest)例2如下给出从二阶到十阶多项式拟合曲线的比较程序，并给出拟合曲线x=0 1 2 3 4 5;y=0 21 62 70 77 110;xi=0:0.2:5;for n=2:10bb=polyfit(x, y,n);yi=polyval(bb,xi); plot(xi,yi,x,y, * )title(int2str(n),次多项式拟合曲线) grid on pause100806040200r-204次多项式拟合曲线12000.511.522.533.544.55end2例3用最小二乘法求一个形如 y =a bx的经验公式，数据如下:x1925313844y19.032.349.073.398.8解用求矩阵除法(因为要拟合的多项式缺了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, -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=y1+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.血液流量问题小哺乳动物与小鸟的心跳速度比大哺乳动物与大鸟的快。如果动物的进化为每种动物确定了最佳心跳速度，为什么各种动物的最佳心跳速度不一样呢？由于热血动物的热量通过身 体表面散失，所以它们要用大量的能量维持体温，而冷血动物在休息时只需要极少的能量， 所以正在休息的热血动物似乎在维持体温。可以认为，热血动物可用的能量与通过肺部的血液流量成正比。（1）试建立一个模型，将体重与通过心脏的基础（即休息时的）血液流量联系起来， 用下面的数据检验你的模型。（2）有许多可得到脉搏数据但没有血液流量数据的动物，建立一个模型将体重与基础 脉搏联系起来，用下面的数据检验你的模型。（3）在检验你在（1 ）和（2）中的模型时会出现不一致，试进行分析。鸟类体重（克）脉搏 （次 /分）蜂鸟4:615鹪鹩11450金丝雀16514麻雀28350鸽子130135表五关于哺乳动物的数据鸟类体重（克）脉搏 （次 /分）海鸥388401鸡1980312秃鹰8310P 199火鸡875093驼鸟8000065表一关于某些哺乳动物的数据哺乳动物名称兔山羊狗狗狗体重（千克）4.12416126.4基础血液流量（分升/分）5.331十221211表二关于人类的数据年龄5101625334760体重（千克）18316668707270基础血液流量（分升/分）23335251434046脉搏（次/分）96906065687280表三 关于小鸟类的数据表四 关于大鸟类的数据哺乳动物名 称体重（千克）脉搏（次/分）哺乳动物名 称体重（千克）脉搏（次/分）小蝙蝠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 = f(W)，可以用表一中的数据来拟合这个函数。函数f (w)是一个什么样的函数呢？由于我们对“动物的基础血液流量与动物的体重” 之间的关系并不清楚，所以只有根据表一中的数据得出函数f (w) 一些性质。先将表一中的数据用 MATLAB 软件作出图形。从图上可以看出，这个函数关系 V = f(w)应当是一个单调增加的函数。因此，拟合的函数如果不具有这一性质的话，就不 能作为是好的选择。一般地，可以假设函数f(w)是一个多项式，通常，这个多项式的次数不要超过3、4次，具体可根据拟合的效果来定。当然也可以用其它函数来拟合。为了提高拟合的效果，函数f(w)还可以用分段函数来拟合。以下是用分段函数拟合的结果：-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,w 18,72拟合函数图形是：人的基砒血液流量与体重之间的函数丫利帥和原始数据图形问题2:写出拟合函数 V二g(w)和作出上面图形的 MATLAB指令。将上面拟合出来的函数 v = f(w)和v =g(w)在它们的公共定义域18,24上的图形画 出来，如下图所示。从图形上可以看出人类与动物之间的差异。人和动物在体重的公共区间上黄于基础血液流量函数的图形问题3:写出作出上面图形的 MATLAB指令。 下面考虑动物、人类的体重与基础脉搏的函数关系。假设人类的体重与基础脉搏之间的函数关系是w = Hi(n)，利用表二中的数据来拟合这个函数。这里用 3次多项式拟合。拟合的结果是：wnHjn) =0.0002566? -0.143347n216.2450n-450.216, n 60,96其图形是：人的基础脉撐与体重之间的函数庐片伸闲I原始数据图形Q I1IIII60 -40+ n=c+c+原始数据 1160657075 冊 669095100问题4:写出拟合函数 wHn)和作出上面图形的 MATLAB指令。假设哺乳动物的体重与基础脉搏之间的函数关系是w = H2(n)，利用表五中的数据来拟合这个函数。这里用分段函数来拟合。由于当n : 100时，w变化激烈，所以用多项式已不能描述其变化的规律，可用其它函数来拟合。拟合的结果是:-0.0345386n 14.7835, n _ 100w = H2 (n) = 7.56493汉 1012 “6,n 4000100.100.30300.150.45500.200.60建立模型的目的选择拖船的船型和船速，使冰山到达目的地后，可得到的每立方米淡水所花的费用最低， 并与海水淡化的费用相比较。根据建模的目的和搜集到的有限的资料，需要对问题作如下的简化和假设：模型假设1. 拖船航行过程中船速不变，航行不考虑天气等任何因素的影响，总航行距离为9600千米。2. 山形状为球形，球面各点融化速率相同。3. 冰山到达目的地后，每立方米冰可融化成0.85立方米水。模型构成首先需要知道冰山体积在运输过程中的变化情况；然后是计算航行过程中的燃料消耗；由此可以算出到达目的地后的冰山体积和运费。在计算过程中需要根据收集到的数据拟合出经验公式。 问题：1.理解建模过程中每一步的作用；2. 用MATLAB件求解这个模型，写出相应的MATLAB旨令。