1、數學建模課程設計 報 告課題名稱：_水塔水流量的估計 系 （院）： 理學院 專 業： 數學與應用數學 班 級： 10122111 學生姓名： 錢宇捷 學 號： 1012211134 指導教師： 陳宏宇 開課時間： 2011-2012 學年 二 學期摘 要由所給的題目可知，本問題是一個關于如何計算居民用水的問題，由題目給出的表格，可知不同時刻的水位，根據所要求的不同時刻水位的不同入手，此計算問題就可以轉化為插值或擬合問題。這里主要考慮采用插值的方法，可以利用MATLAB軟件進行插值和曲線擬合計算并解決一些具體的實際問題。根據題目建立模型并采用插值的方法進行求解，推算出任何時刻（包括水泵正供水時）

2、從水塔流出的水流量，及一天的總用水量。 關鍵詞：建模，流量，插值，擬合，MATLABAbstractKnown by the subject to this problem is a problem on how to calculate residential water table given by the subject, we can see the water level of the different moments, according to the water level of the different requirements at different times to

3、start, this calculation problem can be transformed into interpolation or fitting problems. Here consider using the interpolation method, using MATLAB software interpolation and curve fitting calculations and solve some specific practical problems. The topic model and using the interpolation method t

4、o solve the projected out at any time (including when the pump is water supply) to the flow of water flowing from the water tower, and a day of total water use.Key words: Modeling，Flow，Interpolation,Fitting，MATLAB目 錄一、摘要2二、設計目的4三、設計內容4四、問題重述4五、問題分析5六、模型假設5七、模型的建立與求解6八、實驗結果與分析20九、內容小結20十、模型的分析與改進20十一

5、、參考文獻20水塔水流量的估計一、 設計目的1.掌握四種經典的插值方法：拉格朗日插值法、分段插值法、三次樣條插值法。2.學會用MATLAB軟件進行數據差值計算。3.學會用數據插值，數據擬合方法建立數學模型并求解。二、 設計內容1.數據插值、數據擬合理論方法。2.熟悉使用MATLAB軟件進行數據插值，數據擬合。3.簡單的數據建模實驗：水塔水流量的估計。三、 問題重述問題描述某居民區供水機構有一供居民用水的圓柱形水塔，由于該機構沒有測量流入或流出該水塔的水量的裝置，水塔的管理者只能通過隔一段時間測量一次水塔的水位來估計水的流量。但面臨的困難是，當水塔水位下降到設定的最低水位時，水泵自動啟動向水塔供

6、水，到設定的最高水位時停止供水，這段時間無法測量水塔的水位和水泵的供水量。通常每天水泵供水一到二次，每次約2個小時。水塔是一個高為12.2米，直徑為17.4米的正圓柱。按照設計，水塔水位降至約8.2米時，水泵自動啟動，水位升到約10.8米是水泵停止工作。表1是某一天管理人員測量該水塔的水位測量記錄（符號“/”表示水泵啟動）。請您幫助該管理人員估計一天中任何時刻（包括水泵正在供水時刻）從水塔中流出的水流量及一天的居民用水總量。表1 水位測量記錄時間（h）水位（m）時間（h）水位（m）09.67712.95410.210 0.9219.47913.8759.936 1.8439.30814.982

7、9.653 2.9499.125 15.9039.4093.8718.982 16.8269.180 4.9788.814 17.9318.9215.9008.686 19.0378.6627.0068.52519.9598.4337.9288.388 20.8398.2208.9678.22022.015水泵開啟9.9811水泵開啟22.95810.82010.925水泵開啟23.88010.59110.95410.820 24.98610.35412.03210.50025.90810.180水塔的橫截面積為A=237.8 (平方米)。四、問題分析流量是單位時間內流出水的體積,由于水塔是正

8、圓柱形,橫截面積是常數,所以在水泵不工作時段,流量很容易根據水位相對時間的變化率算出。問題的難點在于如何估計水泵供水時段的流量。 水泵供水時段的流量只能靠供水時段前后的流量經插值或擬合得到。作為用于插值或擬合的原始數據,我們希望水泵不工作時段的流量越準確越好。這些流量大體上可由兩種方法計算,一是直接對表1中的水量用數值微分算出各時段的流量,用它們擬合其它時刻或連續時間的流量;二是先用表中數據擬合水位時間函數,求導數即可得到連續時間的流量。有了任何時刻的流量,就不難計算一天的總用水量。其實,水泵不工作時段的用水量可以由測量記錄直接得到,由表1中下降水位乘以水塔的截面積就是這一時段的水流量這個數值

9、可以用來檢驗數據插值或擬合的結果。五、模型假設1. 影響水從水塔中流出的流量的唯一因素是公眾對水的傳統要求因為表給出的數據沒有提及任何其他的影響因素，我們假定所給數據反映了有代表性的一天，而不包括任何特殊情況，如自然災害、火災、水塔溢水、水塔漏水等對水的特殊要求2. 水塔中的水位不影響水流量的大小，氣候條件、溫度變化等也不影響水流量因為物理學的Torricelli定律指出：水塔的最大水流量與水位高度的平方根成正比，按照設計，水塔水位降至約8.2米時，水泵自動啟動，水位升到約10.8米是水泵停止工作。故最低和最高水位分別是8.2m和10.8m (設出口的水位為零)，因為sqrt(10.8/8.2

10、)1.1476，約為1,所以可忽略水位對流速的影響。3. 水泵工作起止時間由它的水位決定，每次充水時間大約為2個小時水泵工作性能效率總是一定的，沒有工作時需維修、使用次數多影響使用效率問題，水泵充水量遠大于水塔水流量.4. 水泵工作時單位時間的供水量大致為常數，這個常數大于單位時間內從水塔中流出的水流的最大流速，這是因為居民區內一直需要用水，不允許水塔中的水用光。5. 水塔中水流量是時間的連續光滑函數，與水泵工作與否無關。這是因為雖然就個別用戶而言可能用水量有較大的變化，但由于個人的用水量與整個居民區用水量相比是非常小的，從統計意義上來講，不太可能同時整個社區的用水量增長或減少。6. 水塔的水

11、流量曲線可以用一條光滑的曲線來逼近這時，水流量曲線的兩階導數是連續的7. 將流量看作是時間的連續函數，為計算簡單,不妨將流量定義成單位時間流出水的高度,即水位對時間變化率的絕對值(水位是下降的), 水塔截面積為s=237.8 (平方米)得到結果后乘以S即可。 六、模型的建立與求解通過以上對問題的分析，現在的問題已轉化為根據某一天已測量的時刻水塔中水的流速，產生在整個區間（24小時）上的函數或函數值，一般來說插值和擬合是兩種最常用的方法。1）流量估計方法 首先依照表1所給數據,用MATLAB作出時間水位散點圖1。Matlab程序：>> x0=0 0.921 1.843 2.949 3

12、.871 4.978 5.900 7.006 7.928 8.967 10.954 12.032 12.954 13.875 14.982 15.903 16.826 17.931 19.037 19.959 20.839 22.958 23.880 24.986 25.908;>> y0=9.677 9.479 9.308 9.125 8.982 8.814 8.686 8.525 8.388 8.220 10.820 10.500 10.210 9.936 9.653 9.409 9.180 8.921 8.662 8.433 8.220 10.820 10.591 10.35

13、4 10.180;>> plot(x0,y0,'*')圖1下面我們來計算水箱流量與時間的關系。根據表1中數據散點圖1，一種簡單的處理方法為先將表1中的數據分為三段，然后對每一段的數據做如下處理：設某段數據為，相鄰數據中點的平均流速用下面的公式(流速=（左端點的水位-右端點的水位）/區間長度)：，每段數據首尾點的流速用下面的公式計算：，用以上公式求得時間與流速之間的數據表2如下: 表2時間(h)流速(cm/h)時間(h)流速(cm/h)022.951712.49331.45340.460521.498413.414529.75031.38218.546614.4285

14、25.56462.39616.546115.442526.49293.4115.509816.364524.81044.424515.176217.378523.43895.43913.882918.48423.41776.45314.557019.49824.83737.46714.859020.39924.20458.447516.169420.83922.75258.96718.714922.015水泵開動9.9811水泵開動22.95822.189310.925水泵開動23.41924.837310.95433.5024.43321.428611.49329.684625.44718.

15、872025.90814.0533由表2作出圖2時間-流速散點圖如下:Matlab程序：>> x0=0 0.4605 1.382 2.396 3.41 4.4245 5.439 6.453 7.467 8.4475 8.967 10.954 11.493 12.493 13.4145 14.4285 15.4425 16.3645 17.3785 18.484 19.498 20.399 20.839 22.958 23.419 24.433 25.447 25.908;>> y0=22.9517 21.4984 18.5466 16.5461 15.5098 15.1

16、762 13.8829 14.5570 14.8590 16.1694 18.7149 33.50 29.6846 31.4534 29.7503 25.5646 26.4929 24.8104 23.4389 23.4177 24.8373 24.2045 22.7525 22.1893 24.8373 21.4286 18.8720 14.0533;>> plot(x0,y0,'*')圖2(1) 插值法由表2，對水泵不工作時段1,2采取插值方法,可以得到任意時刻的流速,從而可以知道任意時刻的流量.我們分別采取拉格朗日插值法,分段線性插值法及三次樣條插值法；對于水

17、泵工作時段1應用前后時期的流速進行插值，由于最后一段水泵不工作時段數據太少，我們將它與水泵工作時段2合并一同進行插值處理(該段簡稱混合時段)。我們總共需要對四段數據（第1，2未供水時段，第1供水時段，混合時段）進行插值處理。第1未供水時段：首先，對于第1未供水時段，我們分別用三種方法算出流量函數和用水量（用水高度）。下面為實現該過程的MATLAB程序：>> t=0 0.4605 1.382 2.396 3.41 4.4245 5.439 6.453 7.467 8.4475 8.967;>> v=22.9517 21.4984 18.5466 16.5461 15.50

18、98 15.1762 13.8829 14.5570 14.8590 16.1694 18.7149;>> t0=0:0.1:8.967;>> lglr=lglrcz(t,v,t0); >> 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,'

19、*',t0,lglr,'r',t0,fdxx,'g',t0,scyt,'b')>> gtext('lglr')>> gtext('fdxx')>> gtext('scyt')運行結果為：lglrjf = 144.7505fdxxjf =145.4062sancytjf = 144.8002圖3是對第1未供水段數據用三種不同方法得到的插值函數圖，圖中曲線lglr、fdxx和scyt分別表示用拉格朗日插值法,分段線性插值法及三次樣條插值法得到的曲線。圖3分析：

20、由表1知，第1未供水時段的總用水高度為145.7（=967.7-822.0），可見上述三種插值方法計算的結果與實際值（145.7）相比都比較接近。但是考慮到分段線性插值方法具有更加良好的性質，建議采取該方法。第2未供水時段：首先，對于第2未供水時段，我們分別用三種方法算出流量函數和用水量（用水高度）。下面為實現該過程的MATLAB程序：>> t=10.954 11.493 12.493 13.4145 14.4285 15.4425 16.3645 17.3785 18.484 19.498 20.399 20.839;>> v=33.50 29.6846 31.453

21、4 29.7503 25.5646 26.4929 24.8104 23.4389 23.4177 24.8373 24.2045 22.7525;>> t0=10.954:0.1:20.839;>> lglr=lglrcz(t,v,t0); >> 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

22、(scyt)>> plot(t,v,'*',t0,lglr,'r',t0,fdxx,'g',t0,scyt,'b')>> gtext('lglr')>> gtext('fdxx')>> gtext('scyt')運行結果為：lglrjf =258.9555fdxxjf =259.1952sancytjf =258.7886圖4是對第2未供水段數據用三種不同方法得到的插值函數圖，圖中曲線lglr、fdxx和scyt分別表示用拉格朗日插值法

23、,分段線性插值法及三次樣條插值法得到的曲線。圖4分析：由表1可知，第2未供水時段的總用水高度為260（260=1082-822），可見上述三種插值方法計算的結果與實際值260相比都比較接近。但是考慮到分段線性插值方法更加接近實際值，所以建議采取該方法。第1供水時段：首先，對于第1供水時段，我們分別用三種方法算出流量函數和用水量（用水高度）。下面為實現該過程的MATLAB程序：>> t=7.467 8.4475 8.967 10.954 11.493 12.493;>> v=14.8590 16.1694 18.7149 33.50 29.6846 31.4534;>

24、;> t0=8.967:0.1:10.954;>> lglr=lglrcz(t,v,t0); >> 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

25、,scyt,'b')>> gtext('lglr')>> gtext('fdxx')>> gtext('scyt')運行結果為：lglrjf =53.2376fdxxjf =48.9892sancytjf =52.5114圖5是對第1供水段數據用三種不同方法得到的插值函數圖，圖中曲線lglr、fdxx和scyt分別表示用拉格朗日插值法,分段線性插值法及三次樣條插值法得到的曲線。圖5混合時段：首先，對于混合時段，我們分別用三種方法算出流量函數和用水量（用水高度）。下面為實現該過程的MATLAB程序

26、：>>t=19.498 20.399 20.839 22.958 23.419 24.433 25.447 25.908;>> v=24.8373 24.2045 22.7525 22.1893 24.8373 21.4286 18.8720 14.0533;>> t0=20.839:0.1:22.958;>> lglr=lglrcz(t,v,t0); >> lglrjf=0.1*trapz(lglr)>> fdxx=interp1(t,v,t0);>> fdxxjf=0.1*trapz(fdxx)>&g

27、t; 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')>> gtext('fdxx')>> gtext('scyt')運行結果為：lglrjf =39.0871fdxxjf =47.1942sancytjf =

28、43.5828圖6是對混合時段數據用三種不同插值方法得到的插值函數圖。圖6數值結果整理成表格如下：各時段及一天的總用水量(用水高度)第1未供水段第2未供水段第1供水段混合時段全天拉格朗日插值法144.7505258.955553.237639.0871496.0307分段線性插值法145.4062259.195248.989247.1942500.7848三次樣條插值法144.8002258.788652.511443.5828499.6830表三作出分段線性及三次樣條插值方法得到的整個過程的時間-流速函數示意圖圖7。Matlab程序如下：>>t=0 0.4605 1.382 2.

29、396 3.41 4.4245 5.439 6.453 7.467 8.4475 8.967 10.954 11.493 12.493 13.4145 14.4285 15.4425 16.3645 17.3785 18.484 19.498 20.399 20.839 22.958 23.419 24.433 25.447 25.908;>>v=22.9517 21.4984 18.5466 16.5461 15.5098 15.1762 13.8829 14.5570 14.8590 16.1694 18.7149 33.50 29.6846 31.4534 29.7503 2

30、5.5646 26.4929 24.8104 23.4389 23.4177 24.8373 24.2045 22.7525 22.1893 24.8373 21.4286 18.8720 14.0533;>>t0=0:0.1:25.908;>>lglr=lglrcz(t,v,t0); >>lglrjf=0.1*trapz(lglr)>>fdxx=interp1(t,v,t0);>>fdxxjf=0.1*trapz(fdxx)>>scyt=interp1(t,v,t0,'spline');>>s

31、ancytjf=0.1*trapz(scyt);>>plot(t,v,'*' ,t0,fdxx,'g',t0,scyt,'b')>>gtext('lglr')>>gtext('fdxx')>>gtext('scyt')運行程序，得到分段線性及三次樣條插值函數圖圖7如下：圖7下表是對一天中任取的4個時刻分別用三種方法得到的水塔水流量近似值。6.8510.8515.8522.85拉格朗日15.0233.7726.0721.23分段插值14.6633.722

32、5.7422.22三次樣條14.7432.8526.0321.52表四（2）擬合法1) 擬合水位-時間函數從表1中測量記錄看，一天有兩次供水時段和三次未供水時段，分別對第1,2未供水時段的測量數據直接作多項式擬合，可得到水位函數（注意，根據多項式擬合的特點，此處擬合多項式的次數不宜過高，一般以3-6次為宜）。對第3未供水時段來說,數據過少不能得到很好的擬和。設t,h分別為已輸入的時刻和水位測量記錄(由表1提供，水泵啟動的4個時刻不輸入)（a） 對第1未供水時段的數據進行擬合。Matlab程序如下：>> t=0 0.921 1.843 2.949 3.871 4.978 5.900

33、7.006 7.928 8.967 10.954 12.032 12.954 13.875 14.982 15.903 16.826 17.931 19.037 19.959 20.839 22.958 23.880 24.986 25.908;>> h=9.677 9.479 9.308 9.125 8.982 8.814 8.686 8.525 8.388 8.220 10.820 10.500 10.210 9.936 9.653 9.409 9.180 8.921 8.662 8.433 8.220 10.820 10.591 10.354 10.180;>>

34、c1=polyfit(t(1:10),h(1:10),3);>> tp1=0:0.1:8.9;>> x1=polyval(c1,tp1);>> plot(tp1,x1)圖8 （b）對第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)圖92) 確定流量時間函數 （a）第1，2未供水時間段對于第1,2未供水時段的流量可直接對水位函數求

35、導Matlab程序如下：>>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.967;>>tp2=10.954:0.01:20.839;>>x13=-polyval(a1,tp1);>>x113=-polyval(a1,0:0.01: 8.967);>>wgsysl1=100*trapz(tp1,x113); */

36、該語句計算第1未供水時段的總用水量*/>>x14=-polyval(a1, 7.928, 8.967);*/該語句僅為下面的程序準備數據*/>>x23=-polyval(a2,tp2);>>x114=-polyval(a2, 10.954:0.01: 20.839);>>wgsysl2=100*trapz(tp2,x114); */該語句計算第2未供水時段的總用水量*/>>x24=-polyval(a2, 10.954, 12.032);*/該語句僅為下面的程序準備數據*/>>x25=-polyval(a2, 19.959

37、,20.839); */該語句僅為下面的程序準備數據*/>>subplot(1,2,1)>>plot(tp1,x13*100) */與圖（7）單位保持一致*/>>subplot(1,2,2)>>plot(tp2,x23*100) */與圖（7）單位保持一致*/程序運行，得出用3次多項式擬合的第1,2未供水時段時間-流量圖圖10如下：圖10如果用5次多項式擬合則得下面的圖形。Matlab程序： >> c1=polyfit(t(1:10),h(1:10),5); */與3次多項式的區別*/>> c2=polyfit(t(11:

38、21),h(11:21),5); */與3次多項式的區別*/>> a1=polyder(c1);>> a2=polyder(c2);>> tp1=0:0.01:8.967;>> tp2=10.954:0.01:20.839;>> x13=-polyval(a1,tp1);>> x113=-polyval(a1,0:0.01:8.967);>> wgsysl1=100*trapz(tp1,x113); >> x14=-polyval(a1,7.928,8.967); >> x23=-pol

39、yval(a2,tp2);>> x114=-polyval(a2,10.954:0.01:20.839);>> wgsysl2=100*trapz(tp2,x114); >> x24=-polyval(a2,10.954,12.032); >> x25=-polyval(a2,19.959,20.839); >> subplot(1,2,1)>> plot(tp1,x13*100) >> subplot(1,2,2)>> plot(tp2,x23*100)程序運行，得出用5次多項式擬合的第1,2未供

40、水時段時間-流量圖圖11如下：圖11分析： 相比較而言，對于第1，2未供水時段，5次多項式比3次多項式擬合的效果更好，更佳。故建議用5次多項式進行擬合（b）第1供水時間段第1供水時段的流量則用前后時期的流量進行擬合得到。為使流量函數在t=9和t=11連續，我們只取4個點，用3次多項式擬合得到第1供水時段的時間-流量圖形如下，可以看到與圖7中的相應部分比較吻合。Matlab程序如下：>>dygsdsj= 7.928,8.967, 10.954,12.032;>> dygsdls=x14, x24;>> nhjg=polyfit(dygsdsj, dygsdls

41、,3);>> nhsj=7.928:0.1:12.032;>> nhlsjg=polyval(nhjg ,nhsj);>> gssj1=8.967:0.01:10.954;>> gs1=polyval(nhjg,8.967:0.01:10.954);>> gsysl1=100*trapz(gssj1,gs1) */該語句計算第1供水時段的總用水量*/>>plot(nhsj, 100*nhlsjg)圖12用5次多項式擬合得到第1供水時段的時間-流量圖形如下，可以看到與圖7中的相應部分比較吻合。Matlab程序如下：>&

42、gt;dygsdsj= 7.928,8.967, 10.954,12.032;>> dygsdls=x14, x24;>> nhjg=polyfit(dygsdsj, dygsdls,5); */與3次多項式的區別*/>> nhsj=7.928:0.1:12.032;>> nhlsjg=polyval(nhjg ,nhsj);>> gssj1=8.967:0.01:10.954;>> gs1=polyval(nhjg,8.967:0.01:10.954);>> gsysl1=100*trapz(gssj1,gs

43、1) */該語句計算第1供水時段的總用水量*/>>plot(nhsj, 100*nhlsjg)圖13分析： 相比較而言，對于第1供水時段，3次多項式比5次多項式擬合的效果更好，更佳。故建議用3次多項式進行擬合（c）混合時段在第2供水時段之前取t=19.037，19.959，20.839三點的流量，用第三未供水時段的4個記錄做差分得到三個流量數據24.8373，21.4286，18.8720，然后用這6個數據做3次多項式擬合得到第2供水時段與第3未供水時段即混合時段的時間-流量圖圖14Matlab程序如下：>>t3=19.037,19.959,20.839,t(22),t

44、(23); >>ls3=x25*100,24.8373,21.4286,18.8720; >> nhhddxsxs=polyfit(t3,ls3,3); >>tp3=19.96:0.01:25.91; >> xx3=polyval(nhhddxsxs, tp3);>>gssj2=20.839:0.01:23; >>gs2=polyval(nhhddxsxs, 20.839:0.01:23);>>gsysl2=trapz(gssj2,gs2); */該語句計算第2供水時段的總用水量*/ >>plot(

45、tp3,xx3)用上述6個數據做3次多項式擬合得到第2供水時段與第3未供水時段即混合時段的時間-流量圖如下，圖145次多項式擬合Matlab程序如下：>>t3=19.037,19.959,20.839,t(22),t(23);>> ls3=x25*100,24.8373,21.4286,18.8720;>> nhhddxsxs=polyfit(t3,ls3,5); */與混合時段3次多項式擬合的區別/*>> tp3=19.96:0.01:25.91;>> xx3=polyval(nhhddxsxs, tp3);>>gssj2=20.84:0.01:23;>>gs2=polyval(nhhddxsxs, 20.839:0.01:23);>>gsysl2=trapz(gssj2,gs2);>>plot(tp3,xx3)用上述6個數據做5次多項式擬合得到第2供水時段與第3未供水時段即混合時段的時間-流量圖如下，圖153) 一天總用水量的估計分別對供水的兩個時段和不供水的兩個時段積分(流量對時間)并求和得到一天的總用水量約