水塔水流量估計(jì)模型與數(shù)據(jù)插值計(jì)算分析

1、Mathematics Laboratory阮小娥博士Experiments in Mathematics趙小艷數(shù)學(xué)實(shí)驗(yàn)辦公地址：理科樓218實(shí)驗(yàn)14 水塔水流量估計(jì)模型與數(shù)據(jù)插值 2、學(xué)會用matlab軟件進(jìn)行數(shù)據(jù)插值計(jì)算1、掌握四種經(jīng)典的插值方法：拉格朗日插值法、牛頓插值法、分段插值法、三次樣條插值法3、學(xué)會用數(shù)據(jù)插值、數(shù)據(jù)擬合方法建立數(shù)學(xué)模型并求解實(shí)驗(yàn)問題 一條100米寬的河道見圖1所示,為了測量其流量需要知道河道的截面積。為此從一端開始每隔5米測量出河床的深度如表1 試根據(jù)以上數(shù)據(jù)，估計(jì)出河道的截面積，進(jìn)而在已知河流速度（設(shè)為1m/s）的情況下計(jì)算出流量。若沿河床鋪設(shè)一條光纜，試估計(jì)光

2、纜的長度。一 數(shù)據(jù)插值 給定n個(gè)數(shù)據(jù)點(diǎn)稱P(x)為插值函數(shù)。也即求解一條嚴(yán)格通過各數(shù)據(jù)點(diǎn)的曲線，用它來進(jìn)行分析研究和預(yù)測，這種方法常稱為數(shù)據(jù)插值法。對于這類問題，選取一條何種類型的曲線作為插值函數(shù)是求解的關(guān)鍵。P(x)不同，就得到不同的插值函數(shù)。1 拉格朗日插值 令 可以證明，對于n+1個(gè)不同結(jié)點(diǎn)，必存在唯一的次數(shù)不超過n的滿足條件的多項(xiàng)式，這個(gè)多項(xiàng)式稱為插值多項(xiàng)式，這種方法稱為n次多項(xiàng)式插值（或代數(shù)插值。為了以后使用方便，先編制一個(gè)Lagrange插值函數(shù)程序：function p=lagrange(x,y)L=length(x);A=ones(L);for j=2:L A(:,j)=A(:

3、,j-1).*x;endX=inv(A)*y;for i=1:L p(i)=X(L-i+1);end例4.1 已知觀測數(shù)據(jù)x 1 2 3 4 5y-1 1.5 2.1 3.6 4.9求其插值多項(xiàng)式曲線。利用Matlab提供的計(jì)算以向量p為系數(shù)的多項(xiàng)式值的命令polyval , 可以求得多項(xiàng)式函數(shù)在任意一點(diǎn)的值。其使用格式為 y0=polyval(p,x0)其中p為多項(xiàng)式系數(shù)（從高次到底次排列）向量，x0為所求值的點(diǎn)或向量，y0為返回x0處的函數(shù)值，從而可以繪制多項(xiàng)式曲線。程序如下：x=1 2 3 4 5;y=-1 1.5 2.1 3.6 4.9;plot(x,y,b.,markersize,3

4、0)axis(1 5 -1 5)gridhold onp=lagrange(x,y);t=1:0.1:5;u=polyval(p,t);plot(t,u,r-,linewidth,3)從結(jié)果可以看到，所插值的4次多項(xiàng)式曲線較好地連接了5個(gè)數(shù)據(jù)點(diǎn)，從而可以用此多項(xiàng)式曲線作為這5個(gè)數(shù)據(jù)的一個(gè)近似變化。2 牛頓插值法 拉格朗日插值法的最大缺點(diǎn)在于當(dāng)增加差值點(diǎn)時(shí)，需要重新計(jì)算多項(xiàng)式的系數(shù)，沒有承接性。下面的牛頓插值法就避免了這個(gè)問題。差商具有如下性質(zhì)：（1）m階差商是零階差商的線性組合；（2）差商與插值結(jié)點(diǎn)的次序無關(guān)；（3）若f(x)是m次多項(xiàng)式，則 fx0,x1,xm=0由差商公式，逐次代入得稱為牛

5、頓插值公式，最后一項(xiàng)稱為牛頓插值余項(xiàng)，記為Rn(x)，余項(xiàng)前的多項(xiàng)式稱為插值多項(xiàng)式，記為Pn(x)。牛頓插值多項(xiàng)式具有以下特點(diǎn)：（1）在插值結(jié)點(diǎn)處與拉格朗日插值一樣，誤差為零；（2）多項(xiàng)式k次項(xiàng)的系數(shù)是f(x)的k階差商；（3）增加插值節(jié)點(diǎn)時(shí)，只增加最后一項(xiàng)，不必像拉格朗日插值公式那樣需要重新計(jì)算系數(shù)。 在做牛頓插值時(shí)，一般先做出差商表，然后套用公式。3 樣條插值法 特別地，如果m=1，則稱為分段線性插值，即在每一個(gè)子區(qū)間上S(x)是線性函數(shù)，在整個(gè)區(qū)間上是分段線性函數(shù)。這種情況一般很難滿足實(shí)際要求，通常使用最多的是3次樣條插值函數(shù)，即S(x)在每一個(gè)子區(qū)間上是三次多項(xiàng)式函數(shù)，而且在每個(gè)結(jié)點(diǎn)處

6、滿足二階導(dǎo)數(shù)連續(xù)和相關(guān)的邊界條件。二 MATLAB軟件實(shí)現(xiàn)數(shù)據(jù)插值1 一維插值命令 yb=interp1(x,y,xb,method)x,y是同維數(shù)據(jù)向量，表示插值結(jié)點(diǎn)的橫坐標(biāo)和縱坐標(biāo)。若x是向量，y是矩陣，則對y的每一列與x配對進(jìn)行插值；xb表示待求函數(shù)值的插值結(jié)點(diǎn)向量，可以缺?。籱ethod是可選項(xiàng)，說明插值使用的方法。缺省時(shí)為線性插值，也可選擇：nearest（最近插值），linear（線性），spline（三次樣條），cubic（三次插值）。命令返回值yb是插值曲線在xb處的縱坐標(biāo)值。2 二維插值命令zb=interp2(x,y,z,xb,yb,method) 根據(jù)同維數(shù)據(jù)向量x，y，

7、z，按照指定的方法做插值運(yùn)算，然后返回插值函數(shù)的豎坐標(biāo)。3 三維插值命令vb=interp3(x,y,z,v,xb,yb,zb,method)4 樣條插值命令yb=spline(x,y,xb)該命令等同于yb=interp1(x,y,xb,cubic)例4.2 已知觀測數(shù)據(jù) 分別用拉格朗日、分段線性、3次樣條進(jìn)行插值，并繪出插值多項(xiàng)式曲線圖。x0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1y-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2(1)拉格朗日插值法x=0:0.1:1；y=-0.447 1.

8、978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2;plot(x,y,b.,markersize,30)axis(0 1 -2 16)gridhold onp=lagrange(x,y);t=0:0.01:1;u=polyval(p,t);plot(t,u,r-,linewidth,3) 結(jié)果顯示，所插值出的10次多項(xiàng)式曲線，在數(shù)據(jù)點(diǎn)之間產(chǎn)生較大的起伏波動，與數(shù)據(jù)點(diǎn)的變化趨勢有明顯的偏離，這時(shí)曲線并不能很好地反映數(shù)據(jù)點(diǎn)的變化規(guī)律。而且進(jìn)一步實(shí)驗(yàn)發(fā)現(xiàn)，隨著分點(diǎn)的增加，Lagrange插值出現(xiàn)大的起伏波動越明顯，這就是插值問題中典型的“龍格現(xiàn)象”。 針對

9、高次多項(xiàng)式插值時(shí)容易發(fā)生“龍格現(xiàn)象”，在實(shí)際插值時(shí)，常常采用分段低次插值方法，即在相鄰兩個(gè)數(shù)據(jù)點(diǎn)構(gòu)成的子區(qū)間上分別進(jìn)行低次插值, 整個(gè)區(qū)間上的插值函數(shù)將是一個(gè)分段的多項(xiàng)式函數(shù)。 (2)分段線性插值x=0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1;y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2;plot(x,y,b.,markersize,30)；axis(0 1 -1 15)；gridhold ont=0:0.01:1;u=interp1(x,y,t); plot(t,u,r-,linewi

10、dth,3) 結(jié)果分析 分段線性插值有效地回避了插值問題中的“龍格現(xiàn)象”，結(jié)果連線也大致描述了已知數(shù)據(jù)點(diǎn)的變化規(guī)律。但很明顯，由分段直線連接的插值曲線在節(jié)點(diǎn)處不光滑，不可導(dǎo)。(3)3次樣條插值x=0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1;y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2;plot(x,y,b.,markersize,30)axis(0 1 -1 16)gridhold onpause(1)t=0:0.01:1;u=spline(x,y,t);plot(t,u,r-,line

11、width,3)結(jié)果分析 從圖可以看出，樣條插值所得的曲線比較好地連接了已知的數(shù)據(jù)點(diǎn)，既有效地回避了插值問題中的“龍格現(xiàn)象”，又是連續(xù)光滑的，用此曲線來近似描述已知數(shù)據(jù)點(diǎn)的變化規(guī)律，能比較好地進(jìn)行數(shù)據(jù)點(diǎn)之間的預(yù)測分析和求值。三 實(shí)驗(yàn)問題求解1 畫出河床觀測點(diǎn)的散點(diǎn)圖x=0:5:100;y=0 2.41 2.96 2.15 2.65 3.12 4.23 5.12 6.21 5.68 4.22 3.91 3.26 2.85 2.35 3.02 3.63 4.12 3.46 2.08 0;axis(0 100 0 10)gridy1=10-y;plot(x,y1,b.,markersize,30)2

12、 利用分段線性插值繪制河床曲線 根據(jù)已給的數(shù)據(jù)可以進(jìn)行分段線性插值, 在此基礎(chǔ)上，利用梯形法求積分命令trapz來計(jì)算河床截面積，同時(shí)計(jì)算每一段連接線長度之和來近似河床曲線長度。x=0:5:100;y=0 2.41 2.96 2.15 2.65 3.12 4.23 5.12 6.21 5.68 4.22 3.91 3.26 2.85 2.35 3.02 3.63 4.12 3.46 2.08 0;y1=10-y;plot(x,y1,b.,markersize,30);axis(0 100 0 10)gridhold ont=0:100;u=interp1(x,y1,t);plot(t,u,r,

13、linewidth, 3);S=100*10-trapz(x,y1);p=sqrt(diff(x).2+diff(y1).2);L=sum(p);fprintf(S=%.2f,L=%.2fn,S,L)從運(yùn)行結(jié)果得到：河床截面面積：S=337.15，河床曲線長度： L=102.09(3) 樣條插值 另一方面, 為了提高河床曲線的模擬精度，可以根據(jù)已給的數(shù)據(jù)進(jìn)行3次樣條插值, 在此基礎(chǔ)上，利用梯形法求積分命令trapz來計(jì)算河床截面積，同時(shí)對樣條曲線加密，計(jì)算每一段連接線長度之和來近似河床曲線長度。x=0:5:100;y=0 2.41 2.96 2.15 2.65 3.12 4.23 5.12 6

14、.21 5.68 4.22. 3.91 3.26 2.85 2.35 3.02 3.63 4.12 3.46 2.08 0;y1=10-y;plot(x,y1,b.,markersize,30)axis(0 100 2 10)gridhold ont=0:0.5:100;u=spline(x,y1,t);plot(t,u,r,linewidth, 3)；S=100*10-trapz(t,u);p=sqrt(diff(t).2+diff(u).2);L=sum(p);fprintf(S=%.2f,L=%.2fn,S,L)河床截面面積：S=339.50，河床曲線長度：L=102.23.數(shù)據(jù)插值模型

15、試驗(yàn)：水塔水流量估計(jì)一實(shí)驗(yàn)問題 某一地區(qū)的用水管理機(jī)構(gòu)要求各社區(qū)提供每小時(shí)以加侖計(jì)的用水量以及每天所用水的總量。由于許多社區(qū)沒有測量流入或流出水塔的水量裝置，所以他們只能通過測量水塔每小時(shí)的水位高度來代替每小時(shí)的用水量(誤差不超過0.05)。當(dāng)水塔中的水位下降到最低水位L (27.00ft)時(shí)水泵就自動啟動向水塔輸水, 當(dāng)水塔水位達(dá)到限定最高水位H (35.00ft)時(shí)，水泵停止供水。水泵供水期間無法測量水泵的供水量。水泵每天輸水一次或兩次,每次約二小時(shí)。 已知某一小鎮(zhèn)的水塔是一高為40ft，直徑為57ft的正圓柱，下表記錄的是某一天水塔水位的真實(shí)數(shù)據(jù)(1ft=0.3024米(m)。試估計(jì)任何

16、時(shí)刻(包括水泵正在輸水時(shí)間)從水塔流出的水流量和一天的用水總量。二、問題分析 流量是單位時(shí)間內(nèi)流出水的體積，由于水塔是正圓柱形，截面面積是常數(shù)，所以，流量很容易根據(jù)水位相對時(shí)間的變化率算出。 水泵不工作時(shí)段的用水量可以由測量記錄直接得到:由表1中記錄的水位下降高度乘以水塔的截面面積就是這一時(shí)段的用水量。 問題的難點(diǎn)在于:如何估計(jì)水泵供水時(shí)段的流量。 水泵供水時(shí)段的流量只能靠供水時(shí)段前后的流量經(jīng)插值或擬合得到。而作為用于插值或擬合的原始數(shù)據(jù)，我們希望水泵不工作時(shí)段的流量數(shù)據(jù)越準(zhǔn)確越好。這個(gè)數(shù)值可以用來檢驗(yàn)數(shù)據(jù)插值或擬合結(jié)果的準(zhǔn)確性。二、問題分析流量大體上可由兩種方法計(jì)算: 直接對表1中的用水量用

17、數(shù)值插值方法算出各時(shí)段的流量，用它們擬合其它時(shí)刻或連續(xù)時(shí)間的流量; 先用表中數(shù)據(jù)擬合水位-時(shí)間函數(shù),求導(dǎo)數(shù)即可得到連續(xù)時(shí)間的流量。 有了任何時(shí)刻的流量，就不難計(jì)算一天的總用水量。三、問題求解 為了表示方便，將所給表1中的數(shù)據(jù)全部化為國際標(biāo)準(zhǔn)單位，時(shí)間用小時(shí)(h),高度用米(m)。3.1 模型假設(shè)(1) 流量只取決于水位差，與水位本身無關(guān)。(2) Torricelli定律：從小孔流出液體的流速正比于液面高度的平方根。 在假設(shè)出口的水位為零的前提下，題目給出水塔的最低和最高水位分別是8.1648m(270.3024) 和10.7352m(35.500.3024)。 因?yàn)椋簊qrt(10.7352/

18、8.1648)=1.14671所以，可忽略水位對流速的影響。(3) 流量可以看作是時(shí)間的連續(xù)函數(shù)。因此，為了計(jì)算簡單，不妨將流量定義成單位時(shí)間流出水的高度，即水位對時(shí)間變化率的絕對值。(4) 水塔截面面積為3.2 流量估計(jì)方法由表2給出的數(shù)據(jù)，用MTALAB軟件做出時(shí)間水位散點(diǎn)圖1：為了計(jì)算水箱水流量與時(shí)間的關(guān)系，最簡單的方法就是： 根據(jù)給出的散點(diǎn)數(shù)據(jù)圖，將數(shù)據(jù)分為三段，然后對每一段數(shù)據(jù)用數(shù)據(jù)插值或者數(shù)據(jù)擬合的方法得到流量與時(shí)間的近似函數(shù)關(guān)系。具體方法如下：中點(diǎn)流速=(左端點(diǎn)的水位-右端點(diǎn)的水位)/時(shí)間區(qū)間長度每段數(shù)據(jù)首尾點(diǎn)的流速用下面的公式計(jì)算：經(jīng)過計(jì)算，得到時(shí)間與流速之間的關(guān)系數(shù)據(jù)表3 數(shù)

19、據(jù)表3對應(yīng)的時(shí)間-流速散點(diǎn)圖如下：下面分別用數(shù)據(jù)插值和數(shù)據(jù)擬合兩種方法來估計(jì)水塔水流量。(1) 數(shù)據(jù)插值法 采用拉格朗日、分段線性和三次樣條三種數(shù)據(jù)插值方法。根據(jù)表3的數(shù)據(jù)，對水泵不工作時(shí)段1,2通過采取數(shù)據(jù)插值方法可以得到任意時(shí)刻的流速，從而知道任意時(shí)刻的流量。對于水泵工作時(shí)段1，應(yīng)用前后時(shí)期的流速進(jìn)行插值。由于最后水泵工作時(shí)段2的數(shù)據(jù)太少，我們將它與水泵不工作時(shí)段3合并，一同進(jìn)行數(shù)據(jù)插值處理（簡稱混合時(shí)段）。以第1未供水時(shí)段數(shù)據(jù)為例分別用三種方法算出流量函數(shù)和用水量（用水高度）。t=0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97

20、;v=29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.27;t0=0:0.1:8.97;lglr=lglrcz(t,v,t0); %拉格朗日插值，需要定義m文件lglrjf=0.1*trapz(lglr) %梯形法算數(shù)值積分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

21、,b)gtext(lglr)gtext(fdxx)gtext(scyt)計(jì)算結(jié)果：lglrjf=145.6231，fdxxjf=147.1430 ，sancytjf=145.6870function y=lglrcz(x0,y0,x)n=length(x0);m=length(x);for i=1:m z=x (i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end由表2知，第1未供水時(shí)段的總用水高度為146(=968-822)。上述三

22、種插值方法計(jì)算的結(jié)果與實(shí)際值146相比都比較接近。計(jì)算結(jié)果：lglrjf=145.6231，fdxxjf=147.1430 ，sancytjf=145.6870(2) 數(shù)據(jù)擬合法1) 擬合水位-時(shí)間函數(shù)2) 確定流量時(shí)間函數(shù)3) 一天總用水量的估計(jì)4) 流量及總用水量的檢驗(yàn)從表2中測量記錄看，一天有兩次供水時(shí)段和三次未供水時(shí)段，分別對第1,2未供水時(shí)段的測量數(shù)據(jù)直接作多項(xiàng)式擬合，可得到水位函數(shù)。第3未供水時(shí)段數(shù)據(jù)過少,擬合效果不是很理想。應(yīng)當(dāng)注意:根據(jù)多項(xiàng)式擬合的特點(diǎn)，此處擬合多項(xiàng)式的次數(shù)不宜過高，一般以3-6次為宜）。1) 擬合水位-時(shí)間函數(shù)t=0,0.92,1.84,2.95,3.87,4

23、.98,5.90,7.00,7.93,8.97,10.9512.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.9620.84,23.88,24.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.8210.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.5910.35,10.18;figure(1)c1=polyfit(t(1:10),h(1:10),3);tp1=0:0.1:8.9;x1=polyval(c1,t

24、p1);%變量x1中存放了以0.1為步長算出的各個(gè)時(shí)刻的水位高度。plot(tp1,x1)figure(2)c2=polyfit(t(11:21),h(11:21),3);tp2=10.9:0.1:20.9;x2=-polyval(c2,tp2); plot(tp2,x2)第1、2未供水時(shí)段時(shí)間與水位圖程序：2) 確定流量時(shí)間函數(shù)第1,2未供水時(shí)段第1供水時(shí)段第2供水時(shí)段第3未供水時(shí)段混合時(shí)段第1,2未供水時(shí)段的流量可直接對水位函數(shù)求導(dǎo)得到，用三次多項(xiàng)式擬合曲線，程序如下：c1=polyfit(t(1:10),h(1:10),3); c2=polyfit(t(11:21),h(11:21),3

25、);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);x14=-polyval(a1,7.93,8.97);% (為后面的計(jì)算準(zhǔn)備數(shù)據(jù))x23=-polyval(a2,tp2); x114=-polyval(a2,10.95:0.01:20.84);wgsysl2=100*trapz(tp2,x114);x24=-polyval(a2,10.95,1

26、2.03); %(為后面的計(jì)算準(zhǔn)備數(shù)據(jù))x25=-polyval(a2,19.96,20.84); %(為后面的計(jì)算準(zhǔn)備數(shù)據(jù))subplot(1,2,1)plot(tp1,x13*100)subplot(1,2,2)plot(tp2,x23*100)上下曲線比較發(fā)現(xiàn)，用三次多項(xiàng)式擬合效果不是很好，故改用五次多項(xiàng)式重新擬合。 對于第1供水時(shí)段的流量,則用前后時(shí)期的流量進(jìn)行擬合得到。為使流量函數(shù)在t=8.97和t=10.93連續(xù)，只取4個(gè)點(diǎn)，用3次多項(xiàng)式擬合得到第1供水時(shí)段的時(shí)間-流量圖。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);plot(nhsj, 100*nhlsjg) 在第2供水時(shí)段之前取t=19.96