居民區(qū)供水問(wèn)題

1、居民區(qū)供水問(wèn)題的研究與分析摘要：本文針對(duì)居民區(qū)供水問(wèn)題進(jìn)行研究與討論。根據(jù)擬合原理對(duì)水塔的流量進(jìn)行模擬，運(yùn)用數(shù)值積分對(duì)一天的總用水量進(jìn)行估計(jì)。本文建立了水塔流量模擬的數(shù)學(xué)模型，并對(duì)模型進(jìn)行了分析與評(píng)價(jià)。 對(duì)于問(wèn)題一，流量是單位時(shí)間流出的水的體積，由于水塔是正圓柱形，橫截面積是常數(shù)，在水泵不工作的時(shí)段，流量很容易從水位對(duì)時(shí)間的變化率算出，水泵供水時(shí)段的流量只能靠供水時(shí)段前后的流量擬合得到，在本文的計(jì)算中，我們先表中所給數(shù)據(jù)擬合水位時(shí)間函數(shù)，對(duì)此函數(shù)進(jìn)行求倒，得到連續(xù)時(shí)間的流量函數(shù)。對(duì)于問(wèn)題二，運(yùn)用問(wèn)題一得出的結(jié)果，分別對(duì)第1、2時(shí)段和第1、2供水時(shí)段流量進(jìn)行積分，然后求和，得出的結(jié)果即是一天的總

2、用水量。對(duì)于問(wèn)題三,根據(jù)水泵的功率（單位時(shí)間泵如的水量）等于供水時(shí)段的用水量加上水位上升值除以時(shí)段長(zhǎng)度，可以很容易的得出水泵的功率。關(guān)鍵詞：水塔 流量 擬合 模擬 數(shù)值積分一、問(wèn)題的提出    某居民區(qū)的民用自來(lái)水是由圓柱形水塔提供,水塔高12.2米,直徑17.4米.水塔是由水泵根據(jù)水塔內(nèi)水位高低自動(dòng)加水,一般每天水泵工作兩次.現(xiàn)在需要了解居民區(qū)用水規(guī)律與水泵的工作功率.按照設(shè)計(jì),當(dāng)水塔的水位降至最低水位,約8.2米,水泵自動(dòng)啟動(dòng)加水;當(dāng)水位升高到一個(gè)最高水位, 約10.8米,水泵停止工作.    可以考慮采用用水率(單位時(shí)

3、間的用水量)來(lái)反映用水規(guī)律,并通過(guò)間隔一段時(shí)間測(cè)量水塔里的水位來(lái)估算用水率,表1是某一天的測(cè)量記錄數(shù)據(jù),測(cè)量了28個(gè)時(shí)刻,但是由于其中有3個(gè)時(shí)刻遇到水泵正在向水塔供水,而無(wú)水位記錄(表1中用/表示).    試建立合適的數(shù)學(xué)模型,推算任意時(shí)刻的用水率,一天的總用水量和水泵工作功率.    表1 原始數(shù)據(jù)(單位:時(shí)刻(小時(shí)),水塔中水位(米)時(shí)刻t00.9211.8432.9493.8714.9785.900水位9.6779.4799.3089.1258.9828.8148.686時(shí)刻t7.0067.9288.9679.981110.9

4、2510.95412.032水位8.5258.3888.220/10.82010.500時(shí)刻t12.95413.87514.98215.90316.82617.93119.037水位10.2109.9369.6539.4099.1808.9218.662時(shí)刻t19.95920.83922.01522.95823.88024.98625.908水位8.4338.220/10.82010.59710.35410.180二、模型的基本假設(shè)1、 流量只取決于水位差，與水位本身無(wú)關(guān)，由于從小孔流出的流體的流速正比于水面高度的平房根，由題目給出的水塔的最低和最高水位分別為8.2m和10.8m（設(shè)出口的水位

5、為0），所以可以忽略水位對(duì)流速的影響。因?yàn)?、 水泵第一次供水時(shí)段為到，第二次供水時(shí)段為到，這是根據(jù)最低和最高水位分別是8.2m和10.8m，及表一的水位測(cè)量記錄做出的假設(shè)，其中前3個(gè)時(shí)刻直接取自實(shí)測(cè)數(shù)據(jù)（精確到0.1h），最后一個(gè)時(shí)刻來(lái)自每次供水約兩個(gè)小時(shí)的已知條件（從記錄看，大約在20.84h之后開(kāi)始供水，大約在22.9h之前結(jié)束供水）。3、 假設(shè)水泵工作時(shí)單位時(shí)間的供水量是常數(shù)，此常數(shù)大于單位時(shí)間的平均流量。4、 假設(shè)流量是關(guān)于時(shí)間的連續(xù)函數(shù)。5、 假設(shè)流量與水泵是否工作無(wú)關(guān)。6、 由于水塔截面積是常數(shù)，計(jì)算中將流量定義為單位時(shí)間流出的水的高度，即水位對(duì)時(shí)間變化率的絕對(duì)值，將最后的乘以S

6、即可。即：水位是時(shí)間的連續(xù)函數(shù) 水位對(duì)時(shí)間的變化率（流量） 任何時(shí)刻的流量： 三：符號(hào)說(shuō)明：測(cè)量時(shí)刻：水塔截面積：水位測(cè)量記錄：第1時(shí)段的用水量：第2時(shí)段的用水量：第1供水時(shí)段用水量：第2供水時(shí)段用水量：一天總用水量：第1時(shí)段水位下降高度：第2時(shí)段水位下降高度：第1供水時(shí)段水位下降高度：第2供水時(shí)段和第3時(shí)段水位下降高度四：?jiǎn)栴}的分析與準(zhǔn)備1、 擬合水位時(shí)間函數(shù)從表1測(cè)量記錄看，一天有兩個(gè)供水時(shí)段（簡(jiǎn)稱(chēng)第1供水時(shí)段和第2供水時(shí)段），和三個(gè)水泵不工作時(shí)間段（以下簡(jiǎn)稱(chēng)第1時(shí)間段t=0到t=8.967，第2時(shí)間段t=10.954到t=20.839和第3時(shí)段t=23以后）。對(duì)第1，2時(shí)段的測(cè)量數(shù)據(jù)直接

7、分別作多項(xiàng)式擬合，得到水位函數(shù)。為使擬合曲線(xiàn)比較光滑，多項(xiàng)式次數(shù)不宜太高，一般在36.由于第三時(shí)段只有三個(gè)測(cè)量記錄，無(wú)法對(duì)這一時(shí)段的水位做出較好的擬合。2、 確定流量時(shí)間函數(shù)對(duì)于第1、2時(shí)段只要將水位函數(shù)求導(dǎo)數(shù)即可，對(duì)于兩個(gè)供水時(shí)段的流量，則用供水時(shí)段前后（水泵不工作時(shí)段）的流量擬合得到，并且將擬合得到的第二供水時(shí)段流量外推，將第3時(shí)段流量包含在第2供水時(shí)段內(nèi)，需要擬合四個(gè)流量函數(shù)。3、 一天總用水量的估計(jì)總用水量等于兩個(gè)水泵不工作時(shí)段和兩個(gè)供水時(shí)段用水量之和，他們都可以由流量對(duì)時(shí)間的積分得到：五：模型的建立與求解1、擬合第1、2時(shí)段的水位，并導(dǎo)出流量 設(shè)t，h為已輸入的時(shí)刻和水位測(cè)量記錄（水

8、泵啟動(dòng)的4個(gè)時(shí)刻不輸入），為了計(jì)算方便，也為了提高精度，將水位的單位化為cm，時(shí)間精確到0.01h。t=0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 0 0 1082 1050 1021 994 965 941 918 892 866 843 822 0

9、 0 1059 1035 1018;第1時(shí)段各時(shí)刻的流量可如下得：1）c1=polyfit（t（1：10），h（1：10），3）； %用3次多項(xiàng)式擬合第1時(shí)段水位，c1輸出3次多項(xiàng)式的系數(shù)2）a1=polyder（c1）；% a1輸出多項(xiàng)式（系數(shù)為c1）導(dǎo)數(shù)的系數(shù)3）tp1=0：0.1：9；x1=-polyval（a1，tp1）；% x1輸出多項(xiàng)式（系數(shù)a1）在tp1點(diǎn)的函數(shù)值（取負(fù)后邊為正值），即tp1時(shí)刻的流量。4）流量函數(shù)為：具體程序見(jiàn)附錄shuita1.m用程序進(jìn)行計(jì)算，求得在0，9內(nèi)各時(shí)刻的流量值(水位變化率)如下表：ti00.921.842.953.874.985.907.017.

10、938.97h122.1119.8117.9116.1415.1214.4214.2814.6415.3816.69fig.1 第一用水時(shí)段流量曲線(xiàn)圖第2時(shí)段各時(shí)刻的流量可如下得： 1）c2=polyfit(t(13:23),h(13:23),4) %用4次多項(xiàng)式擬合第2時(shí)段水位，c2輸出4次多項(xiàng)式的系數(shù)2）a2=polyder(c2) % a2輸出多項(xiàng)式（系數(shù)為c2）導(dǎo)數(shù)的系數(shù)3）tp2=10.9:0.1:21 x2=-polyval(a2,tp2) % x2輸出多項(xiàng)式（系數(shù)為a2）在tp2點(diǎn)的函數(shù)值（取負(fù)后邊為正值），即tp2時(shí)刻的流量4）流量函數(shù)為：具體程序見(jiàn)附錄shuita2.m用程序

11、進(jìn)行計(jì)算，求得在11,20.8內(nèi)各時(shí)刻的流量值(水位變化率)如下表：ti10.9512.0312.9513.8814.9815.9016.8317.9319.0419.9620.84h131.8030.5329.2927.9926.5025.4024.5223.9123.9024.4625.56fig.2 第二用水時(shí)段流量曲線(xiàn)圖2. 擬合供水時(shí)段的流量 在第1供水時(shí)段（t=911）之前（即第1時(shí)段）和之后（即第2時(shí)段）各取幾點(diǎn)，其流量已經(jīng)得到，用它們擬合第1供水時(shí)段的流量為使流量函數(shù)在t=9和t=11連續(xù)，我們簡(jiǎn)單地只取4個(gè)點(diǎn)，擬合3次多項(xiàng)式（即曲線(xiàn)必過(guò)這4個(gè)點(diǎn)），實(shí)現(xiàn)如下：xx1=-poly

12、val(a1,8 9);%取第1時(shí)段在t=8,9的流量 xx2=-polyval(a2,11 12);%取第2時(shí)段在t=11,12的流量 xx12=xx1 xx2; c12=polyfit(8 9 11 12,xx12，3);%擬合3次多項(xiàng)式 tp12=9：0.1：11; x12=polyval(c12,tp12);%x12輸出第1供水時(shí)段各時(shí)刻的流量 擬合的流量函數(shù)為： 具體程序見(jiàn)附錄shuita3.m用程序進(jìn)行計(jì)算，求得在9,11內(nèi)各時(shí)刻的流量值(水位變化率)如下表：ti911h116.73831.72 fig.3 第一供水時(shí)段流量曲線(xiàn)圖在第2供水時(shí)段之前取t=20，20.8兩點(diǎn)的流水量，

13、在該時(shí)刻之后（第3時(shí)段）僅有3個(gè)水位記錄，我們用差分得到流量，然后用這4個(gè)數(shù)值擬合第2供水時(shí)段的流量如下： dt3=diff（t(22：24)）； %最后3個(gè)時(shí)刻的兩兩之差 dh3=diff（h(22：24)）； %最后3個(gè)水位的兩兩之差dht3=-dh3./dt3； %t(22)和t(23)的流量t3=20 20.8 t(22) t(23)； xx3=-polyval(a2，t3(1：2)，dht3)； %取t3各時(shí)刻的流量 c3=polyfit（t3，xx3，3）； %擬合3次多項(xiàng)式 t3=20.8：0.1：24； x3=polyval（c3，tp3）； % x3輸出第2供水時(shí)段（外推至t

14、=24）各時(shí)刻的流量擬合的流量函數(shù)為： 具體程序見(jiàn)shuita4.m 用程序進(jìn)行計(jì)算，求得在20.8,24內(nèi)各時(shí)刻的流量值(水位變化率)如下表： ti22.0122.9623.8824.9925.91h125.57423.93721.62218.478 16.08 fig.4 第2供水時(shí)段和第3用水時(shí)段流量曲線(xiàn)圖3. 一天總用水量的估計(jì)第1、2時(shí)段和第1、2供水時(shí)段流量的積分之和，就是一天總用水量雖然諸時(shí)段的流量已表為多項(xiàng)式函數(shù)，積分可以解析地算出，這里仍用數(shù)值積分計(jì)算：（1）第一用水時(shí)段的用水量其中積分值通過(guò)梯形公式計(jì)算 （2）第二用水時(shí)段的用水量 （3）第一供水時(shí)段的用水量 (4). 第二

15、供水和第三用水時(shí)段的用水量 y1=0.1*trapz(x1)； %第1時(shí)段用水量（仍按高度計(jì)），0.1為積分步長(zhǎng)y2=0.1*trapz(x2)； %第2時(shí)段用水量y12=0.1*trapz(x12)； %第1供水時(shí)段用水量y3=0.1*trapz(x3)； %第2供水時(shí)段用水量y=(y1+y2+y12+y3)*237.8*0.01; %一天總用水量計(jì)算結(jié)果：y1=146.18, y2=266.59, y12=50.25, y3=70.00，y=1267.52具體程序見(jiàn)shuita5.m通過(guò)程序可以計(jì)算出各時(shí)段的用水量及一天的總用水量如下：各時(shí)段用水量第一用水時(shí)段用量第二用水時(shí)段用量第一供水時(shí)

16、段用量第二供水時(shí)段用量一天總用水量ViV1V2V12V3V347.62633.95119.49166.461267.52其中四個(gè)時(shí)段都用三次多項(xiàng)式進(jìn)行擬合，如果第一用水時(shí)段、第一、二供水時(shí)段用三次多項(xiàng)式，第二用水時(shí)段用四次多項(xiàng)式擬合，得另一計(jì)算結(jié)果：y1=146.18 y2=266.75 y12=49.00 y3= 77.34 程序見(jiàn)shuita6.m各時(shí)段的用水量及一天的總用水量如下:各時(shí)段用水量第一用水時(shí)段用量第二用水時(shí)段用量第一供水時(shí)段用量第二供水時(shí)段用量一天總用水量ViV1V2V12V3V347.62613.76115.33187.021263.7 fig.5一天內(nèi)用水流量曲線(xiàn)圖4.

17、流量及總用水量的檢驗(yàn)與功率的計(jì)算 計(jì)算出的各時(shí)刻的流量可用水位記錄的數(shù)值微分來(lái)檢驗(yàn)用水量可用第1時(shí)段水位測(cè)量記錄中下降高度來(lái)計(jì)算并檢驗(yàn)，在第一用水時(shí)段水的實(shí)際用量為 類(lèi)似地，用來(lái)計(jì)算和檢驗(yàn)，在第二用水時(shí)段實(shí)際用水量為：而計(jì)算值為：=347.62, =613.767比較接近。供水時(shí)段流量的一種檢驗(yàn)方法如下：供水時(shí)段的用水量加上水位上升值260是該時(shí)段泵入的水量，除以時(shí)段長(zhǎng)度得到水泵的功率（單位時(shí)間泵入的水量），而兩個(gè)供水時(shí)段水泵的功率應(yīng)大致相等第1、2時(shí)段水泵的功率可計(jì)算如下：p1=(y12+260)/2； %第1供水時(shí)段水泵的功率（水量仍以高度計(jì)） tp2=20.8：0.1：23；xp2=po

18、lyval（c3，tp4）； % xp2輸出第2供水時(shí)段各時(shí)刻的流量p2=(0.1*trapz(xp2)+260)/2.2；%第2供水時(shí)段水泵的功率（水量仍以高度計(jì)）計(jì)算結(jié)果：p1=153.7 ，p2=140.1具體程序見(jiàn)附錄shuita7.m.六：模型的分析與改進(jìn)從計(jì)算結(jié)果可以看出，第1時(shí)段用水量與水位測(cè)量記錄中的下降高度146相差無(wú)幾，第2時(shí)段用水量與記錄中的下降高度260相差無(wú)幾，所以數(shù)據(jù)擬合、數(shù)值積分的精度是足夠的。對(duì)不同次數(shù)的擬合多項(xiàng)式，第1、2供水時(shí)段用水量相差稍大，兩供水時(shí)段水泵的功率也有差別，這都說(shuō)明供水時(shí)段用3次曲線(xiàn)通過(guò)四點(diǎn)的做法不夠好，應(yīng)該多取幾點(diǎn)做擬合，但要注意讓流量曲線(xiàn)

19、在不同的時(shí)段相接處保持連續(xù)。下面給出擬合第1、2用水時(shí)段水位函數(shù)，采用不同次數(shù)的多項(xiàng)式所得流量及總用水量的結(jié)果，用n = (n1,n2)表示這兩個(gè)時(shí)段所用多項(xiàng)式的次數(shù)，下表是各時(shí)段的用水量、一天總用水量及兩個(gè)供水時(shí)段水泵的功率。(n1,n2)V1V2V12V3Vp1p2(3,4)347.62613.76115.33187.021263.7366.81366.05(5,6)348.41612.95108.79181.761251.9363.53363.66兩個(gè)供水時(shí)段水泵的功率相當(dāng)，說(shuō)明該算法切實(shí)可行！另外可以得知各個(gè)時(shí)段的用水量及一天的總用水量，便于掌握用水情況。 fig.6 n=(3,4)流

20、量曲線(xiàn) fig.7 n=(5,6)流量曲線(xiàn)由上圖可以看出，流量曲線(xiàn)與原始記錄基本上那個(gè)相吻合，0點(diǎn)到10點(diǎn)流量很低，10點(diǎn)到下午3點(diǎn)是用水高峰，全天流量平均在22（）左右，若按這個(gè)平均流量計(jì)算，一天的總用水量應(yīng)為與表二的結(jié)果很接近。七：模型的推廣與應(yīng)用通過(guò)對(duì)水塔流量的問(wèn)題進(jìn)行研究，掌握其工作和控制原理并對(duì)模型進(jìn)行評(píng)價(jià)，此模型不僅適用于水塔的供水水位控制系統(tǒng)，而且還可以應(yīng)用于其他的領(lǐng)域，像自動(dòng)蓄水裝置、水壩閘門(mén)控制蓄水系統(tǒng)。通過(guò)數(shù)據(jù)統(tǒng)計(jì)，建立模型，掌握水量的需求情況，合理控制閘門(mén)，使水量保持在一定的范圍內(nèi)，從而達(dá)到最優(yōu)的效果。八：參考文獻(xiàn)1 趙靜 但琦，數(shù)學(xué)建模與數(shù)學(xué)實(shí)驗(yàn)，北京：高等教育出版社，

21、2000.112 陳杰，Matlab寶典，北京：電子工業(yè)出版社，2007.13 朱仁峰，精通7，北京清華大學(xué)學(xué)研大廈：北京大學(xué)出版社，2006.5九：附錄程序Shuita1.mt=0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 0 0 1082 1050

22、1021 994 965 941 918 892 866 843 822 0 0 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3)a1=polyder(c1)tp1=0:0.1:9x1=-polyval(a1,tp1)plot(tp1,x1,'r.')程序shuita2.mt=0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.96 20.84 22.01 22.

23、96 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 0 0 1082 1050 1021 994 965 941 918 892 866 843 822 0 0 1059 1035 1018;c2=polyfit(t(13:23),h(13:23),4)a2=polyder(c2)tp2=10.9:0.1:21;x2=-polyval(a2,tp2);plot(tp2,x2,'b.')hold onxlabel('hour');ylabel('cm/hour');程序Shu

24、ita3.mt=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3)a1=polyder(c1);c2=polyfit(t(10

25、.9:21),h(10.9:21),4)a2=polyder(c2);xx1=-polyval(a1,8 9);xx2=-polyval(a2,11 12);xx12=xx1 xx2;c12=polyfit(8 9 11 12,xx12,3)a12=polyder(c12)tp12=9:0.02:11;x12=polyval(c12,tp12)plot(tp12,x12,'r.')hold onaxis(9 11 16 32)xlabel('hour')ylabel('cm/houor')程序shuita4.mt=0 0.92 1.84 2.95

26、 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c2=polyfit(t(11:21),h(11:21),4);a2=polyder(c2);dt3=diff(t(22:24);dh3=diff(h(22:24);

27、dht3=-dh3./dt3;t3=20 20.8 t(22) t(23)xx3=-polyval(a2,t3(1:2),dht3;c3=-polyfit(t3,xx3,3);a3=polyder(c3)tp3=20.8:0.05:24;x3=-polyval(c3,tp3);plot(tp3,x3,'b.')hold on axis(20.5 24.5 18 27)xlabel('hour')ylabel('cm/houor')title('流量曲線(xiàn)圖')程序shuita5.mt=0 0.92 1.84 2.95 3.87 4.

28、98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 931 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3);a1=polyder(c1);tp1=0:0.1:9;x1=-polyval(a1,tp1);c2=polyfit(t(1

29、1:21),h(11:21),3);a2=polyder(c2);tp2=10.9:0.1:21;x2=-polyval(a2,tp2);xx1=-polyval(a1,8 9);xx2=-polyval(a2,11 12);xx12=xx1 xx2;c12=polyfit(8 9 11 12,xx12,3);tp12=9:0.1:11;x12=polyval(c12,tp12);dt3=diff(t(22:24);dh3=diff(h(22:24);dht3=-dh3./dt3;t3=20 20.8 t(22) t(23);xx3=-polyval(a2,t3(1:2),dht3;c3=-p

30、olyfit(t3,xx3,3);tp3=20.8:0.1:24;x3=-polyval(c3,tp3);y1=0.1*trapz(x1)y2=0.1*trapz(x2)y12=0.1*trapz(x12) y3=0.1*trapz(x3)y=(y1+y2+y12+y3)*237.8*0.01程序shuita6.mt=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.9 20.84 23.88 24.99 25.91;h=968 948 9

31、31 913 898 881 869 852 839 822 1082 1050 1021 994 965 941 918 892 866 843 822 1059 1035 1018;c1=polyfit(t(1:10),h(1:10),3);a1=polyder(c1);tp1=0:0.1:9;x1=-polyval(a1,tp1);c2=polyfit(t(11:21),h(11:21),4);a2=polyder(c2);tp2=10.9:0.1:21;x2=-polyval(a2,tp2);xx1=-polyval(a1,8 9);xx2=-polyval(a2,11 12);xx12=xx1 xx2;c12=polyfit(8 9 11 12,xx12,3);tp12=9:0.1:11;x12=polyval(c12,tp12);dt3=diff(t(22:24);dh3=diff(h(22:24);dht3=-dh3./dt3;t3=20 20.8 t(22) t(23);