水流量的估計(jì)課件

水流量的估計(jì)12.1實(shí)驗(yàn)?zāi)康?/p>

12.2實(shí)驗(yàn)內(nèi)容

美國某州的用水管理機(jī)構(gòu)要求各社區(qū)提供以每小時(shí)多少加侖計(jì)的用水量以及每天所用的總水量。許多社區(qū)沒有測(cè)量流入或流出水塔水量的裝置，只能代之以每小時(shí)測(cè)量水塔中的水位，其誤差不超過5%。需要注意的是，當(dāng)水塔中的水位下降到最低水位L時(shí)，水泵就自動(dòng)向水塔輸水直到最高水位H，此期間不能測(cè)量水泵的供水量，因此，當(dāng)水泵正在輸水時(shí)不容易建立水塔中水位和用水量之間的關(guān)系。水泵每天輸水一次或兩次，每次約2小時(shí)。

本實(shí)驗(yàn)的主要目的是使學(xué)生學(xué)會(huì)用MATLAB軟件進(jìn)行插值和曲線擬合計(jì)算并解決一些具體的實(shí)際問題。通過實(shí)際問題的解決，使學(xué)生了解如何利用曲線插值及曲線擬合解決實(shí)際問題的全過程。2.1實(shí)驗(yàn)問題1

試估計(jì)任何時(shí)刻（包括水泵正在輸水時(shí)間）從水塔流出的水流量f(t)，并估計(jì)一天的總用水量。已知該水塔是一個(gè)高為40ft（英尺），直徑為57ft（英尺）的正圓柱，表5-1給出了某個(gè)小鎮(zhèn)一天水塔水位的真實(shí)數(shù)據(jù)，水位降至約27.00ft水泵開始工作，水位升到35.50ft停止工作。（注：1ft（英尺）=0.3048m（米））2表12.1某小鎮(zhèn)某天水塔水位時(shí)間/s水位/0.01ft時(shí)間/s水位/0.01ft03175466363350331631104995332606635305453936316710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254水泵開動(dòng)35932水泵開動(dòng)82649水泵開動(dòng)39332水泵開動(dòng)8596834753943535508995333974331834459327033403有了任何時(shí)刻的流量，就不難計(jì)算一天的總用水量。其實(shí)，水泵不工作時(shí)段的用水量可以由測(cè)量記錄直接得到，由表12-1中下降水位乘以水搭的截面積就是這一時(shí)段的用水量。這個(gè)數(shù)值可以用來檢驗(yàn)數(shù)據(jù)插值或擬合的結(jié)果。在具體給出本問題的解答之前，先介紹一個(gè)簡(jiǎn)單的數(shù)據(jù)插值方法。52.3拉格朗日插值1、線性插值假設(shè)已知在區(qū)間上的兩個(gè)結(jié)點(diǎn)和它們的函數(shù)值

求一個(gè)一次多項(xiàng)式，使得多項(xiàng)式在結(jié)點(diǎn)上滿足條件

這種插值方法稱為線性插值方法（也稱兩點(diǎn)插值）?？梢郧蟪觯?2、拋物插值已知在區(qū)間上的三個(gè)結(jié)點(diǎn)和它們的函數(shù)值

求一個(gè)次數(shù)不超過2的多項(xiàng)式，使得它在結(jié)點(diǎn)上滿足條件這種插值方法稱為拋物線插值法，可求出：7上述多項(xiàng)式稱為n次拉格朗日（Lagrange）插值多項(xiàng)式，函數(shù)稱為拉格朗日插值基函數(shù)。當(dāng)n=1,2時(shí)，n次拉格朗日（Lagrange）插值多項(xiàng)式即為線性插值多項(xiàng)式和拋物插值多項(xiàng)式。9例12.1已知函數(shù)發(fā)f(x)的函數(shù)表如下：

求其拉格朗日插值多項(xiàng)式，并求的近似值。解由于給出了4個(gè)插值結(jié)點(diǎn)，所以可做出次數(shù)不超過3的拉格朗日插值多項(xiàng)式。10將上列4式代入n=3的拉格朗日插值公式，可得所要求的插值多項(xiàng)式為將x=2.5代入可得f(2.5)的近似值為1.8496。拉格朗日插值法適合節(jié)點(diǎn)較少的情況，當(dāng)節(jié)點(diǎn)較多的大范圍高次插值的逼近效果往往并不理想且當(dāng)插值結(jié)點(diǎn)增加時(shí)，計(jì)算越來越繁。為了提高精度和減少計(jì)算，還有牛頓插值法下、三次樣條插值等，具體可參閱有關(guān)書籍。112．二維插值命令interp2的具體使用格式zz=interp2(x,y,z,xx,yy,’method’)該指令的意思是根據(jù)數(shù)據(jù)向量x，y，z按method指定的方法來做插值，然后將xx，yy處插值函數(shù)的插值結(jié)點(diǎn)向量，如果xx，yy在插值范圍之內(nèi)，則返回值在zz中，否則返回值為空——NaN?！痬ethod’是插值方法可選項(xiàng)，具體要求同一維插值的情況。該命令還有以下幾種省略格式：zz=interp2(z,xx,yy)zz=interp2(x,y,z,xx,yy)zz=interp2(z,ntimes)133．三維插值命令interp3的具體使用格式vi=interp3(x,y,z,v,xi,yi,zi,’method’)它的具體含義跟前面的一、二維插值是相似的，在此不作解釋，讀者可在MATLAB工作空間中用helpinterp3命令獲得。4．樣條插值命令spline的具體使用格式y(tǒng)y=spline(x,y,xx)它的意思等同于命令yy=interp1(x,y,xx,’cubic’)14例12.2在用外接電源給電容器充電時(shí)，電容器兩端的電壓V將會(huì)隨著充電時(shí)間t發(fā)生變化，已知在某一次實(shí)驗(yàn)時(shí)，通過測(cè)量得到下列觀測(cè)值，分別用拉格朗日插值法、分段線性插值法、三次樣條插值法畫出V隨著時(shí)間t變化的曲線圖，分別計(jì)算當(dāng)時(shí)間t=7s時(shí)，三種插值法各自算得電容器兩端電壓的近似值。

解由于MATLAB沒有提供現(xiàn)成的拉格朗日插值命令，我們可以編寫一個(gè)函數(shù)lglrcz.m來完成，其他兩種插值法可用現(xiàn)成的命令。用MATLAB軟件進(jìn)行三種插值計(jì)算的程序?yàn)閟zczqx.m。15程序?yàn)閍zczf.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=interp1(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’)17執(zhí)行結(jié)果是laglr=9.52988980716254fendxx=9.700sancyt=9.671圖形如圖12.1所示。圖中曲線lglr表示拉格朗日插值曲線，scyt表示三次樣條插值曲線，fdxx表示分段線性插值曲線。18

2.5問題求解

為了表示方便，我們將2.1節(jié)問題中所給表12-1中的數(shù)據(jù)全部化為國際標(biāo)準(zhǔn)單位（表12-2），時(shí)間用小時(shí)（h），高度用米（m）：表12-2一天內(nèi)水塔水位記錄191．模型假設(shè)（1）流量只取決于水位差，與水位本身無關(guān)，故由物理學(xué)中Torriceli定律：小孔流出的液體的流速正比于水面高度的平方根。題目給出水塔的最低和最高水位分別是8.1648m和10.7352m（設(shè)出口的水位為零），因?yàn)閟qrt，約為1，所以可忽略水位對(duì)流速的影響。（2）將流量看作是時(shí)間的連續(xù)函數(shù)，為計(jì)算簡(jiǎn)單，不妨將流量定義成單位時(shí)間流出水的高度，即水位對(duì)時(shí)間變化率的絕對(duì)值（水位是下降的），水塔截面積為（m2），得到結(jié)果后乘以s即可。212．流量估計(jì)方法首先依照表12-2所給數(shù)據(jù)，用MATLAB作出時(shí)間—水位散點(diǎn)圖（圖12.2）。22下面來計(jì)算水箱流量與時(shí)間的關(guān)系。

根據(jù)圖12..，一種簡(jiǎn)單的處理方法為，將表12-2中的數(shù)據(jù)分為三段，然后對(duì)每一段的數(shù)據(jù)做如下處理：設(shè)某段數(shù)據(jù)，相鄰數(shù)據(jù)中點(diǎn)的平均流速用下面的公式（流速=（右端點(diǎn)的水位－右端點(diǎn)的水位）/區(qū)間長(zhǎng)度）：每段數(shù)據(jù)首尾點(diǎn)的流速用下面的公式計(jì)算：用以上公式求得時(shí)間與流速之間的數(shù)據(jù)如表12-3。23由表12-3作出時(shí)間—流速散點(diǎn)圖如圖12.3。1）插值法

由表12-3，對(duì)水泵不工作時(shí)段1，2采取插值方法，可以得到任意時(shí)刻的流速，從而可以知道任意時(shí)刻的流量。25我們分別采取拉格朗日插值法，分段線性插值法及三次樣條插值法；對(duì)于水泵工作時(shí)段1應(yīng)用前后時(shí)期的流速進(jìn)行插值，由于最后一段水泵不工作時(shí)段數(shù)據(jù)太少，我們將它與水泵工作時(shí)段2合并一同進(jìn)行插值處理（該段簡(jiǎn)稱混合時(shí)段）。

我們總共需要對(duì)四段數(shù)據(jù)（第1，2未供水時(shí)段，第1供水時(shí)段，混合時(shí)段）進(jìn)行插值處理，下面以第1未供水時(shí)段數(shù)據(jù)為例分別用三種方法算出流量函數(shù)和用水量（用水高度）。下面是用MATLAB實(shí)現(xiàn)該過程的程序。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；26lglr=lglrcz(t,v,t0)；/*注：lglrcz為一函數(shù)，程序同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’)運(yùn)行結(jié)果lglrjf=145.6231；fdxxjf=147.1430；sancytjf=1412.687027

由表12-2知，第1未供水時(shí)段的總用水高度為146(=968－822)，可見上述三種插值方法計(jì)算的結(jié)果與實(shí)際值（146）相比都比較接近。考慮到三次樣條插值方法具有更加良好的性質(zhì)，建議采取該方法。

其他三段的處理方法與第1未供水時(shí)段的處理方法類似，這里不再詳細(xì)敘述，只給出數(shù)值結(jié)果和函數(shù)圖像（圖12.5～圖12.7），圖中曲線標(biāo)記同圖12.4。29圖12.5第一供水段時(shí)間—流速示意圖30圖12.6第2未供水段時(shí)間—流速示意圖31圖12.7混合時(shí)段時(shí)間—流速示意圖32圖12.8是用分段線性及三次樣條插值方法得到的整個(gè)過程的時(shí)間—流速函數(shù)示意圖。33表12-4各時(shí)段及一天的總用水量（用水高度）第1未供水段第2未供水段第3供水段混合時(shí)段全天拉格朗日插值法145.6231258.866454.068992.1337550.6921分段線性插值法147.1430258.969749.605176.4688532.1866三次樣條插值法145.6870258.654753.333481.7699539.445034表12-5是對(duì)一天中任取的4個(gè)時(shí)刻分別用3種方法得到的水塔水流量近似值。時(shí)間6.8810.8815.8822.88①15.9826671234851433.7426009085346325.5662241818047734.794②14.8272413793103432.9976262626262725.4465591397849525.47③15.82033.75925.54929.41733175863551注：①拉格朗日插值法②分段線性插值法③三次樣條插值法352）擬合法（1）擬合水位—時(shí)間函數(shù)

從表12-2中的測(cè)量記錄看，一天有兩次供水時(shí)段和三次未供水時(shí)段，分別對(duì)第1，2未供水時(shí)段的測(cè)量數(shù)據(jù)直接作多項(xiàng)式擬合，可得到水位函數(shù)（注意，根據(jù)多項(xiàng)式擬合的特點(diǎn)，此處擬合多項(xiàng)式的次數(shù)不宜過高，一般以3～6次為宜）。對(duì)第3未供水時(shí)段來說，數(shù)據(jù)過少不能得到很好的擬合。設(shè)t，h分別為已輸入的時(shí)刻和水位測(cè)量記錄（由表12.2提供，水泵啟動(dòng)的4個(gè)時(shí)刻不輸入），這樣第1未供水時(shí)段各時(shí)刻的水位可由如下MATLAB程序完成：36t=[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)37圖12.9給出的是第1未供水時(shí)段的時(shí)間—水位擬合函數(shù)圖形。38變量x1中存放了以0.1為步長(zhǎng)算出的各個(gè)時(shí)刻的水位高度。同樣地，第2未供水時(shí)段時(shí)間—水位圖可由如下MATLAB程序完成，讀者可自己上機(jī)運(yùn)行查看。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í)間函數(shù)對(duì)于第1，2未供水時(shí)段的流量可直接對(duì)水位函數(shù)求導(dǎo)，程序如下：39c1=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)；*/計(jì)算第1未供水時(shí)段的總用水量/*x14=-polyval(a1,[7.93,8.97])；*/為下面的程序準(zhǔn)備數(shù)據(jù)/*x23=-polyval(a2,tp2)；x114=-polyval(a2,[10.95:0.01:20.84])wgsysl2=100*trapz(tp2,x114)；*/計(jì)算第2未供水時(shí)段的總用水量/*x24=-polyval(a2,[10.95,12.03])；*/為下面的程序準(zhǔn)備數(shù)據(jù)/*x25=-polyval(a2,[19.96,20.84])；*/為下面的程序準(zhǔn)備數(shù)據(jù)/*subplot(1,2,1)plot(tp1,x13*100)*/與圖12.10單位保持一致/*subplot(1,2,2)plot(tp2,x23*100)*/與圖12.10單位保持一致/*40程序運(yùn)行得到第1，2未供水時(shí)段的時(shí)間—流量圖如圖12.10，可以看到與圖12.8中用插值給出的曲線比較吻合。41如果用5次多項(xiàng)式擬合則得圖12.11的圖形，顯然較三次擬合的效果好。42

而第1供水時(shí)段的流量則用前后時(shí)期的流量進(jìn)行擬合得到。為使流量函數(shù)在t=9和t=11連續(xù)，我們只取4個(gè)點(diǎn)，用三次多項(xiàng)式擬合得到第1供水時(shí)段的時(shí)間—流量圖形如圖12.12，可以看到與圖12.8中的相應(yīng)部分比較吻合。圖12.12圖12.843dygsdsy=[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)；*/該語句計(jì)算第1供水時(shí)段的總用水量/*plot(nhsj,100*nhlsjg)程序如下：44在第2供水時(shí)段之前取t=19.96,20.84兩點(diǎn)的流量，用第3未供水時(shí)段的3個(gè)記錄做差分得到兩個(gè)流量數(shù)據(jù)21.62，18.48，然后用這4個(gè)數(shù)據(jù)做三次多項(xiàng)式擬合得到第2供水時(shí)段與第3未供水時(shí)段的時(shí)間—流量圖如圖12.13，可以看到與圖12.8中的相應(yīng)部分也比較吻合。圖12.13，圖12.845程序如下：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=polyval(nhhddxsxs,[20.84:0.01:24])；gsysl2=trapz(gssj2,gs2)；*/該語句計(jì)算第2供水時(shí)段的總用水量/*plot(tp3,xx3)46（3）一天總用水量的估計(jì)

分別對(duì)供