水塔水流量估計2009

1、實驗 水塔水流量的估計 實驗?zāi)康谋敬螌嶒灥闹饕康氖亲寣W(xué)生會用數(shù)學(xué)軟件進行插值計算并解決一些具體的實際問題。介紹一些經(jīng)典的插值方法，包括拉格朗日插值法、埃爾米特插值法、分段插值法、三次樣條插值法等等。實驗內(nèi)容1實驗問題 美國某州的各用水管理機構(gòu)要求各社區(qū)提供以每小時多少加侖計的用水率以及每天所用的總水量。許多社區(qū)沒有測量流入或流出水塔的水量裝置,他們只能代之以每小時測量水塔中的水位,其誤差不超過5%。更重要的是,當(dāng)水塔中的水位下降到最低水位L時水泵就啟動向水塔輸水直到最高水位H,期間不能測量水泵的供水量。因此,當(dāng)水泵正在輸水時不容易建立水塔中水位和用水量之間的關(guān)系。水泵每天輸水一次或兩次,每次

2、約二小時。試估計任何時刻(包括水泵正在輸水時間)從水塔流出的水流量f(t),并估計一天的總用水量。已知該水塔是一個高為40英尺（ft），直徑為57英尺(ft)的正圓柱，表12.1給出了某個小鎮(zhèn)一天水塔水位的真實數(shù)據(jù)，水位降至約27.00ft水泵開始工作,水位升到35.50ft停止工作。(注：1英尺（ft）=0.3024米（m）)表12-1 某小鎮(zhèn)某天水塔水位時間/s水位/0.01ft時間/s水位/0.01ft0317546636335033163110499533260663530545393631671061929945725430871393729476057430121792128926

3、4554292721240285068535284225223279571854276728543275275021269732284269779254水泵開動35932水泵開動82649水泵開動39332水泵開動8596834753943535508995333974331834459327033402 問題分析 流量是單位時間內(nèi)流出水的體積,由于水塔是正圓柱形,橫截面積是常數(shù),所以在水泵不工作時段,流量很容易根據(jù)水位相對時間的變化率算出。問題的難點在于如何估計水泵供水時段的流量。 水泵供水時段的流量只能靠供水時段前后的流量經(jīng)插值或擬合得到。作為用于插值或擬合的原始數(shù)據(jù),我們希望水泵不工作時

4、段的流量越準確越好。這些流量大體上可由兩種方法計算,一是直接對表12.1中的水量用數(shù)值微分算出各時段的流量,用它們擬合其它時刻或連續(xù)時間的流量;二是先用表中數(shù)據(jù)擬合水位時間函數(shù),求導(dǎo)數(shù)即可得到連續(xù)時間的流量。有了任何時刻的流量,就不難計算一天的總用水量。其實,水泵不工作時段的用水量可以由測量記錄直接得到,由表12.1中下降水位乘以水塔的截面積就是這一時段的量這個數(shù)值可以用來檢驗數(shù)據(jù)插值或擬合的結(jié)果。模型建立為了表示方便,我們將問題中所給表12-1中的數(shù)據(jù)全部化為國際標準單位（表12-5）,時間用小時(h),高度用米(m):表12-5 一天內(nèi)水塔水位記錄時間(h)水位(m)時間(h)水位(m)0

5、9.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模型假設(shè)(1) 流量只取決于水位差,與水位本身無關(guān)，故由物理學(xué)中Torricelli定律：從小孔流出的液體的流速正

6、比于水面高度的平方根。題目給出水塔的最低和最高水位分別是8.1648m（27×0.3024）和10.7352m（35.50×0.3024 ）(設(shè)出口的水位為零)，因為sqrt(10.7352/8.1648)1.1467，約為1,所以可忽略水位對流速的影響。（2） 將流量看作是時間的連續(xù)函數(shù)，為計算簡單,不妨將流量定義成單位時間流出水的高度,即水位對時間變化率的絕對值(水位是下降的), 水塔截面積為S=(57*0.3048)2*/4=237.8(m2).得到結(jié)果后乘以S即可。流量估計方法首先依照表12-5所給數(shù)據(jù),用MATLAB作出時間水位散點圖12.7。 Fig.12.7下

7、面我們來計算水箱流量與時間的關(guān)系。根據(jù)表12-5中數(shù)據(jù)散點圖12.7，一種簡單的處理方法為先將表12-5中的數(shù)據(jù)分為三段，然后對每一段的數(shù)據(jù)做如下處理：設(shè)某段數(shù)據(jù)為，相鄰數(shù)據(jù)中點的平均流速用下面的公式(流速=（左端點的水位-右端點的水位）/區(qū)間長度)：，每段數(shù)據(jù)首尾點的流速用下面的公式計算：，用以上公式求得時間與流速之間的數(shù)據(jù)表12-6如下: 表12-6時間(h)流速(cm/h)時間(h)流速(cm/h)029.8912.4931.520.4621.7413.4229.031.3818.4814.4326.362.39516.2215.4426.093.4116.3016.3724.734.4

8、2515.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由表12-6作出時間-流速散點圖如下: Fig.12.8 (1) 插值法由表12-6,對水泵不工作時段1,2采取插值方法,可以得到任意時刻的流速,從而可以知道任意時刻的流量.我們分別采取拉格朗日插值法,分

9、段線性插值法及三次樣條插值法；對于水泵工作時段1應(yīng)用前后時期的流速進行插值，由于最后一段水泵不工作時段數(shù)據(jù)太少，我們將它與水泵工作時段2合并一同進行插值處理(該段簡稱混合時段)。我們總共需要對四段數(shù)據(jù)（第1，2未供水時段，第1供水時段，混合時段）進行插值處理,下面以第1未供水時段數(shù)據(jù)為例分別用三種方法算出流量函數(shù)和用水量（用水高度）。下面是實現(xiàn)該過程的MATLAB程序。t=0,0.46,1.38,2.395,3.41,4.425,5.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,1

10、6.35,19.27;t0=0:0.1:8.97;lglr=lglrcz(t,v,t0); /*注:lglrcz為一函數(shù),程序同lglrcz.m*/ lglrjf=0.1*trapz(lglr)fdxx=interp1(t,v,t0);fdxxjf=0.1*trapz(fdxx)scyt=interp1(t,v,t0,'spline');sancytjf=0.1*trapz(scyt)plot(t,v,'*',t0,lglr,'r',t0,fdxx,'g',t0,scyt,'b')gtext('lglr&#

11、39;)gtext('fdxx')gtext('scyt')運行結(jié)果為 lglrjf=145.6231 fdxxjf=147.1430 sancytjf=145.6870圖12.9是對第1未供水段數(shù)據(jù)用三種不同方法得到的插值函數(shù)圖，圖中曲線lglr、fdxx和scyt分別表示用拉格朗日插值法,分段線性插值法及三次樣條插值法得到的曲線。 Fig.12.9由表12-5知，第1未供水時段的總用水高度為146（=968-822），可見上述三種插值方法計算的結(jié)果與實際值（146）相比都比較接近?？紤]到三次樣條插值方法具有更加良好的性質(zhì)，建議采取該方法。其它三段的處理方法與

12、第1未供水時段的處理方法類似，這里不再詳細敘述，只給出數(shù)值結(jié)果和函數(shù)圖象，圖中曲線標記同F(xiàn)ig.12.9。表12-7 各時段及一天的總用水量(用水高度)第1未供水段第2未供水段第1供水段混合時段全天拉格朗日插值法145.6231258.866454.068992.1337550.6921分段線性插值法147.1430258.969749.605176.4688532.1866三次樣條插值法145.6870258.654753.333481.7699539.4450Fig.12.10 (第一供水段時間-流速示意圖)Fig.12.11(第2未供水段時間-流速示意圖)Fig.12.12(混合時段段時

13、間-流速示意圖)下圖12.13是用分段線性及三次樣條插值方法得到的整個過程的時間-流速函數(shù)示意圖。Fig.12.13下表是對一天中任取的4個時刻分別用三種方法得到的水塔水流量近似值。表12-86.8810.8815.8822.88拉格朗日15.9826671234851433.7426009085346325.5662241818047734.70996210551694分段插值14.8272413793103432.9976262626262725.44655913978495三次樣條15.0527858196582033.7089553595325925.54908920557423(2)

14、 擬合法1) 擬合水位-時間函數(shù)從表12-5中測量記錄看，一天有兩次供水時段和三次未供水時段，分別對第1,2未供水時段的測量數(shù)據(jù)直接作多項式擬合，可得到水位函數(shù)（注意，根據(jù)多項式擬合的特點，此處擬合多項式的次數(shù)不宜過高，一般以3-6次為宜）。對第3未供水時段來說,數(shù)據(jù)過少不能得到很好的擬和。設(shè)t,h分別為已輸入的時刻和水位測量記錄(由表12-5提供，水泵啟動的4個時刻不輸入)，這樣第1未供水時段各時刻的水位可由如下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,1

15、5.90,16.83,17.93,19.04,19.96,20.84,23.88,24.99,25.66h=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未供水時段的時間-水位擬合函數(shù)圖形。變量x1中存放了以0.1為步長算出的各

16、個時刻的水位高度。 同樣地，第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) 確定流量時間函數(shù)對于第1,2未供水時段的流量可直接對水位函數(shù)求導(dǎo)，程序如下：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.

17、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); */該語句僅為下面的程序準備數(shù)據(jù)*/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); */該語句僅為下

18、面的程序準備數(shù)據(jù)*/x25=-polyval(a2,19.96,20.84); */該語句僅為下面的程序準備數(shù)據(jù)*/subplot(1,2,1)plot(tp1,x13*100) */與Fig.12.13單位保持一致*/subplot(1,2,2)plot(tp2,x23*100) */與Fig.12.13單位保持一致*/程序運行得到第1,2未供水時段的-時間流量圖如下，可以看到與Fig.12.13中用插值給出的曲線比較吻合。如果用5次多項式擬合則得下面的圖形，顯然較3次擬合的效果好。而第1供水時段的流量則用前后時期的流量進行擬合得到。為使流量函數(shù)在t=9和t=11連續(xù)，我們只取4個點，用3次

19、多項式擬合得到第1供水時段的時間-流量圖形如下，可以看到與Fig.12.13中的相應(yīng)部分比較吻合。 dygsdsj= 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(nhjg,8.97:0.01:10.95);gsysl1=100*trapz(gssj1,gs1); */該語句計算第1供水時段的總用水量*/plot(nhsj, 1

20、00*nhlsjg)在第2供水時段之前取t=19.96,20.84兩點的流量，用第三未供水時段的3個記錄做差分得到兩個流量數(shù)據(jù)21.62,18.48，然后用這4個數(shù)據(jù)做3次多項式擬合得到第2供水時段與第3為供水時段的時間-流量圖如下，可以看到與Fig.12.13中的相應(yīng)部分也比較吻合。t3=19.96,20.84,t(22),t(23); ls3=x25*100,21.62,18.48; nhhddxsxs=polyfit(t3,ls3,3); tp3=19.96:0.01:25.91; xx3=polyval(nhhddxsxs, tp3);gssj2=20.84:0.01:24;gs2=p

21、olyval(nhhddxsxs, 20.84:0.01:24);gsysl2=trapz(gssj2,gs2); */該語句計算第2供水時段的總用水量*/ plot(tp3,xx3)3) 一天總用水量的估計分別對供水的兩個時段和不供水的兩個時段積分(流量對時間)并求和得到一天的總用水量約為526.8935(此數(shù)據(jù)是總用水高度cm)。下表列出各段用水量，與插值法算得的表12-7相比，二者較為吻合。表12-9 時段第1未供水段第2未供水段第1供水段第2供水段全天用水用水高度145.67260.6646.6073.9635526.89354） 流量及總用水量的檢驗計算出各時刻的流量可用水位記錄的數(shù)值微分來檢驗，各時段的用水高度可以用實際記錄的水位下降高度來檢驗。例如，算得第1未供水段的用水量高度是145.67，而實際記錄的水位下降高度為968-822=146，兩者是吻合的；同樣地，算得第2未供水段的用水量高度是260.66，而實際記錄的水位下降高度為1082-822=260，兩者也是吻合的。從算法設(shè)計和分析可知,計算結(jié)果與各時段所用的擬合多項式的次數(shù)有關(guān)。下表給出的是對第1,2未供水時段分別用5、6次多