数值与数据拟合模型.doc

第二讲 插值与数据拟合模型函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟用插值还是拟合，有时容易确定，有时则并不明显。 在数学建模过程中，常常需要确定一个变量依存于另一个或更多的变量的关系，即函数。但实际上确定函数的形式（线性形式、乘法形式、幂指形式或其它形式）时往往没有先验的依据。只能在收集的实际数据的基础上对若干合乎理论的形式进行试验，从中选择一个最能拟合有关数据，即最有可能反映实际问题的函数形式，这就是数据拟合问题。一、插值方法简介插值问题的提法是，已知个节点，其中互不相同，不妨设，求任一插值点处的插值。可以看成是由某个函数产生的，的解析表达式可能十分复杂，或不存在封闭形式。也可以未知。求解的基本思路是，构造一个相对简单的函数，使通过全部节点，即，再由计算插值，即。1拉格朗日多项式插值插值多项式从理论和计算的角度看，多项式是最简单的函数，设是n次多项式，记作 （1）对于节点应有 （2）为了确定插值多项式中的系数，将（1）代入（2），有 （3）记方程组（3）简写成 （4）注意是Vandermonde行列式，利用行列式性质可得因互不相同，故，于是方程（4）中A有唯一解，即根据个节点可以确定唯一的n次插值多项式。拉格朗日插值多项式实际上比较方便的做法不是解方程（4）求A，而是先构造一组基函数： （5）是n次多项式，满足 （6）令 （7）显然是满足（2）的n次多项式，由方程（4）解的唯一性，（7）式表示的的解与（1）式相同。（5）、（7）称拉格朗日插值多项式，用计算插值称拉格朗日多项式插值。误差估计插值的误差通过插值多项式与产生节点的之差来估计，记作。虽然我们可能不知道的解析表达式，但不妨设充分光滑，具有阶导数。利用泰勒展开可以推出，对于任意。 （8）若可以估计 （9）则 （10）实际上因为常难以确定，所以（10）式并不能给出精确的误差估计。但是可能看出，n增加，减少；越光滑，越小，越小；x越接近，越小。例 将区间n等分，用产生个节点，然后作拉格朗日插值多项式。用计算（取4位有效数字）。估计（取）。解 若，则, 。由（5）、（7）式若，则，由（5）、（7）式。估计：对于可设，记节点间隔。当时于是（10）式给出可以算出n12340.30.04的精确值是0.8660（4位有效数字）的误差在范围内。插值多项式的振荡用拉格朗日插值多项式近似虽然随着节点个数的增加，的次数变大，多数情况下误差会变小，但n增加时，的光滑性变坏，有时会出现很大的振荡。理论上，当时，在内并不能保证处处收敛于。Runge给出了一个有名的例子：取。对于作，会得到如下图所示的结果。可以看出，对于较大的，随着n的增加，的振荡越来越大，事实上可以证明，仅当时，才有，而在此区间外，是发散的。高次插值多项式的这些缺陷，促使人们转而寻求简单的低次数多项式插值。2. 分段线性插值简单地说，将每两个相邻的节点用直线连起来，如此形成的一条折线就是分段线性插值函数，记作，它满足，且在每个小区间上是线性函数。可以表示为 （12） （13）有良好的收敛性，即对于有，。用计算点的插值时，只用到左右的两个节点，计算量与节点个数无关。但越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了，如数学、物理中用的特殊函数表，数理统计中用的概率分布表等。3. 三次样条插值样条函数的由来分段线性插值虽然简单，足够大时精度也相当高。但是折线在节点处显然不光滑，即在节点处导数不连续。这影响了它在诸如机械加工等领域(希望插值曲线光滑)中的应用。所谓样条(Spline)，来源于船舶、飞机等设计中描绘光滑外形曲线用的绘图工具。一根有弹性的细长木条用压铁固定在节点上，其它地方让它自然弯曲，如此画出的曲线称为样条曲线。因为这种曲线的曲率是处处连续的，所以要求样条函数的二阶导数连续。人们普遍使用的样条函数是分段三次多项式。三次样条函数三次样条函数 记作。要求它满足以下条件：a) 在每个小区间上是3次多项式；b) 在上二阶导数连续；c) 。 (14)由条件a，不妨将记为 （15）其中为待定系数，共4个。由条件b， （16）容易看出，(14)、(16)式共含有4-2个方程，为确定的4个待定参数，尚需再给出2个条件。最常用的是所谓自然边界条件： （17）可以证明，4阶线性方程组(14)、(16)、(17)有唯一解，即被唯一确定。但是，这种解法的工作量太大，方程组又常呈病态，所以实际上要设计简便的解法。另外，像分段线性函数一样，三次样条函数也有良好的收敛性，即在相当一般的条件下，。4. 用MATLAB作插值计算拉格朗日插值需先按照(5)、(7)式编写一个程序。设个节点以数组x0,y0输入(注意：程序中用个节点，而不是(5)、(7)式中的+1个节点)，m个插值点以数组x输入。输出数组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同长度)。例 对，用(=11)个节点(等分)作上述三种插值，用m(=21)个插值点(等分)作图，比较结果。插值方法小结拉格朗日插值是高次多项式插值(+1个节点上用不超过次的多项式)，插值曲线光滑，误差估计有表达式。但有振荡现象，收敛性不能保证。这种插值主要用于理论分析，实际意义不大。分段线性和三次样条插值是低次多项式插值，简单实用，收敛性有保证，但不光滑，三次样条插值的整体光滑性已大有提高，应用广泛，唯误差估计较困难。 二、最小二乘法简介下面先看一个例子。例1 “人口问题”是我国最大社会问题之一，估计人口数量和发展趋势是我们制定一系列相关政策的基础。有人口统计年鉴，可查的我国从1949年至1994年人口数据智料如下：年份1949195419591964196919741979198419891994人口数 (百万)541.67602.66672.09704.99806.71908.59975.421034.751106.761176.74分析：（1） 在直角坐标系上作出人口数的图象。（2） 估计出这图象近似地可看做一条直线。（3） 用以下几种方法（之一）确定直线方程，并算出1999年人口数。方法一：先选择能反映直线变化的两个点，如（1949，541.67），（1984，1034.75）二点确定一条直线，方程为：N = 14.088 t 26915.842，代入t =1999,得N 12.46亿。方法二：可以多取几组点对，确定几条直线方程，将t = 1999代入，分别求出人口数，再取其算数平值。方法三：可采用“最小二乘法”求出直线方程。最小二乘法简介设是直角平面坐标系下给出的一组数据，设，我们可以把这组数据看作是一个离散的函数。根据观察，如果这组数据图象“很象”一条直线（不是直线），我们的问题是确定一条直线，使得它能最好的反映出这组数据的变化。对个别观察值来说，用直线的值来近似代替其观察值时，所产生的误差可能是正的，也可能是负的。为了不使它们相加彼此抵消，可用来表示用直线来近似代原来实验数据时所产生的误差。为了在数学上处理方便，又把上式改成也就是说，我们选取常数，使得总误差达到最小。这就是所谓的最小二乘法。用微分法不难求出上面最小值问题的驻点，这里不列出其结果。事实上，在MATLAB中已有现成的求最小二乘问题的函数polyfit，称为多项式拟合函数，并且这个函数允许多项式的次数可以是任意次的。除外，还可以用解线性方程组中的除法运算（矩阵除法）来求解。这两个方法的区别在于：用polyfit函数求拟合问题时，多项式的次数必须从0次到最高次数之间每个次数都要出现。而如果需要选择一些次数进行拟合时，就可用矩阵除法运算来进行。矩阵除法还可以求一般的线性拟合问题，例如拟合函数不是多项式的线性拟合问题。上面例1中的问题就可以用polyfit来求解。例2 用最小二乘法求一个形如的经验公式，数据如下：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.2;x1=ones(5,1),x1;ab=x1yx0=19:0.2:44;y0=ab(1)+ab(2)*x0.2;clfplot(x,y,o, x0,y0,-r);多项式拟合是线性拟合问题（注意：无论拟合的多项式次是多少，多项式拟合都是线性拟合!）。但在实际应用中，有时还需要作非线性拟合问题。所谓线性拟合问题是指：需要拟合的函数中的未知常数都线性的。如函数中，常数是线性的。但、中的常数都是非线性。这种函数的拟合问题称为非线性拟合问题。有的非线性拟合问题可以化为线性拟合问题。例如在函数中，两边取对数，得，再令，则要拟合的函数就成，这样就变成线性拟合问题了。但也有不能化成线性拟合问题的情况，如函数就是这样。在MATLAB5.3中求非线性拟合问题的函数是curvefit。例3 在区间内拟合函数。解 用非线性拟合函数curvefit来拟合。先建立拟合函数。% 建立拟合函数，文件名是nxxyhhx.m，必须与函数名相同。% 要拟合的函数中参数用x表示，即x(1)=a x(2)=b；% 而拟合函数中x的值则用xdata表示。function v=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)hold ony=y1+1.2*(rand(size(x)-0.5); %将原型函数加一些扰动plot(x,y,*g)x0=2.5,-0.5; a=curvefit(nxxyhhx,x0,x,y) %用原始实验数据拟合函数nxxyhhx (x),vpa(a(1),a(2),8) % nxxyhhx (t)表达式中各项的系数。y2=nxxyhhx(a,x);plot(x,y2,-r) legend(原型函数,原始数据,用原始数据拟合的结果,4);三、血液流量问题小哺乳动物与小鸟的心跳速度比大哺乳动物与大鸟的快。如果动物的进化为每种动物确定了最佳心跳速度，为什么各种动物的最佳心跳速度不一样呢？由于热血动物的热量通过身体表面散失，所以它们要用大量的能量维持体温，而冷血动物在休息时只需要极少的能量，所以正在休息的热血动物似乎在维持体温。可以认为，热血动物可用的能量与通过肺部的血液流量成正比。（1）试建立一个模型，将体重与通过心脏的基础（即休息时的）血液流量联系起来，用下面的数据检验你的模型。（2）有许多可得到脉搏数据但没有血液流量数据的动物，建立一个模型将体重与基础脉搏联系起来，用下面的数据检验你的模型。（3）在检验你在（1）和（2）中的模型时会出现不一致，试进行分析。表一 关于某些哺乳动物的数据哺乳动物名称兔山羊狗狗狗体重（千克）4.12416126.4基础血液流量（分升/分）5.331221211表二 关于人类的数据年龄5101625334760体重（千克）18316668707270基础血液流量（分升/分）23335251434046脉搏（次/分）96906065687280表三 关于小鸟类的数据 表四 关于大鸟类的数据鸟类体重（克）脉搏（次/分）鸟类体重（克）脉搏（次/分）蜂鸟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 这里只对该问题作一些拟合方面的练习。其它问题读者可自己进行讨论。符号用表示动物的体重，单位：千克用表示动物的基础血液流量，单位：公升/分用表示动物的年龄，单位：岁用表示动物的脉搏，单位：次/分假设动物的基础血液流量与动物的体重之间存在一定的函数关系，可以用表一中的数据来拟合这个函数。函数是一个什么样的函数呢？由于我们对“动物的基础血液流量与动物的体重”之间的关系并不清楚，所以只有根据表一中的数据得出函数一些性质。先将表一中的数据用MATLAB软件作出图形。从图上可以看出，这个函数关系应当是一个单调增加的函数。因此，拟合的函数如果不具有这一性质的话，就不能作为是好的选择。一般地，可以假设函数是一个多项式，通常，这个多项式的次数不要超过3、4次，具体可根据拟合的效果来定。当然也可以用其它函数来拟合。为了提高拟合的效果，函数还可以用分段函数来拟合。以下是用分段函数拟合的结果：拟合函数图形是：问题1：写出拟合函数和作出上面图形的MATLAB指令。同样可以拟合人的基础血液流量与体重之间的函数关系，可以用表二中的数据来拟合这个函数。这里用4次多项式来拟合，拟合的结果是：拟合函数图形是：问题2：写出拟合函数和作出上面图形的MATLAB指令。将上面拟合出来的函数和在它们的公共定义域上的图形画出来，如下图所示。从图形上可以看出人类与动物之间的差异。问题3：写出作出上面图形的MATLAB指令。下面考虑动物、人类的体重与基础脉搏的函数关系。假设人类的体重与基础脉搏之间的函数关系是，利用表二中的数据来拟合这个函数。这里用3次多项式拟合。拟合的结果是：其图形是：问题4：写出拟合函数和作出上面图形的MATLAB指令。假设哺乳动物的体重与基础脉搏之间的函数关系是，利用表五中的数据来拟合这个函数。这里用分段函数来拟合。由于当时，变化激烈，所以用多项式已不能描述其变化的规律，可用其它函数来拟合。拟合的结果是： 图形如下。问题5：写出拟合函数和作出上面图形的MATLAB指令。假设小鸟类、大鸟类的体重与基础脉搏之间的函数关系分别是和，利用表三、四中的数据来拟合这两个函数。拟合的结果是：其图形如下。问题6：写出拟合函数和和作出上面图形的MATLAB指令。下面考虑人类的基础血液流量与基础脉搏的函数关系。假设人类的基础血液流量与基础脉搏之间的函数关系是，利用表二中的数据来拟合这个函数。拟合结果是：考虑复合函数=-.2489e-4*(.256557e-3*n3-.143347*n2+16.2450*n-450.216)4+.3506e-2*(.256557e-3*n3-.143347*n2+16.2450*n-450.216)3-.1728*(.256557e-3*n3-.143347*n2+16.2450*n-450.216)2+.1113457380e-2*n3-.622125980*n2+70.5033000*n-1970.917440;用MATLAB软件画出上面两个函数的图形：由上图可以看出，两者有较大的差异。原因就是前面的假设不合理，或不够完善。下面我们假设人类的的基础血液流量是体重和年龄的二元函数，用表二中的数据来拟合该函数。用二元二次多项式来拟合。结果是：问题7：写出以上相应的MATLAB程序和作出上面图形的MATLAB指令。四、冰山的运输在以盛产石油著称的波斯湾地区，浩瀚的沙漠覆盖着大地，淡水资源十分贫乏，不得不采用淡化海水的办法为国民提供用水，成本大约是每立方米淡水0.05英镑，有些专家提出从相距9600千米之遥的南极用拖船运送冰山到波斯湾，以取代淡化海水的办法。这个模型要从经济的角度研究冰山运输的可行性。为了计算用拖船运送冰山所获得的每立方米淡水所花的费用，我们需要关于拖船的租金、运量、燃料消耗及冰山运输过程中融化速度等方面的数据，以此作为建模所必须的准备工作。模型准备1.三种拖船的日租金和最大运量：船型小中大日租金(英镑)4.06.28.0最大运量(米3)51061071082.燃料消耗(英镑/千米)。主要依赖于船速和所运冰山的体积，船形的影响可以忽略：冰山的体积(米3)船速（千米/小时）10510610718.410.512.6310.813.516.2513.216.519.83.冰山在运输过程中的融化速率（米/天）。指在冰山与海水接触处每天融化的深度，融化速率除与船速有关外，还和运输过程中冰山到达处与南极的距离有关，这是由于冰山要从南极运往赤道附近的缘故。与南极的距离（千米）船速（千米/小时）0100040001350000.100.150.200.300.450.60建立模型的目的选择拖船的船型和船速，使冰山到达目的地后，可得到的每立方米淡水所花的费用最低，并与海水淡化的费用相比较。根据建模的目的和搜集到的有限的资料，需要对问题作如下的简化和假设：模型假设1. 拖船航行过程中船速不变，航行不考虑天气等任何因素的影响，总航行距离为9600千米。2. 山形状为球形，球面各点融化速率相同。3. 冰山到达目的地后，每立方米冰可融化成0.85立方米水。 模型构成首先需要知道冰山体积在运输过程中的变化情况;然后是计算航行过程中的燃料消耗；由此可以算出到达目的地后的冰山体积和运费。在计算过程中需要根据收集到的数据拟合出经验公式。问题：1.理解建模过程中每一步的作用；2.用MATLAB软件求解这个模型，写出相应的MATLAB指令。五、水塔流量估计某居民区有一供居民用水的圆柱形水塔，一般可以通过测量其水位来估计水的流量。但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量。通常水泵每天供水一两次，每次约两个小时。水塔是一个高12.2米、直径17.4米的正圆柱。按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作。下表是某一天的水位测量记录、试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量。时刻（h）水位（cm）0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 968 948 931 913 898 881 869 852 839 822时刻（h）水位（cm）9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93/ / 1082 1050 1021 994 965 941 918 892时刻（h）水位（cm）19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91866 843 822 / / 1059 1035 1018六、加工工序的问题设有14件工件等待在一台机床上加工，某些工件的加工必须安排在另一些工件加工完后才能进行，第j 件工件加工时间及先期必须完工的工件号由下表给出:工件号j12345678910111213142028251642123210242040243616前期工件号i3，45，7，85，9-10，113，8，943，5，74-4，76，7，145，121，2，6(1).若给出一个加工顺序，则确定了每个工件完工的时间(包括等待与加工两个阶段)，试设计一个满足条件的加工顺序，使各个工件完工的时间之和最小。(2).若第j号工件紧接着第i号工件完工后开工，机床需要准备时间是试设计一个满足条件的加工顺序，使机床花费的总时间最小。(3).若工件完工时间超过一定时限u，则需支付一定的补尝费，其数额等于超过的时间与费用率之积:工件号j1234567891011121314费用121015161011108541010812对于u = 100，安排一个满足条件的加工顺序，使总的补尝费用最少。上面的三个问题，它们都有一定的实际意义。对于问题(1)可以看着是:如果在某个时刻同时有很多客户需要加工自己的工件，而你用来加工这些工件的机床只有一台，这里假设机床不能同时加工两个及其以上的工件，这样客户们就需要排队等待加工，这时对于客户们来说当然希望等待的时间越少越好。对于问题(2)则可以看着是