微积分模型ppt

1、模型分析和建立模型分析和建立先看如下三个例子先看如下三个例子1. 卫星轨道的长度2. 射击命中概率3. 人口增长模型1、卫星轨道的长度 人造地球卫星的轨道可以视为平面上的椭圆，中国第一颗人造卫星近地点距离地球表面439km，远地点距离地球表面2384km，地球半径为6371m，试求该卫星的轨道长度？分析：2、射击命中概率 炮弹射击的目标为一正椭圆形区域，当瞄准目标的中心发射时，在纵多因素的影响下，弹着点与目标中心有随机偏差。可以合理地假设弹着点围绕中心呈二维正态分布，且偏差在x方向和y方向相互独立。若椭圆区域在x方向半轴长120m,y方向半轴80m，设弹着点偏差的均方差在x方向和y方向均为10

2、0m,试求 炮 弹 落 在 椭 圆 形 区 域 内 的 概 率 ？分析：3、人口增长率20世纪美国人口统计数据如表6.5所示，计算表中这些年份的人口增长率。又已知某地区20世纪70年代的人口增长率如表6.6所示，且1970年人口为210万，试估计该地区1980年的人口？分析:小小结结模型求解u 数值微分相关知识回顾u 数值积分相关知识u Matlab 相关命令u 程序实现u 数值微分相关知识回顾hhxfhxfhhxfxfhxfhxfxfhhh2)()(lim)()(lim)()(lim)( 0001. 函数f(x)以离散点列给出时，而要求我们给出导数值，2. 函数f(x)过于复杂这两种情况都要

3、求我们用数值的方法求函数的导数值微积分中，关于导数的定义如下：自然，而又简单的方法就是，取极限的近似值，即差商 向前差商公式()( )( ),0,1iiif xhf xfxih 向后差商公式 中心差商公式( )()( ),1,2iiif xf xhfxih111()()()2f xhf xhfxh给定点列20)(,(iiixfx且2110,xxxxh)( ),( ),( 012xfxfxf求两点差商公式：两点差商公式：三点差商公式三点差商公式0000012201()3 ()4 ()(2 )21 3 ()4 ()()fxf xf xhf xhhf xf xf xxx1002201()()(2 )

4、21 ()()fxf xf xhhf xf xxx2000012201()()4 ()3 (2 )21 =()4 ()3 ()fxf xf xhf xhhf xf xf xxxu 数值积分相关知识回顾)()()(aFbFdxxfba关于积分，有Newton-Leibniz公式但是，在很多情况下，还是要数值积分：1、函数有离散数据组成2、F(x)求不出3、F(x)非常复杂定义数值积分如下：是离散点上的函数值的线性组合)()(0iniinxfafI称为求积系数，与f(x)无关，与积分区间和积分点有关矩形公式矩形公式左矩形公式左矩形公式( )() ( )baf x dxba f a右矩形公式右矩形公

5、式( )() ( )baf x dxb a f b中矩形公式中矩形公式( )() ()2babaf x dxba f梯形公式梯形公式( )( )( )2babafx dxf af b辛普森辛普森(Simpson)公式公式( ) ( )4 ()( )62babaabf x dxf aff b复合梯形公式复合梯形公式11011101( )( ) ( )2()( )22kknbxaxknnkkkkkfx dxfx dxhhfxfxf afxf b复合辛普森复合辛普森(Simpson)公式公式111012112( ) ( )4()2()( )6=2nnbkakkkkkkbaf x dxf af xf

6、xf bxxx其中 将将 分为分为 等份，步长等份，步长 ，分点，分点, a bbahnn,0,1,kxakh knMatlab相关命令 x=linspace(0,pi,50);y=sin(x);R1=trapz(x,y)；% 梯形公式计算R2=quad(sin,0,pi,1.0e06)；% 辛普森公式计算结果R1=1.999314849324062R2=1.9999934965349640sin( )x dx程序实现模型实例讲解 森林救火问题 水箱的水流问题 森林救火问题森林救火问题森林失火后，要确定派出消防队员的数量。森林失火后，要确定派出消防队员的数量。队员多，森林损失小，救援费用大；队

7、员多，森林损失小，救援费用大；队员少，森林损失大，救援费用小。队员少，森林损失大，救援费用小。所以需要综合考虑森林损失费和救援费与消防队员所以需要综合考虑森林损失费和救援费与消防队员之间的关系，以总费用最小来决定派出队员的数量。之间的关系，以总费用最小来决定派出队员的数量。问题问题分析分析问题问题记队员人数记队员人数x, 失火时刻失火时刻t=0, 开始救火时刻开始救火时刻t1, 灭火时刻灭火时刻t2, 时刻时刻t森林烧毁面积森林烧毁面积B(t). 损失费损失费f1(x)是是x的减函数的减函数, 由烧毁面积由烧毁面积B(t2)决定决定. 救援费救援费f2(x)是是x的增函数的增函数, 由队员人数

8、和救火时间决定由队员人数和救火时间决定.找到找到恰当的恰当的x，使，使f1(x), f2(x)之和最小之和最小 关键是对关键是对B(t)作出合理的简化假设作出合理的简化假设.问题问题分析分析失火时刻失火时刻t=0, 开始救火时刻开始救火时刻t1, 灭火时刻灭火时刻t2, 画出时刻画出时刻 t 森林烧毁面积森林烧毁面积B(t)的大致图形的大致图形t1t20tBB(t2)分析分析B(t)比较困难比较困难,转而讨论森林烧毁转而讨论森林烧毁速度速度dB/dt. （其中（其中dB/dt表示表示单位时间烧毁的面单位时间烧毁的面积）积）模型假设模型假设 3）f1(x)与与B(t2)成正比，系数成正比，系数c

9、1 (烧毁单位面积损失费）烧毁单位面积损失费） 1）0 t t1, dB/dt 与与 t成正比，系数成正比，系数 (火势蔓延速度）火势蔓延速度） 2）t1 t t2, 降为降为 - x ( 为队员的平均灭火为队员的平均灭火速度）速度） 4）每个）每个队员的单位时间灭火费用队员的单位时间灭火费用c2, 一次性费用一次性费用c3 (运运送消防队员和器材等一次性支出送消防队员和器材等一次性支出)假设假设1）的解释的解释 rB火势以失火点为中心，火势以失火点为中心，均匀向四周呈圆形蔓延，均匀向四周呈圆形蔓延，半径半径 r与与 t 成正比成正比面积面积 B与与 t2成正比，成正比， dB/dt与与 t成

10、正比成正比.xbtt12202)()(tdttBtB模型建立模型建立dtdBb0t1tt2x假设假设1）,1tbxcttxcxftBcxf31222211)()(),()(目标函数目标函数总费用总费用)()()(21xfxfxC假设假设3）4）xttt112假设假设2）)(222212212xttbt0dxdCxcxxtcxtctcxC3122121211)(22)(模型建立模型建立目标函数目标函数总费用总费用模型求解模型求解求求 x使使 C(x)最小最小231221122ctctcx结果解释结果解释 / 是火势不继续蔓延的最少队员数是火势不继续蔓延的最少队员数dtdBb0t1t2tx其中其中

11、 c1,c2,c3, t1, , 为已知参数为已知参数模型模型应用应用c1,c2,c3已知已知, t1可估计可估计, c2 x c1, t1, x c3 , x 结果结果解释解释231221122ctctcxc1烧毁单位面积损失费烧毁单位面积损失费, c2每个每个队员单位时间灭火费队员单位时间灭火费, c3每个每个队员一次性费用队员一次性费用, t1开始救火时刻开始救火时刻, 火火势蔓延速度势蔓延速度, 每个每个队员平均灭火队员平均灭火速度速度.为什么为什么? ? , 可可设置一系列数值设置一系列数值由模型决定队员数量由模型决定队员数量x练习题 在森林救火模型中，如果考虑消防队员的灭火速度 与

12、开始救火时的火势b有关，试假设一个合理的函数关系，重新求解模型。程序实现 美国某州的用水管理机构要求各社区提供以每小时多少加美国某州的用水管理机构要求各社区提供以每小时多少加仑计的用水量以及每天所用的总水量。许多社区没有测量流入仑计的用水量以及每天所用的总水量。许多社区没有测量流入或流出水塔水量的装置，只能代之以每小时测量水塔中的水位，或流出水塔水量的装置，只能代之以每小时测量水塔中的水位，其误差不超过其误差不超过5%。需要注意的是，当水塔中的水位下降到最低。需要注意的是，当水塔中的水位下降到最低水位水位L时，水泵就自动向水塔输水直到最高水位时，水泵就自动向水塔输水直到最高水位H，此期间不能，

13、此期间不能测量水泵的供水量，因此，当水泵正在输水时不容易建立水塔测量水泵的供水量，因此，当水泵正在输水时不容易建立水塔中水位和用水量之间的关系。水泵每天输水一次或两次，每次中水位和用水量之间的关系。水泵每天输水一次或两次，每次约约2小时。小时。水箱的水流问题水箱的水流问题 试估计任何时刻（包括水泵正在输水时间）从水塔试估计任何时刻（包括水泵正在输水时间）从水塔流出的水流量流出的水流量f(t)，并估计一天的总用水量。已知该水塔，并估计一天的总用水量。已知该水塔是一个高为是一个高为40ft（英尺），直径为（英尺），直径为57ft（英尺）的正圆（英尺）的正圆柱，表柱，表6-1给出了某个小镇一天水塔水

14、位的真实数据，给出了某个小镇一天水塔水位的真实数据，水位降至约水位降至约27.00ft水泵开始工作，水位升到水泵开始工作，水位升到35.50ft停止停止工作。（注：工作。（注：1ft（英尺）（英尺）=0.3024m（米）（米）表表6.1 某小镇某天水塔水位某小镇某天水塔水位时间时间/s水位水位/0.01ft时间时间/s水位水位/0.01ft031754663633503316311049953326066353054539363167106192994572543087139372947605743012179212892645542927212402850685352842252232795

15、71854276728543275275021269732284269779254水泵开动水泵开动35932水泵开动水泵开动82649水泵开动水泵开动39332水泵开动水泵开动8596834753943535508995333974331834459327033402.2 问题分析问题分析 模型中所指流量可视为单位时间内流出水的模型中所指流量可视为单位时间内流出水的体积体积。由于水塔是正圆柱形，横截面积是常数，所以在水泵由于水塔是正圆柱形，横截面积是常数，所以在水泵不工作时段，流量很容易根据水位相对时间的变化算不工作时段，流量很容易根据水位相对时间的变化算出。问题的难点在于如何估计水泵供水时段

16、的流量。出。问题的难点在于如何估计水泵供水时段的流量。 水泵供水时段的流量只能靠供水时段前后的流量经水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到。作为用于插值或拟合的原始数据，我插值或拟合得到。作为用于插值或拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。这此流量大们希望水泵不工作时段的流量越准确越好。这此流量大体上可由两种方法计算，一是直接对表体上可由两种方法计算，一是直接对表6-1中的水量用中的水量用数值微分算出各时段的流量，用它们拟合其它时刻或连数值微分算出各时段的流量，用它们拟合其它时刻或连续时间的流量；二是先用表中数据拟合水位一时间函数，续时间的流量；二是先用表

17、中数据拟合水位一时间函数，求导数即可得到连续时间的流量。求导数即可得到连续时间的流量。 有了任何时刻的流量，就不难计算一天的总用水有了任何时刻的流量，就不难计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录量。其实，水泵不工作时段的用水量可以由测量记录直接得到，由表直接得到，由表6-1中下降水位乘以水搭的截面积就是中下降水位乘以水搭的截面积就是这一时段的用水量。这个数值可以用来检验数据插值这一时段的用水量。这个数值可以用来检验数据插值或拟合的结果。或拟合的结果。 在具体给出本问题的解答之前，先介绍一个简单在具体给出本问题的解答之前，先介绍一个简单的数据插值方法。的数据插值方法。2.

18、3 拉格朗日插值拉格朗日插值1、线性插值、线性插值)(xf假设已知假设已知在区间在区间,ba上的两个结点上的两个结点10, xx和它们的函数值和它们的函数值 10 xxx )()()(10 xfxfxf求一个一次多项式求一个一次多项式xaaxP101)(，使得多项式，使得多项式)(1xP在结点上满足条件在结点上满足条件 1, 0),()(1ixfxPii 这种插值方法称为线性插值方法（也称两点插值）。这种插值方法称为线性插值方法（也称两点插值）。 可以求出：可以求出： )()()(101001011xfxxxxxfxxxxxP2、抛物插值、抛物插值)(xf已知已知在区间在区间,ba上的三个结点

19、上的三个结点210,xxx和它们的函数值和它们的函数值 210 xxxx )()()()(210 xfxfxfxf求一个次数不超过求一个次数不超过2的多项式的多项式22102)(xaxaaxP，使得它在结点上满足条件，使得它在结点上满足条件2, 1, 0),()(2ixfxPii这种插值方法称为抛物线插值法，这种插值方法称为抛物线插值法， 可求出：可求出：)()( )()( )()()( )()( )()(121012002010212xfxxxxxxxxxfxxxxxxxxxP)()( )()( )(2021210 xfxxxxxxxx3、n次拉格朗日插值次拉格朗日插值,ba假设取区间假设取

20、区间 上的上的1n个结点个结点nxxx,10，并且已知函数，并且已知函数)(xf在这此点的函数值在这此点的函数值 nxxxxx210 )()()()()(210nxfxfxfxfxf现在求一个次数不超过现在求一个次数不超过n的多项式的多项式nnnxaxaxaaxP2210)(，使得满足条件，使得满足条件 nixfxPiin, 1, 0),()( 这种插值方法称为这种插值方法称为n次多项式插值（或称代数插值），次多项式插值（或称代数插值）， 利用拉格朗日插值插值方法可得利用拉格朗日插值插值方法可得)()()()()()()()()(,2, 21, 10, 0nnnnnnnxfxlxfxlxfxl

21、xfxlxP0( )0,1,njknjkjj kxxlxknxx其中上述多项式称为上述多项式称为 n次拉格朗日（次拉格朗日（Lagrange）插值多项式，）插值多项式， 函数函数nkxlnk, 1, 0),(,称为拉格朗日插值基函数。称为拉格朗日插值基函数。当当n=1,2时，时，n次拉格朗日（次拉格朗日（Lagrange）插值多项式即）插值多项式即为线性插值多项式和抛物插值多项式。为线性插值多项式和抛物插值多项式。4321x例例6.1 已知函数发已知函数发f(x)的函数表如下：的函数表如下： 31 . 28 . 15 . 2)(xf求其拉格朗日插值多项式，并求求其拉格朗日插值多项式，并求)5

22、. 2(f的近似值。的近似值。 解解 由于给出了由于给出了4个插值结点，所以可做出次数不超个插值结点，所以可做出次数不超过过3的拉格朗日插值多项式。的拉格朗日插值多项式。43132361)41 ( )31 ( )21 ()4( )3( )2()(233 , 0 xxxxxxxl6219421)42( )32( )22()4( )3( )2()(233 , 1xxxxxxxl472721)43( )33( )23()4( )3( )2()(233 , 2xxxxxxxl161161)44( )34( )24()4( )3( )2()(233 , 3xxxxxxxl 将上列将上列4式代入式代入n=

23、3 的拉格朗日插值公式，可得所的拉格朗日插值公式，可得所要求的插值多项式为要求的插值多项式为6 . 4933. 29 . 00667. 0)(233xxxxP将将x=2.5代入可得代入可得f(2.5)的近似值为的近似值为1.8496。 拉格朗日插值法适合节点较少的情况，当节点较多拉格朗日插值法适合节点较少的情况，当节点较多的大范围高次插值的逼近效果往往并不理想且当插值的大范围高次插值的逼近效果往往并不理想且当插值结点增加时，计算越来越繁。为了提高精度和减少计结点增加时，计算越来越繁。为了提高精度和减少计算，还有牛顿插值法下、三次样条插值等，具体可参算，还有牛顿插值法下、三次样条插值等，具体可参

24、阅有关书籍。阅有关书籍。1一维插值命令一维插值命令interp1的具体使用格式的具体使用格式yy=interp1(x,y,xx,method) 其中其中x,y是插值结点的数据向量，如果是插值结点的数据向量，如果y是矩阵，是矩阵，则对矩阵则对矩阵y的每一列相对的每一列相对x进行插值，进行插值，xx是待求函数值是待求函数值的插值结点向量，可以缺省。的插值结点向量，可以缺省。method是可选项，是可选项，说明插值使用的方法。对于一维插值，说明插值使用的方法。对于一维插值，MATLAB提供提供可选的方法有：可选的方法有：nearest，linear，spline，cubic，它们分别表示最近插值、线

25、性插值、三次样条插值和它们分别表示最近插值、线性插值、三次样条插值和三次插值。三次插值。2.4 MATLAB软件实现插值法软件实现插值法 MATLAB软件提供了专门做各种插值的命令：软件提供了专门做各种插值的命令：interp1（一维插值），（一维插值），interp2（二维插值），（二维插值），interp3（三维插值），（三维插值），interpn（n维插值），维插值），spline（样条插值）等。（样条插值）等。 2二维插值命令二维插值命令interp2的具体使用格式的具体使用格式zz=interp2(x,y,z,xx,yy,method) 该指令的意思是根据数据向量该指令的意思是根据数

26、据向量x，y，z按按method指定的方法来做插值，然后将指定的方法来做插值，然后将xx，yy处插值处插值函数的插值结点向量，如果函数的插值结点向量，如果xx，yy在插值范围之内，在插值范围之内，则返回值在则返回值在zz中，否则返回值为空中，否则返回值为空NaN。method是插值方法可选项，具体要求同一是插值方法可选项，具体要求同一维插值的情况。维插值的情况。 该命令还有以下几种省略格式：该命令还有以下几种省略格式：zz=interp2(z,xx,yy)zz=interp2(x,y,z,xx,yy)zz=interp2(z,ntimes)3三维插值命令三维插值命令interp3的具体使用格式

27、的具体使用格式vi=interp3(x,y,z,v,xi,yi,zi,method) 它的具体含义跟前面的一、二维插值是相似的，在它的具体含义跟前面的一、二维插值是相似的，在此不作解释，读者可在此不作解释，读者可在MATLAB工作空间中用工作空间中用help interp3命令获得。命令获得。4样条插值命令样条插值命令spline的具体使用格式的具体使用格式yy=spline(x,y,xx)它的意思等同于命令它的意思等同于命令yy=interp1(x,y,xx,spline)例例6.2 在用外接电源给电容器充电时，电容器两端在用外接电源给电容器充电时，电容器两端的电压的电压V将会随着充电时间将

28、会随着充电时间t发生变化，已知在某一次发生变化，已知在某一次实验时，通过测量得到下列观测值，分别用拉格朗日实验时，通过测量得到下列观测值，分别用拉格朗日插值法、分段线性插值法、三次样条插值法画出插值法、分段线性插值法、三次样条插值法画出V随随着时间着时间t变化的曲线图，分别计算当时间变化的曲线图，分别计算当时间t=7s时，三时，三种插值法各自算得电容器两端电压的近似值。种插值法各自算得电容器两端电压的近似值。1295 . 64321st 4 .101 .106 . 90 . 92 . 83 . 72 . 6vV解解:由于由于MATLAB没有提供现成的拉格朗日插值命令，没有提供现成的拉格朗日插值

29、命令，我们可以编写一个函数我们可以编写一个函数lglrcz.m来完成，其他两种插值来完成，其他两种插值法可用现成的命令。法可用现成的命令。 用用MATLAB软件进行三种插值计算的程序为软件进行三种插值计算的程序为Three_interpmethods.m。程序程序lagrange.m：function y=lagrange(x0,y0,x) n=length(x0)； m=length(x) ； for i=1:m z=x(i) ； s=0.0； for k=1:n p=1.0； for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j) ； end end s=p*

30、y0(k)+s； end y(i)=s； end程序为程序为Three_interpmethods.mt=1,2,3,4,6.5,9,12 ；v=6.2,7.3,8.2,9.0,9.6,10.1,10.4 ；t0=0.2:0.1:12.5；lglr=lagrange(t,v,t0) ；laglr=lagrange(t,v,7) ；fdxx=interp1(t,v,t0) ；fendxx=interp1(t,v,7) ；scyt=interp1(t,v,t0,spline) ；sancyt=interp1(t,v,7,spline)plot(t,v,*,t0,lglr,r,t0,fdxx,g,t

31、0,scyt,b)gtext(lglr)gtext(fdxx)gtext(scyt)执行结果是执行结果是laglr= 9.52988980716254fendxx= 9.70000000000000sancyt= 9.67118039327444图形如图图形如图6.1所示。所示。图中曲线图中曲线lglr表示拉格朗日插值曲线，表示拉格朗日插值曲线，scyt表示三次样表示三次样条插值曲线，条插值曲线，fdxx表示分段线性插值曲线。表示分段线性插值曲线。 2.5 问题求解问题求解 为了表示方便，我们将为了表示方便，我们将2.1节问题中所给表节问题中所给表6-1中的中的数据全部化为国际标准单位（表数据

32、全部化为国际标准单位（表6-2），时间用小时），时间用小时（h），高度用米（），高度用米（m）：）：表表6-2 一天内水塔水位记录一天内水塔水位记录(单位：米单位：米)时间时间/h水位水位/m时间时间/h水位水位/m09.6812.9510.210.929.4513.889.941.849.3114.989.652.959.1315.909.413.878.9816.839.184.988.8117.938.925.908.6919.048.667.008.5219.968.437.938.3920.848.228.978.2222.02水泵开动水泵开动9.98水泵开动水泵开动22.96水泵开

33、动水泵开动10.93水泵开动水泵开动23.8810.5910.9510.8224.9910.3512.0310.5025.9110.181模型假设模型假设（1）流量只取决于水位差，与水位本身无关，故由物）流量只取决于水位差，与水位本身无关，故由物理学中理学中Torriceli定律：小孔流出的液体的流速正比于定律：小孔流出的液体的流速正比于水面高度的平方根。问题中给出水塔的最低和最高水水面高度的平方根。问题中给出水塔的最低和最高水位分别是位分别是8.1648m)3024. 027(和和10.7352m)3024. 050.35(（设出口的水位为零），因为（设出口的水位为零），因为sqrt1467

34、. 1)1648. 87352.10(，约为，约为1，所以可忽略水位对流速的影响。，所以可忽略水位对流速的影响。（2）将流量看作是时间的连续函数，为计算简单，）将流量看作是时间的连续函数，为计算简单，不妨将不妨将流量定义成单位时间流出水的高度流量定义成单位时间流出水的高度，即水位对，即水位对时间变化率的绝对值（水位是下降的），水塔截面积时间变化率的绝对值（水位是下降的），水塔截面积为为3475.2334*)3024. 0*57(2S（m2），得到结果），得到结果后乘以后乘以s即可得真实的流量。即可得真实的流量。 2流量估计方法流量估计方法 首先依照表首先依照表6-2所给数据，用所给数据，用MA

35、TLAB作出时间作出时间水水位散点图（图位散点图（图6.2）。）。下面来计算水箱下面来计算水箱流量流量与与时间时间的关系。的关系。 根据图根据图6.2，一种简单的处理方法为，将表，一种简单的处理方法为，将表6-2中的中的数据分为三段，然后对每一段的数据做如下处理：数据分为三段，然后对每一段的数据做如下处理： 设某段数据设某段数据),(,),(),(1100nnyxyxyx，相邻数据中点的平均流速用下面的公式（流速，相邻数据中点的平均流速用下面的公式（流速=（左（左端点的水位右端点的水位）端点的水位右端点的水位）/区间长度）：区间长度）：iiiiiixxyyxxv111)2(每段数据首尾点的流速

36、用下面的公式计算：每段数据首尾点的流速用下面的公式计算：)()43()(022100 xxyyyxv)()43()(221nnnnnnxxyyyxv用以上公式求得时间与流速之间的数据如表用以上公式求得时间与流速之间的数据如表6-3。时间/h流速/cmh-1时间/h流速/ cmh-1029.8912.4931.520.4621.7413.4229.031.3818.4814.4326.362.39516.2215.4426.093.4116.3016.3724.734.42515.3217.3823.645.4413.0418.4923.426.4515.4519.5025.007.46513.

37、9820.4023.868.4516.3520.8422.178.9719.2922.02水泵开动9.98水泵开动22.96水泵开动10.93水泵开动23.8827.0910.9533.5024.4321.6211.4929.6325.4518.4825.9113.30由表由表6-3作出时间作出时间流速散点图如图流速散点图如图6.3。1）插值法）插值法 由表由表6-3，对水泵不工作时段，对水泵不工作时段1，2采取插值方法，可以得采取插值方法，可以得到任意时刻的流速，从而可以知道任意时刻的流量。到任意时刻的流速，从而可以知道任意时刻的流量。 我们分别采取拉格朗日插值法，分段线性插值法及我们分别采

38、取拉格朗日插值法，分段线性插值法及三次样条插值法；对于水泵工作时段三次样条插值法；对于水泵工作时段1应用前后时期应用前后时期的流速进行插值，由于最后一段水泵不工作时段数据的流速进行插值，由于最后一段水泵不工作时段数据太少，我们将它与水泵工作时段太少，我们将它与水泵工作时段2合并一同进行插值合并一同进行插值处理（该段简称混合时段）。处理（该段简称混合时段）。 我们总共需要对四段数据（第我们总共需要对四段数据（第1，2未供水时段，第未供水时段，第1供水时段，混合时段）进行插值处理，下面以第供水时段，混合时段）进行插值处理，下面以第1未供未供水时段数据为例分别用三种方法算出水时段数据为例分别用三种方

39、法算出流量函数流量函数和和用水用水量量（用水高度）。（用水高度）。 下面是用下面是用MATLAB实现该过程的程序实现该过程的程序(第第1未供水时段未供水时段)。 t=0,0.46,1.38,2.395,3.41,4.425,12.44,6.45,7.465,8.45,8.97；v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.29； t0=0:0.1:8.97；lglr=lagrange(t,v,t0); /*注：注：lagrange为一函数，程序同为一函数，程序同 lagrange.m*/lglrjf=trapz

40、(t0,lglr)fdxx=interpl(t,v,t0)；fdxxjf=trapz(t0,fdxx)scyt=interpl(t,v,t0,spline)；sancytjf=trapz(t0,scyt)plot(t,v,o,t0,lglr,r,t0,fdxx,g,t0,scyt,b)gtext(lglr)gtext(fdxx)gtext(spline)运行结果lglrjf=145.6231 ； fdxxjf=147.1430； sancytjf=145.6870 图图6.4是对第是对第1未供水段数据用三种不同方法得到的未供水段数据用三种不同方法得到的插值函数图，插值函数图， 图中曲线图中曲线

41、lglr、fdxx和和spline分别表示用拉格朗日插分别表示用拉格朗日插值法，分段线性插值法及三次样条插值法得到的曲线。值法，分段线性插值法及三次样条插值法得到的曲线。 由表由表6-2知，第知，第1未供水时段的总用水高度为未供水时段的总用水高度为146(=968822)，可见上述三种插值方法计算的结果，可见上述三种插值方法计算的结果与实际值（与实际值（146）相比都比较接近。考虑到三次样条插）相比都比较接近。考虑到三次样条插值方法具有更加良好的性质，建议采取该方法。值方法具有更加良好的性质，建议采取该方法。 其他三段的处理方法与第其他三段的处理方法与第1未供水时段的处理方法未供水时段的处理方

42、法类似，这里不再详细叙述，只给出数值结果和函数图类似，这里不再详细叙述，只给出数值结果和函数图像（图像（图6.5图图6.7），图中曲线标记同图），图中曲线标记同图6.4。 图图6.5 第一供水段时间第一供水段时间流速示意图流速示意图图图12.6 第第2未供水段时间未供水段时间流速示意图流速示意图图图6.7 混合时段时间混合时段时间流速示意图流速示意图图图6.8是用分段线性及三次样条插值方法得到的整个过是用分段线性及三次样条插值方法得到的整个过程的时间程的时间流速函数示意图。流速函数示意图。表表6-4 各时段及一天的总用水量（用水高度）各时段及一天的总用水量（用水高度）第第1未供水段未供水段第第

43、2未供水段未供水段第第3供水段供水段混合时段混合时段全天全天拉格朗日拉格朗日插值法插值法145.6231258.866454.068992.1337550.6921分段线性分段线性插值法插值法147.1430258.969749.605176.4688532.1866三次样条三次样条插值法插值法145.6870258.654753.333481.7699539.4450 表表6-5是对一天中任取的是对一天中任取的4个时刻分别用个时刻分别用3种方法得到种方法得到的水塔水流量近似值。的水塔水流量近似值。时间时间6.8810.8815.8822.8815.9826671234851433.74260

44、09085346325.5662241818047734.7099621055169414.8272413793103432.9976262626262725.4465591397849525.471578947368415.0527858196582033.7089553595325925.5490892055742329.41733175863551注：注：拉格朗日插值法拉格朗日插值法分段线性插值法分段线性插值法三次样条插值法三次样条插值法2）拟合法）拟合法（1）拟合水位）拟合水位时间函数时间函数 从表从表6-2中的测量记录看，一天有两次供水时段和三中的测量记录看，一天有两次供水时段和三次

45、未供水时段，分别对第次未供水时段，分别对第1，2未供水时段的测量数据未供水时段的测量数据直接作多项式拟合，可得到水位函数（注意，根据多直接作多项式拟合，可得到水位函数（注意，根据多项式拟合的特点，此处拟合多项式的次数不宜过高，项式拟合的特点，此处拟合多项式的次数不宜过高，一般以一般以36次为宜）。对第次为宜）。对第3未供水时段来说，数据未供水时段来说，数据过少不能得到很好的拟合。过少不能得到很好的拟合。 设设t，h分别为已输入的时刻和水位测量记录（由表分别为已输入的时刻和水位测量记录（由表6.2提供，水泵启动的提供，水泵启动的4个时刻不输入），这样第个时刻不输入），这样第1未未供水时段各时刻的

46、水位可由如下供水时段各时刻的水位可由如下MATLAB程序完成：程序完成：t=0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,23.88, 24.99,25.91h=9.68,9.45,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,10.35,10.18；c1=

47、polyfit(t(1:10),h(1:10),3)；% 用用3次多项式拟合次多项式拟合tp1=0:0.1:8.9；x1=polyval(c1,tp1)；plot(tp1,x1) 图图6.9给出的是第给出的是第1未供水时段的时间未供水时段的时间水位拟合水位拟合函数图形。函数图形。变量变量x1中存放了以中存放了以0.1为步长算出的各个时刻的水位高为步长算出的各个时刻的水位高度。同样地，第度。同样地，第2未供水时段时间未供水时段时间水位图可由如下水位图可由如下MATLAB程序完成，读者可自己上机运行查看。程序完成，读者可自己上机运行查看。c2=polyfit(t(11:21),h(11:21),3

48、)；tp2=10.9:0.1:20.9；x2=-polyval(c2,tp2)；plot(tp2,x2)（2）确定流量）确定流量时间函数时间函数对于第对于第1，2未供水时段的流量可直接对水位函数求导未供水时段的流量可直接对水位函数求导，程序如下：，程序如下：c1=polyfit(t(1:10),h(1:10),3)；c2=polyfit(t(11:21),h(11:21),3)；a1=polyder(c1)；a2=polyder(c2)； % 对拟合得到的水位函数求导，得流量对拟合得到的水位函数求导，得流量tp1=0:0.01:8.97；tp2=10.95:0.01:20.84；x13=-po

49、lyval(a1,tp1)； % Water_total1=100*trapz(tp1,x13); % 计算第计算第1未供水时段的总用水量未供水时段的总用水量x14=-polyval(a1,7.93,8.97)； % 求某些点的流量求某些点的流量 x23=-polyval(a2,tp2)；Water_total2=100*trapz(tp2,x114); %计算第计算第2未供水时段的总用水量未供水时段的总用水量x24=-polyval(a2,10.95,12.03)； x25=-polyval(a2,19.96,20.84)；subplot(1,2,1)plot(tp1,x13*100) %与

50、图与图6.10单位保持一致单位保持一致 % 倍数倍数100是使单位由米变成厘米是使单位由米变成厘米subplot(1,2,2)plot(tp2,x23*100) %与图与图6.10单位保持一致单位保持一致 程序运行得到第程序运行得到第1，2未供水时段的时间未供水时段的时间流量图如图流量图如图6.10，可以看到与图，可以看到与图6.8中用插值给出的曲线比较吻合。中用插值给出的曲线比较吻合。 如果用如果用5次多项式拟合则得图次多项式拟合则得图6.11的图形，显然较的图形，显然较三次拟合的效果好。三次拟合的效果好。 而第而第1供水时段的流量则用前后时期的流量进行拟供水时段的流量则用前后时期的流量进行

51、拟合得到。为使流量函数在合得到。为使流量函数在t=9和和t=11连续，我们只取连续，我们只取4个点，用三次多项式拟合得到第个点，用三次多项式拟合得到第1供水时段的时间供水时段的时间流量图形如图流量图形如图6.12，可以看到与图，可以看到与图6.8中的相应部分中的相应部分比较吻合。比较吻合。图图6.12图图6.8程序如下：程序如下： 在第在第2供水时段之前取供水时段之前取t=19.96,20.84两点的流量，两点的流量，用第用第3未供水时段的未供水时段的3个记录做差分得到两个流量数个记录做差分得到两个流量数据据21.62，18.48，然后用这，然后用这4个数据做三次多项式拟个数据做三次多项式拟合

52、得到第合得到第2供水时段与第供水时段与第3未供水时段的时间未供水时段的时间流量流量图如图图如图6.13，可以看到与图，可以看到与图6.8中的相应部分也比较中的相应部分也比较吻合吻合。图图6.13，图图6.8程序如下：程序如下：t3=19.96,20.84,t(22),t(23)；x3=x25*100,21.62,18.48；x_flow3=polyfit(t3,x3,3)；tp3=19.96:0.01:25.91；xx3=polyval(x_flow3,tp3)；tp4=20.84:0.01:23.88；gs2=polyval(x_flow3,tp4)；gsysl2=trapz(tp4,gs2); %该语句计算第该语句计算第2供水时段的总用水量供水时段的总用水量figure;plot(tp3,xx3)（3）一天总用水量的估计）一天总用水量的估计 分别对供水的两个时段和不供水的两个时段积分分别对供水的两个时段和不供水的两个时段积分（流量对时间）并求和得到一天的总用水量约为（流量对时间）并求和得到一天的总用水量约为526.8935（此数据是总用水高度，单位为（此数据是总用水高度，单位为cm）。表）。表6-6列出各段用水量，与插值法算得的表列出各段用水量，与插值法算得的表6-4相比，二相比，二者较为吻合。者较为吻合。表表6-6时段时段第第