数学建模

1水流量的估计美国某州的用水管理机构要求各社区提供以每小时多少加仑计的用水量以及每天所用的总水量。许多社区没有测量流入或流出水塔水量的装置，只能代之以每小时测量水塔中的水位，其误差不超过5%。需要注意的是，当水塔中的水位下降到最低水位L时，水泵就自动向水塔输水直到最高水位H，此期间不能测量水泵的供水量，因此，当水泵正在输水时不容易建立水塔中水位和用水量之间的关系。水泵每天输水一次或两次，每次约2小时。试估计任何时刻（包括水泵正在输水时间）从水塔流出的水流量f(t)，并估计一天的总用水量。已知该水塔是一个高为40ft（英尺），直径为57ft（英尺）的正圆柱，表5-1给出了某个小镇一天水塔水位的真实数据，水位降至约27.00ft水泵开始工作，水位升到35.50ft停止工作。（注：1ft（英尺）=0.3048m（米））题目：2某小镇某天水塔水位（表一）时间/s水位/0.01ft时间/s水位/0.01ft03175466363350331631104995332606635305453936316710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254水泵开动35932水泵开动82649水泵开动39332水泵开动8596834753943535508995333974331834459327033403一、问题分析本题所指流量可视为单位时间内流出水的体积。由于水塔是正圆柱形，横截面积是常数，所以在水泵不工作时段，流量很容易根据水位相对时间的变化算出。问题的难点在于如何估计水泵供水时段的流量。水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到。作为用于插值或拟合的原始数据，我们希望水泵不工作时段的流量越准确越好。在此流量大体上可由两种方法计算，一是直接对表1中的水量用数值微分算出各时段的流量，用它们拟合其它时刻或连续时间的流量；二是先用表中数据拟合水位一时间函数，求导数即可得到连续时间的流量。4有了任何时刻的流量，就能计算一天的总用水量。其实，水泵不工作时段的用水量可以由测量记录直接得到，由表1中下降水位乘以水搭的截面积就是这一时段的用水量。这个数值可以用来检验数据插值或拟合的结果（并且用MATLAB软件实现插值法）。为了表示方便，我们将问题中所给表1中的数据全部化为国际标准单位（表2），时间用小时（h），高度用米（m）：问题求解5时间/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.5910.9510.8224.9910.3512.0310.5025.9110.18表2一天内水塔水位记录6二、模型假设（1）流量只取决于水位差，与水位本身无关，故由物理学中Torriceli定律：小孔流出的液体的流速正比于水面高度的平方根。题目给出水塔的最低和最高水位分别是8.1648m)3024.027(和10.7352m)3024.050.35(（设出口的水位为零），因为sqrt1467.1)1648.87352.10(约为1，所以可忽略水位对流速的影响。（2）将流量看作是时间的连续函数，为计算简单，不妨将流量定义成单位时间流出水的高度，即水位对时间变化率的绝对值（水位是下降的），水塔截面积为3475.2334*)3024.0*57(2S（m2），最后公式得流量=H×S。72．流量估计方法首先依照表2所给数据，用MATLAB作出时间—水位散点图。（图1）8下面来计算水箱流量与时间的关系。根据图1，一种简单的处理方法为，将表2中的数据分为三段，然后对每一段的数据做如下处理：设某段数据)},(,),,(),,{(1100nnyxyxyx，相邻数据中点的平均流速用下面的公式（流速=（右端点的水位－右端点的水位）/区间长度）：iiiiiixxyyxxv111)2(每段数据首尾点的流速用下面的公式计算：)()43()(022100xxyyyxv)()43()(221nnnnnnxxyyyxv9时间/h流速/cm·h-1时间/h流速/cm·h-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.2922.02水泵开动9.98水泵开动22.96水泵开动10.93水泵开动23.8827.0910.9533.5024.4321.6211.4929.6325.4518.4825.9113.30用以上公式求得时间与流速之间的数据（表3）10由表3作出时间—流速散点图如图21插值法由表3，对水泵不工作时段1，2采取插值方法，可以得到任意时刻的流速，从而可以知道任意时刻的流量。11我们分别采取拉格朗日插值法，分段线性插值法及三次样条插值法；对于水泵工作时段1应用前后时期的流速进行插值，由于最后一段水泵不工作时段数据太少，我们将它与水泵工作时段2合并一同进行插值处理（该段简称混合时段）。我们总共需要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理，下面以第1未供水时段数据为例分别用三种方法算出流量函数和用水量（用水高度）。12运行结果lglrjf=145.6231；fdxxjf=147.1430；sancytjf=1412.6870下面是用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；lglr=lglrcz(t,v,t0)；/*注：lglrcz为一函数，程序同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’)13图3是对第1未供水段数据用三种不同方法得到的插值函数图，图中曲线lglr、fdxx和scyt分别表示用拉格朗日插值法，分段线性插值法及三次样条插值法得到的曲线。14由表2知，第1未供水时段的总用水高度为146(=968－822)，可见上述三种插值方法计算的结果与实际值（146）相比都比较接近。考虑到三次样条插值方法具有更加良好的性质，我们采取该方法。其他三段的处理方法与第1未供水时段的处理方法类似，这里不再详细叙述，只给出数值结果和函数图像（图4～图6），图中曲线标记同图3。15图4第一供水段时间—流速示意图图5第2未供水段时间—流速示意图图6混合时段时间—流速示意图图7是用分段线性及三次样条插值方法得到的整个过程的时间—流速函数示意图。16表4各时段及一天的总用水量（用水高度）第1未供水段第2未供水段第3供水段混合时段全天拉格朗日插值法145.6231258.866454.068992.1337550.6921分段线性插值法147.1430258.969749.605176.4688532.1866三次样条插值法145.6870258.654753.333481.7699539.445017表5是对一天中任取的4个时刻分别用3种方法得到的水塔水流量近似值。时间6.8810.8815.8822.88①15.9826671234851433.7426009085346325.5662241818047734.70996210551694②14.8272413793103432.9976262626262725.4465591397849525.4715789473684③15.0527858196582033.7089553595325925.5490892055742329.41733175863551注：①拉格朗日插值法②分段线性插值法③三次样条插值法18拟合法（1）拟合水位—时间函数从表2中的测量记录看，一天有两次供水时段和三次未供水时段，分别对第1，2未供水时段的测量数据直接作多项式拟合，可得到水位函数（另外根据多项式拟合的特点，拟合多项式的次数不宜过高，所以我们们用了3～6次）。对第3未供水时段来说，数据过少不能得到很好的拟合。设t，h分别为已输入的时刻和水位测量记录（由表2提供，水泵启动的4个时刻不输入）19t=[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.66]h=[9.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)这样第1未供水时段各时刻的水位可由如下MATLAB程序完成：20图8给出的是第1未供水时段的时间—水位拟合函数图形。21变量x1中存放了以0.1为步长算出的各个时刻的水位高度。同样地，第2未供水时段时间—水位图可由如下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未供水时段的流量可直接对水位函数求导，程序如下：22c1=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(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])；*/为下面的程序准备数据/*x25=-polyval(a2,[19.96,20