水流量的估计最新课件

1、水流量的估计最新课件1 水流量的估计12.1 实验目的实验目的 12.2 实验内容实验内容 美国某州的用水管理机构要求各社区提供以每小时多少加仑计的用水量以及每天所用的总水量。许多社区没有测量流入或流出水塔水量的装置，只能代之以每小时测量水塔中的水位，其误差不超过5%。需要注意的是，当水塔中的水位下降到最低水位L时，水泵就自动向水塔输水直到最高水位H，此期间不能测量水泵的供水量，因此，当水泵正在输水时不容易建立水塔中水位和用水量之间的关系。水泵每天输水一次或两次，每次约2小时。 本实验的主要目的是使学生学会用本实验的主要目的是使学生学会用MATLAB软件进行插值软件进行插值和曲线拟合计算并解决

2、一些具体的实际问题。通过实际问题的和曲线拟合计算并解决一些具体的实际问题。通过实际问题的解决，使学生了解如何利用曲线插值及曲线拟合解决实际问题解决，使学生了解如何利用曲线插值及曲线拟合解决实际问题的全过程。的全过程。2.1 实验问题实验问题水流量的估计最新课件2 试估计任何时刻（包括水泵正在输水时间）从水塔流出的水流量f(t)，并估计一天的总用水量。已知该水塔是一个高为40ft（英尺），直径为57ft（英尺）的正圆柱，表5-1给出了某个小镇一天水塔水位的真实数据，水位降至约27.00ft水泵开始工作，水位升到35.50ft停止工作。（注：1ft（英尺）=0.3048m（米）水流量的估计最新课件

3、3表12.1 某小镇某天水塔水位时间/s水位/0.01ft时间/s水位/0.01ft03175466363350331631104995332606635305453936316710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254水泵开动35932水泵开动82649水泵开动39332水泵开动859683475394353550899533397433183445932703340水流量的估计最新课件42

4、.2 问题分析 本实验所指流量可视为单位时间内流出水的体积。由于水塔是正圆柱形，横截面积是常数，所以在水泵不工作时段，流量很容易根据水位相对时间的变化算出。问题的难点在于如何估计水泵供水时段的流量。 水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到。作为用于插值或拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。这此流量大体上可由两种方法计算，一是直接对表12-1中的水量用数值微分算出各时段的流量，用它们拟合其它时刻或连续时间的流量；二是先用表中数据拟合水位一时间函数，求导数即可得到连续时间的流量。水流量的估计最新课件5 有了任何时刻的流量，就不难计算一天的总用水量。其实，水泵

5、不工作时段的用水量可以由测量记录直接得到，由表12-1中下降水位乘以水搭的截面积就是这一时段的用水量。这个数值可以用来检验数据插值或拟合的结果。 在具体给出本问题的解答之前，先介绍一个简单的数据插值方法。水流量的估计最新课件62.3 拉格朗日插值拉格朗日插值1、线性插值)(xf假设已知在区间,ba上的两个结点10, xx和它们的函数值 10 xxx )()()(10 xfxfxf求一个一次多项式xaaxP101)(，使得多项式)(1xP在结点上满足条件 1, 0),()(1ixfxPii 这种插值方法称为线性插值方法（也称两点插值）。 可以求出： )()()(101001011xfxxxxxf

6、xxxxxP水流量的估计最新课件72、抛物插值)(xf已知在区间,ba上的三个结点210,xxx和它们的函数值 210 xxxx )()()()(210 xfxfxfxf求一个次数不超过2的多项式22102)(xaxaaxP，使得它在结点上满足条件2, 1, 0),()(2ixfxPii这种插值方法称为抛物线插值法， 可求出：)()( )()( )()()( )()( )()(121012002010212xfxxxxxxxxxfxxxxxxxxxP)()( )()( )(2021210 xfxxxxxxxx水流量的估计最新课件83、n次拉格朗日插值,ba假设取区间 上的1n个结点nxxx,1

7、0，并且已知函数)(xf在这此点的函数值 nxxxxx210 )()()()()(210nxfxfxfxfxf现在求一个次数不超过n的多项式nnnxaxaxaaxP2210)(，使得满足条件 nixfxPiin, 1, 0),()( 这种插值方法称为n次多项式插值（或称代数插值）， 利用拉格朗日插值插值方法可得水流量的估计最新课件9)()()()()()()()()(,2, 21, 10, 0nnnnnnnxfxlxfxlxfxlxfxlxP0( )0,1,njknjkjj kxxlxknxx其中上述多项式称为 n次拉格朗日（Lagrange）插值多项式， 函数nkxlnk, 1, 0),(,

8、称为拉格朗日插值基函数。当n=1,2时，n次拉格朗日（Lagrange）插值多项式即为线性插值多项式和抛物插值多项式。水流量的估计最新课件104321x例12.1 已知函数发f(x)的函数表如下： 31 . 28 . 15 . 2)(xf求其拉格朗日插值多项式，并求)5 . 2(f的近似值。 解 由于给出了4个插值结点，所以可做出次数不超过3的拉格朗日插值多项式。43132361)41 ( )31 ( )21 ()4( )3( )2()(233 , 0 xxxxxxxl6219421)42( )32( )22()4( )3( )2()(233 , 1xxxxxxxl472721)43( )33

9、( )23()4( )3( )2()(233 , 2xxxxxxxl水流量的估计最新课件11161161)44( )34( )24()4( )3( )2()(233 , 3xxxxxxxl 将上列4式代入n=3 的拉格朗日插值公式，可得所要求的插值多项式为6 . 4933. 29 . 00667. 0)(233xxxxP将x=2.5代入可得f(2.5)的近似值为1.8496。 拉格朗日插值法适合节点较少的情况，当节点较多的大范围高次插值的逼近效果往往并不理想且当插值结点增加时，计算越来越繁。为了提高精度和减少计算，还有牛顿插值法下、三次样条插值等，具体可参阅有关书籍。水流量的估计最新课件121

10、一维插值命令一维插值命令interp1的具体使用格式的具体使用格式yy=interp1(x,y,xx,method) 其中x,y是插值结点的数据向量，如果y是矩阵，则对矩阵y的每一列相对x进行插值，xx是待求函数值的插值结点向量，可以缺省。method是可选项，说明插值使用的方法。对于一维插值，MATLAB提供可选的方法有：nearest，linear，spline，cubic，它们分别表示最近插值、线性插值、三次样条插值和三次插值。2.4 MATLAB软件实现插值法软件实现插值法 MATLAB软件提供了专门做各种插值的命令：interp1（一维插值），interp2（二维插值），interp

11、3（三维插值），interpn（n维插值），spline（样条插值）等。 水流量的估计最新课件132二维插值命令二维插值命令interp2的具体使用格式的具体使用格式zz=interp2(x,y,z,xx,yy,method) 该指令的意思是根据数据向量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

12、(z,ntimes)水流量的估计最新课件143三维插值命令三维插值命令interp3的具体使用格式的具体使用格式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,cubic)水流量的估计最新课件15例12.2 在用外接电源给电容器充电时，电容器两端的电压V将会随着充电时间t发生变化，已知在某一次实验时

13、，通过测量得到下列观测值，分别用拉格朗日插值法、分段线性插值法、三次样条插值法画出V随着时间t变化的曲线图，分别计算当时间t=7s时，三种插值法各自算得电容器两端电压的近似值。1295 . 64321st 4 .101 .106 . 90 . 92 . 83 . 72 . 6vV解 由于MATLAB没有提供现成的拉格朗日插值命令，我们可以编写一个函数lglrcz.m来完成，其他两种插值法可用现成的命令。 用MATLAB软件进行三种插值计算的程序为szczqx.m。水流量的估计最新课件16程序程序lglrcz.m：function y=lglrcz(x0,y0,x)n=length(x0)；m=

14、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*y0(k)+s； end y(i)=s； end水流量的估计最新课件17程序为azczf.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=lglrcz(t,v,t0) ；laglr=lglrcz(t,v,7) ；fdxx=interp1(t,v,t0) ；fendxx=int

15、erp1(t,v,7) ；scyt=interp1(t,v,t0,spline) ；sancyt=interp1(t,v,7,spline)plot(t,v,*,t0,lglr,r,t0,fdxx,g,t0,scyt,b)gtext(lglr)gtext(fdxx)gtext(scyt)水流量的估计最新课件18执行结果是laglr=9.52988980716254fendxx=9.700sancyt=9.671图形如图12.1所示。图中曲线lglr表示拉格朗日插值曲线，scyt表示三次样条插值曲线，fdxx表示分段线性插值曲线。水流量的估计最新课件19 2.5 问题求解问题求解 为了表示方便，

16、我们将2.1节问题中所给表12-1中的数据全部化为国际标准单位（表12-2），时间用小时（h），高度用米（m）：表12-2 一天内水塔水位记录水流量的估计最新课件20时间/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水泵开动10.93水泵开动23.8810

17、.5910.9510.8224.9910.3512.0310.5025.9110.18水流量的估计最新课件211模型假设（1）流量只取决于水位差，与水位本身无关，故由物理学中Torriceli定律：小孔流出的液体的流速正比于水面高度的平方根。题目给出水塔的最低和最高水位分别是8.1648m)3024. 027(和10.7352m)3024. 050.35(（设出口的水位为零），因为sqrt1467. 1)1648. 87352.10(，约为1，所以可忽略水位对流速的影响。（2）将流量看作是时间的连续函数，为计算简单，不妨将流量定义成单位时间流出水的高度，即水位对时间变化率的绝对值（水位是下降的

18、），水塔截面积为3475.2334*)3024. 0*57(2S（m2），得到结果后乘以s即可。 水流量的估计最新课件222流量估计方法 首先依照表12-2所给数据，用MATLAB作出时间水位散点图（图12.2）。水流量的估计最新课件23下面来计算水箱流量与时间的关系。 根据图12.，一种简单的处理方法为，将表12-2中的数据分为三段，然后对每一段的数据做如下处理： 设某段数据),(,),(),(1100nnyxyxyx，相邻数据中点的平均流速用下面的公式（流速=（右端点的水位右端点的水位）/区间长度）：iiiiiixxyyxxv111)2(每段数据首尾点的流速用下面的公式计算：)()43()

19、(022100 xxyyyxv)()43()(221nnnnnnxxyyyxv用以上公式求得时间与流速之间的数据如表12-3。水流量的估计最新课件24时间/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.9820.4023.868.4516.3520.8422.178.9719.292

20、2.02水泵开动9.98水泵开动22.96水泵开动10.93水泵开动23.8827.0910.9533.5024.4321.6211.4929.6325.4518.4825.9113.30水流量的估计最新课件25由表12-3作出时间流速散点图如图12.3。1）插值法 由表12-3，对水泵不工作时段1，2采取插值方法，可以得到任意时刻的流速，从而可以知道任意时刻的流量。 水流量的估计最新课件26 我们分别采取拉格朗日插值法，分段线性插值法及三次样条插值法；对于水泵工作时段1应用前后时期的流速进行插值，由于最后一段水泵不工作时段数据太少，我们将它与水泵工作时段2合并一同进行插值处理（该段简称混合时

21、段）。 我们总共需要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理，下面以第1未供水时段数据为例分别用三种方法算出流量函数和用水量（用水高度）。 下面是用MATLAB实现该过程的程序。 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；水流量的估计最新课件27lglr=lglrcz(t,v,t0)； /*注：注：lglrcz为一函数，程序同为一

22、函数，程序同lglrcz.m*/lglrjf=0.1*trapz(lglr)fdxx=interpl(t,v,t0)；fdxxjf=0.1*trapz(fdxx)scyt=interpl(t,v,t0,spline)；sancytjf=0.1*trapz(scyt)plot(t,v,*,t0,lglr,r,t0,fdxx,g,t0,scyt,b)gtext(lglr)gtext(fdxx)gtext(scyt)运行结果lglrjf=145.6231 ； fdxxjf=147.1430； sancytjf=1412.6870水流量的估计最新课件28 图12.4是对第1未供水段数据用三种不同方法得

23、到的插值函数图， 图中曲线lglr、fdxx和scyt分别表示用拉格朗日插值法，分段线性插值法及三次样条插值法得到的曲线。 水流量的估计最新课件29 由表12-2知，第1未供水时段的总用水高度为146(=968822)，可见上述三种插值方法计算的结果与实际值（146）相比都比较接近。考虑到三次样条插值方法具有更加良好的性质，建议采取该方法。 其他三段的处理方法与第1未供水时段的处理方法类似，这里不再详细叙述，只给出数值结果和函数图像（图12.5图12.7），图中曲线标记同图12.4。 水流量的估计最新课件30图12.5 第一供水段时间流速示意图水流量的估计最新课件31图12.6 第2未供水段时

24、间流速示意图水流量的估计最新课件32图12.7 混合时段时间流速示意图水流量的估计最新课件33图12.8是用分段线性及三次样条插值方法得到的整个过程的时间流速函数示意图。水流量的估计最新课件34表12-4 各时段及一天的总用水量（用水高度）第1未供水段第2未供水段第3供水段混合时段全天拉格朗日插值法145.6231258.866454.068992.1337550.6921分段线性插值法147.1430258.969749.605176.4688532.1866三次样条插值法145.6870258.654753.333481.7699539.4450水流量的估计最新课件35 表12-5是对一天

25、中任取的4个时刻分别用3种方法得到的水塔水流量近似值。时间6.8810.8815.8822.8815.9826671234851433.7426009085346325.5662241818047734.79414.8272413793103432.9976262626262725.4465591397849525.4715.82033.75925.54929.41733175863551注：拉格朗日插值法分段线性插值法三次样条插值法水流量的估计最新课件362）拟合法（1）拟合水位时间函数 从表12-2中的测量记录看，一天有两次供水时段和三次未供水时段，分别对第1，2未供水时段的测量数据直接作

26、多项式拟合，可得到水位函数（注意，根据多项式拟合的特点，此处拟合多项式的次数不宜过高，一般以36次为宜）。对第3未供水时段来说，数据过少不能得到很好的拟合。 设t，h分别为已输入的时刻和水位测量记录（由表12.2提供，水泵启动的4个时刻不输入），这样第1未供水时段各时刻的水位可由如下MATLAB程序完成：水流量的估计最新课件37t=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.8824.99,25.66h=9.

27、68,9.48,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=polyfit(t(1:10),h(1:10),3)；tp1=0:0.1:8.9；x1=polyval(c1,tp1)；plot(tp1,x1)水流量的估计最新课件38 图12.9给出的是第1未供水时段的时间水位拟合函数图形。水流量的估计最新课件39变量x1中存放了以0.1为步长算出的各个时刻的水位高度。同样地，第2未供水时段时间水位图可由如下

28、MATLAB程序完成，读者可自己上机运行查看。c2=polyfit(t(11:21),h(11:21),3)；tp2=10.9:0.1:20.9；x2=-polyval(c2,tp2)；plot(tp2,x2)（2）确定流量时间函数对于第1，2未供水时段的流量可直接对水位函数求导，程序如下：水流量的估计最新课件40c1=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=-polyval(

29、a1,tp1)；x113=-polyval(a1,0:0.01:8.97)；wgsysl1=100*trapz(tp1,x113)； */计算第计算第1未供水时段的总用水量未供水时段的总用水量/*x14=-polyval(a1,7.93,8.97)； */为下面的程序准备数据为下面的程序准备数据/*x23=-polyval(a2,tp2)；x114=-polyval(a2,10.95:0.01:20.84)wgsysl2=100*trapz(tp2,x114)； */计算第计算第2未供水时段的总用水量未供水时段的总用水量/*x24=-polyval(a2,10.95,12.03)； */为下面

30、的程序准备数据为下面的程序准备数据/*x25=-polyval(a2,19.96,20.84)； */为下面的程序准备数据为下面的程序准备数据/*subplot(1,2,1)plot(tp1,x13*100) */与图与图12.10单位保持一致单位保持一致/*subplot(1,2,2)plot(tp2,x23*100) */与图与图12.10单位保持一致单位保持一致/*水流量的估计最新课件41 程序运行得到第1，2未供水时段的时间流量图如图12.10，可以看到与图12.8中用插值给出的曲线比较吻合。水流量的估计最新课件42 如果用5次多项式拟合则得图12.11的图形，显然较三次拟合的效果好。

31、水流量的估计最新课件43 而第1供水时段的流量则用前后时期的流量进行拟合得到。为使流量函数在t=9和t=11连续，我们只取4个点，用三次多项式拟合得到第1供水时段的时间流量图形如图12.12，可以看到与图12.8中的相应部分比较吻合。图12.12图12.8水流量的估计最新课件44dygsdsy=7.93,8.97,10.95,12.03；dygsdls=x14,x24；nhjg=polyfit(dygsdsj,dygsdls,3)；nhsj=7.93:0.1:12.03；nhlsjg=polyval(nhjg,nhsj)；gssj1=8.97:0.01:10.95；gs1=polyval(nh

32、jg,8.97:0.01:10.95)；gsysl1=100*trapz(gssj1,gs1)； */该语句计算第1供水时段的总用水量/*plot(nhsj,100*nhlsjg) 程序如下：水流量的估计最新课件45 在第2供水时段之前取t=19.96,20.84两点的流量，用第3未供水时段的3个记录做差分得到两个流量数据21.62，18.48，然后用这4个数据做三次多项式拟合得到第2供水时段与第3未供水时段的时间流量图如图12.13，可以看到与图12.8中的相应部分也比较吻合。图12.13，图12.8水流量的估计最新课件46程序如下：t3=19.96,20.84,t(22),t(23)；ls3=x25*100,21.62,18.48；nhhddxsxs=polyfit(t3,ls3,3