水流量的估計12.1實驗目的

12.2實驗內(nèi)容

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

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

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

求一個一次多項式，使得多項式在結點上滿足條件

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

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

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

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

2.5問題求解

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

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

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

我們總共需要對四段數(shù)據(jù)（第1，2未供水時段，第1供水時段，混合時段）進行插值處理，下面以第1未供水時段數(shù)據(jù)為例分別用三種方法算出流量函數(shù)和用水量（用水高度）。下面是用MATLAB實現(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’)運行結果lglrjf=145.6231；fdxxjf=147.1430；sancytjf=1412.687027

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

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

從表12-2中的測量記錄看，一天有兩次供水時段和三次未供水時段，分別對第1，2未供水時段的測量數(shù)據(jù)直接作多項式擬合，可得到水位函數(shù)（注意，根據(jù)多項式擬合的特點，此處擬合多項式的次數(shù)不宜過高，一般以3～6次為宜）。對第3未供水時段來說，數(shù)據(jù)過少不能得到很好的擬合。設t，h分別為已輸入的時刻和水位測量記錄（由表12.2提供，水泵啟動的4個時刻不輸入），這樣第1未供水時段各時刻的水位可由如下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ù)圖形。38變量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）確定流量—時間函數(shù)對于第1，2未供水時段的流量可直接對水位函數(shù)求導，程序如下：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)；*/計算第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])；*/為下面的程序準備數(shù)據(jù)/*x25=-polyval(a2,[19.96,20.84])；*/為下面的程序準備數(shù)據(jù)/*subplot(1,2,1)plot(tp1,x13*100)*/與圖12.10單位保持一致/*subplot(1,2,2)plot(tp2,x23*100)*/與圖12.10單位保持一致/*40程序運行得到第1，2未供水時段的時間—流量圖如圖12.10，可以看到與圖12.8中用插值給出的曲線比較吻合。41如果用5次多項式擬合則得圖12.11的圖形，顯然較三次擬合的效果好。42

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

分別對供