1、實驗實驗14 14 水塔水流量估計模型與數據插值水塔水流量估計模型與數據插值 2 2、學會用、學會用matlabmatlab軟件進行數據插值計算軟件進行數據插值計算1 1、掌握四種經典的插值方法：拉格朗日插值法、掌握四種經典的插值方法：拉格朗日插值法、牛頓插值法、分段插值法、三次樣條插值法牛頓插值法、分段插值法、三次樣條插值法3、學會用數據插值、數據擬合方法建立數學、學會用數據插值、數據擬合方法建立數學模型并求解模型并求解實驗問題實驗問題 一條一條100100米寬的河道米寬的河道見圖見圖1 1所示所示, ,為了測量為了測量其流量需要知道河道其流量需要知道河道的截面積。為此從一的截面積。為此從一

2、端開始每隔端開始每隔5 5米測量出米測量出河床的深度如表河床的深度如表1 1 010203040506070809010001234567 試根據以上數據，估計出河道的截面積，進而在試根據以上數據，估計出河道的截面積，進而在已知河流速度（設為已知河流速度（設為1m/s1m/s）的情況下計算出流量。）的情況下計算出流量。若沿河床鋪設一條光纜，試估計光纜的長度。若沿河床鋪設一條光纜，試估計光纜的長度。),( ,),(),(2211nnyxyxyx一一 數據插值數據插值 給定給定n n個數據點個數據點稱稱P(x)為插值函數，也即求解一條嚴格通過各數為插值函數，也即求解一條嚴格通過各數據點的曲線，用它

3、來進行分析研究和預測，這種據點的曲線，用它來進行分析研究和預測，這種方法常稱為方法常稱為數據插值法數據插值法。對于這類問題，對于這類問題，選取一條何種類型的曲線作為插選取一條何種類型的曲線作為插值函數是求解的關鍵值函數是求解的關鍵。P(x)不同，就得到不同的不同，就得到不同的插值函數。插值函數。nixPyxPxPii, 2 , 1),()()( 滿足滿足，使得，使得數數試求一個較為簡單的函試求一個較為簡單的函1 拉格朗日插值所以滿足所以滿足個數據點都在曲線上，個數據點都在曲線上，因為這因為這次多項式曲線次多項式曲線總可以唯一確定一條總可以唯一確定一條個數據點個數據點對于已給的對于已給的naxa

4、xaxayniyxnnnnii0111,n, 1, 2 , 1),(1n 111102221011110nnnnnnnnnyxaxaayxaxaayxaxaa， nnnnnxxxxxxA112211111， naaaX10， 121nyyyY 令令.1YAXYAX 解為解為的解。的解。方程組方程組則所求的多項式系數為則所求的多項式系數為 可以證明，對于可以證明，對于n+1n+1個不同結點，必存在唯一的次個不同結點，必存在唯一的次數不超過數不超過n n的滿足條件的多項式，的滿足條件的多項式，這個多項式稱為這個多項式稱為插值多項式，插值多項式，這種方法稱為這種方法稱為n n次多項式插值（或代次多項

5、式插值（或代數插值）。數插值）。為了以后使用方便，先編制一個為了以后使用方便，先編制一個Lagrange插值函數程序：插值函數程序：function p=lagrange(x,y)L=length(x);A=ones(L); %產生元素全是產生元素全是1的的L階方陣階方陣for j=2:L A(:,j)=A(:,j-1).*x;endX=inv(A)*y;for i=1:L p(i)=X(L-i+1);end例例4.1 已知觀測數據已知觀測數據x 1 2 3 4 5y-1 1.5 2.1 3.6 4.9求其插值多項式曲線。求其插值多項式曲線。利用利用Matlab提供的計算以向量提供的計算以向量

6、p p為系數的多項式值為系數的多項式值的命令的命令polyval , , 可以求得多項式函數在任意一點的值。可以求得多項式函數在任意一點的值。其使用格式為其使用格式為 y0=polyval(p,x0)其中其中p p為多項式系數（從高次到底次排列）向量，為多項式系數（從高次到底次排列）向量，x0 x0為所為所求值的點或向量，求值的點或向量，y0y0為返回為返回x0 x0處的函數值，從而可以繪制處的函數值，從而可以繪制多項式曲線。多項式曲線。程序如下：程序如下：x=1 2 3 4 5;y=-1 1.5 2.1 3.6 4.9;plot(x,y,b.,markersize,30)axis(1 5 -

7、1 5)gridhold onp=lagrange(x,y);t=1:0.1:5;u=polyval(p,t);plot(t,u,r-,linewidth,3)從結果可以看到，所插值的從結果可以看到，所插值的4 4次多項式曲線較好地連接了次多項式曲線較好地連接了5 5個個數據點，從而可以用此多項式曲線作為這數據點，從而可以用此多項式曲線作為這5 5個數據的一個近個數據的一個近似變化。似變化。11.522.533.544.55-10123452 牛頓插值法 拉格朗日插值法的最大缺點在于當增加差值點時，需拉格朗日插值法的最大缺點在于當增加差值點時，需要重新計算多項式的系數，沒有承接性。下面的牛頓插

8、要重新計算多項式的系數，沒有承接性。下面的牛頓插值法就避免了這個問題。值法就避免了這個問題。為零階差商。為零階差商。數值數值特別地，稱結點處的函特別地，稱結點處的函的一階差商，記為的一階差商，記為關于關于為為上的平均變化率上的平均變化率在區間在區間稱稱)(.,)()()(,)(11111iiiiiiiiiiixfxxfxxxfxxxfxfxxxf 0101111011021212121,1,)(.,)(,xxxxfxxfxxxfmxxxxfxxxfxxxxfxxxxfxxfmmmmmiiiiiiiiiiii 階差商階差商處的處的在點在點類似定義類似定義的二階差商，記為的二階差商，記為關于點關于

9、點為為率率稱一階差商的平均變化稱一階差商的平均變化差商具有如下性質：差商具有如下性質：（1 1）m m階差商是零階差商的線性組合；階差商是零階差商的線性組合；（2 2）差商與插值結點的次序無關；）差商與插值結點的次序無關；（3 3）若）若f(x)f(x)是是m m次多項式，則次多項式，則 fx0,x1,fx0,x1,xm=0,xm=0由差商公式，由差商公式，),(,),(,)()(1210100000 xxxxxfxxfxxfxxxxfxfxf )()(,)()(,)()(,)(,)(,)()(1010110102103210102100100nnnnxxxxxxxxxxfxxxxxxxxxf

10、xxxxxxxxxxfxxxxxxxfxxxxfxfxf 逐次代入得逐次代入得稱為稱為牛頓插值公式牛頓插值公式，最后一項稱為，最后一項稱為牛頓插值余項牛頓插值余項，記為記為Rn(x)，余項前的多項式稱為，余項前的多項式稱為插值多項式插值多項式，記，記為為Pn(x)。牛頓插值多項式具有以下特點：牛頓插值多項式具有以下特點：（1 1）在插值結點處與拉格朗日插值一樣，誤差）在插值結點處與拉格朗日插值一樣，誤差為零；為零；（2 2）多項式）多項式k k次項的系數是次項的系數是f(x)f(x)的的k k階差商；階差商；（3 3）增加插值節點時，只增加最后一項，不必）增加插值節點時，只增加最后一項，不必像

11、拉格朗日插值公式那樣需要重新計算系數。像拉格朗日插值公式那樣需要重新計算系數。結點結點零階差商零階差商一階差商一階差商二階差商二階差商三階差商三階差商0 x1x3x2x0()f x1( )f x2()f x3( )f x01 ,f x x12 ,f x x23,f x x012 , ,f x x x123 ,f x x x0123 ,f x x x x 在做牛頓插值時，一般先做出差商表，然后在做牛頓插值時，一般先做出差商表，然后套用公式。套用公式。3 3 樣條插值法樣條插值法插值函數。插值函數。次樣條次樣條上關于已知劃分的上關于已知劃分的在區間在區間是是則稱則稱，有，有在各個結點上的函數值在各

12、個結點上的函數值函數函數次樣條插值；若同時對次樣條插值；若同時對是關于已知劃分的是關于已知劃分的則稱則稱）（次多項式函數；次多項式函數；上是上是在每一個子區間在每一個子區間）（上滿足條件上滿足條件在在函數函數的一個劃分，如果的一個劃分，如果是區間是區間設設mbaxfxSxfxSxfxfmxSbaCxSmxxxSbaxSbabxxxaiiimiin,)()()()()()()(;,)(2,)(1,)(,)1(110 特別地，如果特別地，如果m=1m=1，則稱為，則稱為分段線性插值，分段線性插值，即在每一個子區間上即在每一個子區間上S(x)S(x)是線性函數，在整個是線性函數，在整個區間上是分段線

13、性函數。這種情況一般很難滿區間上是分段線性函數。這種情況一般很難滿足實際要求，通常使用最多的是足實際要求，通常使用最多的是3 3次樣條插值次樣條插值函數函數，即，即S(x)S(x)在每一個子區間上是三次多項式在每一個子區間上是三次多項式函數，而且在每個結點處滿足二階導數連續和函數，而且在每個結點處滿足二階導數連續和相關的邊界條件。相關的邊界條件。二二 MATLABMATLAB軟件實現數據插值軟件實現數據插值1 1 一維插值命令一維插值命令 yb=interp1(x,y,xb,method)x,yx,y是同維數據向量，表示插值結點的橫坐標和縱是同維數據向量，表示插值結點的橫坐標和縱坐標。若坐標。

14、若x x是向量，是向量，y y是矩陣，則對是矩陣，則對y y的每一列與的每一列與x x配對進行插值；配對進行插值；xbxb表示待求函數值的插值結點向量，可以缺省；表示待求函數值的插值結點向量，可以缺省；methodmethod是可選項，說明插值使用的方法。是可選項，說明插值使用的方法。缺省缺省時為線性插值，時為線性插值，也可選擇：也可選擇：nearestnearest（最近插（最近插值），值），linearlinear（線性），（線性），splinespline（三次樣條），（三次樣條），cubiccubic（三次插值）。（三次插值）。命令返回值命令返回值ybyb是插值曲線在是插值曲線在xbx

15、b處的縱坐標值。處的縱坐標值。2 2 二維插值命令二維插值命令zb=interp2(x,y,z,xb,yb,method) 根據同維數據向量根據同維數據向量x，y，z，按照指定的方，按照指定的方法做插值運算，然后返回插值函數的豎坐標。法做插值運算，然后返回插值函數的豎坐標。3 3 三維插值命令三維插值命令vb=interp3(x,y,z,v,xb,yb,zb,method)4 4 樣條插值命令樣條插值命令yb=spline(x,y,xb)該命令等同于該命令等同于yb=interp1(x,y,xb,cubic)zb = griddata(x,y,z,xb,yb,method) x,y,zx,y,

16、z是同維數據向量，表示插值結點的橫坐標、縱是同維數據向量，表示插值結點的橫坐標、縱坐標和豎坐標。坐標和豎坐標。xbxb，ybyb表示待求函數值的插值結表示待求函數值的插值結點向量，通常由點向量，通常由meshgridmeshgrid生成。生成。methodmethod是可選項，說明插值使用的方法。是可選項，說明插值使用的方法。缺省時缺省時為線性插值，為線性插值，也可選擇：也可選擇：nearestnearest（最近插值），（最近插值），linearlinear（線性），（線性），cubiccubic（三次插值）（三次插值），v4v4。 v4v4是一種插值算法，沒有具體的名字，原文稱為是一種插值

17、算法，沒有具體的名字，原文稱為“MATLAB 4 griddata methodMATLAB 4 griddata method”，是一種很圓滑的，是一種很圓滑的差值算法，效果不錯。差值算法，效果不錯。 命令返回值命令返回值zbzb是插值曲面在是插值曲面在(xb,yb)(xb,yb)處的豎坐標值。處的豎坐標值。5 用用griddata命令插值生成三維曲面命令插值生成三維曲面 interp2和和griddata都是二維插值，兩都是二維插值，兩者的區別是：者的區別是： interp2的插值數據必須是矩形域，的插值數據必須是矩形域， 即即已知數據點已知數據點(x,y)組成規則的矩陣，或稱之為組成規則

18、的矩陣，或稱之為柵格，可使用柵格，可使用meshgid生成。生成。 而而griddata函數的已知數據點（函數的已知數據點（X，Y）不要求規則排列，特別是對試驗中隨機沒有不要求規則排列，特別是對試驗中隨機沒有規律采取的數據進行插值具有很好的效果。規律采取的數據進行插值具有很好的效果。例例4.2 4.2 已知觀測數據已知觀測數據 分別用拉格朗日、分段線性、分別用拉格朗日、分段線性、3 3次樣條進行次樣條進行插值，并繪出插值多項式曲線圖。插值，并繪出插值多項式曲線圖。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.

19、08 7.34 7.66 9.56 9.48 9.3 11.2(1)拉格朗日插值法拉格朗日插值法x=0:0.1: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 -2 16)gridhold onp=lagrange(x,y);t=0:0.01:1;u=polyval(p,t);plot(t,u,r-,linewidth,3) 結果顯示，所插值出的結果顯示，所插值出的1010次多項式曲線，在次多項式曲線，在數據點之間產生較大的起伏波動，與數據點的變數據點

20、之間產生較大的起伏波動，與數據點的變化趨勢有明顯的偏離，這時曲線并不能很好地反化趨勢有明顯的偏離，這時曲線并不能很好地反映數據點的變化規律。而且進一步實驗發現，隨映數據點的變化規律。而且進一步實驗發現，隨著分點的增加，著分點的增加，LagrangeLagrange插值出現大的起伏波動插值出現大的起伏波動越明顯，這就是插值問題中典型的越明顯，這就是插值問題中典型的“龍格現象龍格現象”。 針對高次多項式插值時容易發生針對高次多項式插值時容易發生“龍格現龍格現象象”，在實際插值時，常常采用分段低次插值方，在實際插值時，常常采用分段低次插值方法，即在相鄰兩個數據點構成的子區間上分別進法，即在相鄰兩個數

21、據點構成的子區間上分別進行低次插值行低次插值, , 整個區間上的插值函數將是一個分整個區間上的插值函數將是一個分段的多項式函數。段的多項式函數。 (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-,linewidth,3) 結果分

22、析結果分析 分段線性插值有效地回避了插值問題中的分段線性插值有效地回避了插值問題中的“龍格現象龍格現象”，結果連線也大致描述了已知數，結果連線也大致描述了已知數據點的變化規律。但很明顯，由分段直線連接據點的變化規律。但很明顯，由分段直線連接的插值曲線在節點處不光滑，不可導。的插值曲線在節點處不光滑，不可導。(3)3(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)a

23、xis(0 1 -1 16)gridhold ont=0:0.01:1;u=spline(x,y,t);plot(t,u,r-,linewidth,3)結果分析結果分析 從圖可以看出，樣條插值所得的曲線比較從圖可以看出，樣條插值所得的曲線比較好地連接了已知的數據點，既有效地回避了插好地連接了已知的數據點，既有效地回避了插值問題中的值問題中的“龍格現象龍格現象”，又是連續光滑的，又是連續光滑的，用此曲線來近似描述已知數據點的變化規律，用此曲線來近似描述已知數據點的變化規律，能比較好地進行數據點之間的預測分析和求值。能比較好地進行數據點之間的預測分析和求值。例例4.3 某海域上頻繁地有各種噸位的船

24、只經過，為某海域上頻繁地有各種噸位的船只經過，為保證船只的航行安全，有關機構在低潮時對水深保證船只的航行安全，有關機構在低潮時對水深進行了測量，得到了如下的測量數據：進行了測量，得到了如下的測量數據： x=129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5;y=7.5 141.5 23.0 147.0 22.5 137.5 85.5 60.5 121.0 3.0 56.5 116.5 84.0 43.5;Z=4 8 6 8 6 8 8 9 9 8 8 9 4 9;x,y為平面坐標，為平

25、面坐標，z為高度。請根據測量的數據描為高度。請根據測量的數據描述該海域的地貌。述該海域的地貌。 程序如下：程序如下：clc;x=129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5;y=7.5 141.5 23.0 147.0 22.5 137.5 85.5 60.5 121.0 3.0 56.5 116.5 84.0 43.5;z=4 8 6 8 6 8 8 9 9 8 8 9 4 9;h=-z;xi=linspace(min(x),max(x),100);yi=linspace(mi

26、n(y),max(y),100);X,Y=meshgrid(xi,yi);H=griddata(x,y,h,X,Y,v4);mesh(X,Y,H);三三 實驗問題求解實驗問題求解1 1 畫出河床觀測點的散點圖畫出河床觀測點的散點圖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)0102030405060708090

27、1000123456789102 2 利用分段線性插值繪制河床曲線利用分段線性插值繪制河床曲線 根據已給的數據可以進行分段線性插值根據已給的數據可以進行分段線性插值, , 在此基礎上，利用梯形法求積分命令在此基礎上，利用梯形法求積分命令trapztrapz來計算河床截面積，同時計算每來計算河床截面積，同時計算每一段連接線長度之和來近似河床曲線長一段連接線長度之和來近似河床曲線長度。度。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.

28、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,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)從運行結果得到：從運行結果得到：河床截面面積：河床截面面積：S=337.15，河床曲線長度：河床曲線長度： L=102.09010203040506070809010001234567891

29、03 3 利用樣條插值計算截面積與曲線長度利用樣條插值計算截面積與曲線長度 另一方面另一方面, , 為了提高河床曲線的模擬為了提高河床曲線的模擬精度，可以根據已給的數據進行精度，可以根據已給的數據進行3 3次樣條次樣條插值插值, , 在此基礎上，利用梯形法求積分在此基礎上，利用梯形法求積分命令命令trapztrapz來計算河床截面積，同時對樣來計算河床截面積，同時對樣條曲線加密，計算每一段連接線長度之條曲線加密，計算每一段連接線長度之和來近似河床曲線長度。和來近似河床曲線長度。x=0:5:100;y=0 2.41 2.96 2.15 2.65 3.12 4.23 5.12 6.21 5.68

30、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

31、.01020304050607080901002345678910數據插值模型試驗：水塔水流量估計數據插值模型試驗：水塔水流量估計一實驗問題一實驗問題 某一地區的用水管理機構要求各社區提供每小時以加某一地區的用水管理機構要求各社區提供每小時以加侖計的用水量以及每天所用水的總量。由于許多社區沒侖計的用水量以及每天所用水的總量。由于許多社區沒有測量流入或流出水塔的水量裝置，所以他們只能通過有測量流入或流出水塔的水量裝置，所以他們只能通過測量水塔每小時的水位高度來代替每小時的用水量測量水塔每小時的水位高度來代替每小時的用水量( (誤誤差不超過差不超過0.05)0.05)。當水塔中的水位下降到最低水位

32、。當水塔中的水位下降到最低水位L L (27.00ft)(27.00ft)時水泵就自動啟動向水塔輸水時水泵就自動啟動向水塔輸水, , 當水塔水位當水塔水位達到限定最高水位達到限定最高水位H (35.00ft)H (35.00ft)時，水泵停止供水。水時，水泵停止供水。水泵供水期間無法測量水泵的供水量。水泵每天輸水一次泵供水期間無法測量水泵的供水量。水泵每天輸水一次或兩次或兩次, ,每次約二小時。每次約二小時。 已知某一小鎮的水塔是一高為已知某一小鎮的水塔是一高為40ft40ft，直徑為，直徑為57ft57ft的正圓柱，的正圓柱，下表記錄的是某一天水塔水位的真實數據下表記錄的是某一天水塔水位的真

33、實數據(1ft=0.3024(1ft=0.3024米米(m)(m)。試估計任何時刻。試估計任何時刻( (包括水泵正在輸水時間包括水泵正在輸水時間) )從水塔流出從水塔流出的水流量和一天的用水總量。的水流量和一天的用水總量。二、問題分析二、問題分析 流量是單位時間內流出水的體積，由于水塔是正圓柱流量是單位時間內流出水的體積，由于水塔是正圓柱形，截面面積是常數，所以，流量很容易根據水位相對形，截面面積是常數，所以，流量很容易根據水位相對時間的變化率算出。時間的變化率算出。 水泵不工作時段水泵不工作時段的用水量可以由測量記錄直接得到的用水量可以由測量記錄直接得到: :由表由表1 1中記錄的水位下降高

34、度乘以水塔的截面面積就是這中記錄的水位下降高度乘以水塔的截面面積就是這一時段的用水量。一時段的用水量。 問題的難點在于問題的難點在于: :如何估計水泵供水時段的流量。如何估計水泵供水時段的流量。 水泵供水時段的流量只能靠供水時段前后的流量經插水泵供水時段的流量只能靠供水時段前后的流量經插值或擬合得到。值或擬合得到。而作為用于插值或擬合的原始數據，我而作為用于插值或擬合的原始數據，我們希望水泵不工作時段的流量數據越準確越好。這個數們希望水泵不工作時段的流量數據越準確越好。這個數值可以用來檢驗數據插值或擬合結果的準確性。值可以用來檢驗數據插值或擬合結果的準確性。二、問題分析二、問題分析流量大體上可

35、由兩種方法計算流量大體上可由兩種方法計算: :(1)(1) 直接對表直接對表1 1中的用水量用數值插值方法算出各時段的中的用水量用數值插值方法算出各時段的流量，用它們擬合其它時刻或連續時間的流量流量，用它們擬合其它時刻或連續時間的流量; ;(2)(2) 先用表中數據擬合水位先用表中數據擬合水位-時間函數時間函數, ,求導數即可得到求導數即可得到連續時間的流量。連續時間的流量。 有了任何時刻的流量，就不難計算一天的總用水量。有了任何時刻的流量，就不難計算一天的總用水量。三、問題求解三、問題求解 為了表示方便，將所給表為了表示方便，將所給表1 1中的數據全部化為國際標準中的數據全部化為國際標準單位

36、，時間用小時單位，時間用小時(h),(h),高度用米高度用米(m)(m)。3.1 3.1 模型假設模型假設(1) (1) 流量只取決于水位差，與水位本身無關。流量只取決于水位差，與水位本身無關。(2) Torricelli(2) Torricelli定律：從小孔流出液體的流速正比于液面高定律：從小孔流出液體的流速正比于液面高度的平方根。度的平方根。 在假設出口的水位為零的前提下，題目給出水塔的最低在假設出口的水位為零的前提下，題目給出水塔的最低和最高水位分別是和最高水位分別是8.1648m(278.1648m(270.3024) 0.3024) 和和10.7352m(35.5010.7352m

37、(35.500.3024)0.3024)。 因為：因為：sqrt(10.7352/8.1648)=1.14671sqrt(10.7352/8.1648)=1.14671所以，可忽略水位對流速的影響。所以，可忽略水位對流速的影響。(3) (3) 流量可以看作是時間的連續函數。因此，為了計算簡單，流量可以看作是時間的連續函數。因此，為了計算簡單，不妨將流量定義成單位時間流出水的高度，即水位對時間不妨將流量定義成單位時間流出水的高度，即水位對時間變化率的絕對值。變化率的絕對值。(4) (4) 水塔截面面積為水塔截面面積為3.2 3.2 流量估計方法流量估計方法由表由表2 2給出的數據，用給出的數據，

38、用MTALABMTALAB軟件做出時間軟件做出時間水位散點圖水位散點圖1 1：為了計算水箱水流量與時間的關系，最簡單的方法就是：為了計算水箱水流量與時間的關系，最簡單的方法就是： 根據給出的散點數據圖，將數據分為三段，然后對每根據給出的散點數據圖，將數據分為三段，然后對每一段數據用數據插值或者數據擬合的方法得到流量與時一段數據用數據插值或者數據擬合的方法得到流量與時間的近似函數關系。間的近似函數關系。具體方法如下：具體方法如下：的平均流速采用公式的平均流速采用公式，相鄰數據中點，相鄰數據中點設某段數據為設某段數據為 ),( ,),(),( 1100nnyxyxyx中點流速中點流速=(=(左端點

39、的水位左端點的水位- -右端點的水位右端點的水位)/)/時間區間長度時間區間長度每段數據首尾點的流速用下面的公式計算：每段數據首尾點的流速用下面的公式計算：經過計算，得到時間與流速之間的關系數據表經過計算，得到時間與流速之間的關系數據表3 3 數據表數據表3 3對應的時間對應的時間-流速散點圖如下：流速散點圖如下：下面分別用數據插值和數據擬合兩種方法下面分別用數據插值和數據擬合兩種方法來估計水塔水流量。來估計水塔水流量。(1) (1) 數據插值法數據插值法 采用拉格朗日、分段線性和三次樣條三種數據采用拉格朗日、分段線性和三次樣條三種數據插值方法。插值方法。根據表根據表3 3的數據，對的數據，對

40、水泵不工作時段水泵不工作時段1,21,2通過采取數據插通過采取數據插值方法可以得到任意時刻的流速，從而知道任意時刻的值方法可以得到任意時刻的流速，從而知道任意時刻的流量。流量。對于對于水泵工作時段水泵工作時段1 1，應用前后時期的流速進行插值。，應用前后時期的流速進行插值。由于最后由于最后水泵工作時段水泵工作時段2 2的數據太少，我們將它與的數據太少，我們將它與水泵水泵不工作時段不工作時段3 3合并，一同進行數據插值處理（簡稱合并，一同進行數據插值處理（簡稱混合混合時段時段）。）。以第以第1 1未供水時段數據為例分別用三種方法算出流量函數和未供水時段數據為例分別用三種方法算出流量函數和用水量（

41、用水高度）。用水量（用水高度）。t=0,0.46,1.38,2.395,3.41,4.425,5.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.27;t0=0:0.1:8.97;lglr=lglrcz(t,v,t0); %拉格朗日插值，需要定義拉格朗日插值，需要定義m文件文件lglrjf=0.1*trapz(lglr) %梯形法算數值積分梯形法算數值積分fdxx=interp1(t,v,t0); %分段線性插值分段線性插值fdxxjf=0.1*trapz(fdxx

42、)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=145.6231lglrjf=145.6231，fdxxjf=147.1430 fdxxjf=147.1430 ，sancytjf=145.6870sancytjf=145.6870function y=lglrcz(x0,y0,x)n=length(x0);m=length(x);

43、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 2知，第知，第1 1未供水時段的總用水高度為未供水時段的總用水高度為146(=968-822)146(=968-822)。上述三種插值方法計算的結果與實際值上述三種插值方法計算的結果與實際值146146相比都比較接近。相比都比較接近。計算結果：計算結果：lglrjf=145.6231lglrjf=145.6231，fdxxjf=147.143

44、0 fdxxjf=147.1430 ，sancytjf=145.6870sancytjf=145.6870(2) (2) 數據擬合法數據擬合法1) 1) 擬合水位擬合水位-時間函數時間函數2) 2) 確定流量確定流量時間函數時間函數3) 3) 一天總用水量的估計一天總用水量的估計4) 4) 流量及總用水量的檢驗流量及總用水量的檢驗從表從表2 2中測量記錄看，一天有兩次供水時段和三次未供水中測量記錄看，一天有兩次供水時段和三次未供水時段，分別對第時段，分別對第1,21,2未供水時段的測量數據直接作多項式未供水時段的測量數據直接作多項式擬合，可得到水位函數。第擬合，可得到水位函數。第3 3未供水時

45、段數據過少未供水時段數據過少, ,擬合擬合效果不是很理想。效果不是很理想。應當注意應當注意: :根據多項式擬合的特點，此處擬合多項式的次根據多項式擬合的特點，此處擬合多項式的次數不宜過高，一般以數不宜過高，一般以3-63-6次為宜）。次為宜）。1) 1) 擬合水位擬合水位-時間函數時間函數t=0,0.92,1.84,2.95,3.87,4.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.

46、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,tp1);%變量變量x1中存放了以中存放了以0.1為步長算出的各個時刻為步長算出的各個時刻的水位高度。的水位高度。plot(tp1,x1)figure(2)c2=polyfit(t(11:21),h(11:21),3);tp2=10.9:0.1:20.9;x2=

47、-polyval(c2,tp2); plot(tp2,x2)第第1 1、2 2未供水時段時間與水位圖程序：未供水時段時間與水位圖程序：2) 2) 確定流量確定流量時間函數時間函數第第1,21,2未供水時段未供水時段第第1 1供水時段供水時段第第2 2供水時段第供水時段第3 3未供水時段混合時段未供水時段混合時段第第1,21,2未供水時段的流量可直接對水位函數求導得到，未供水時段的流量可直接對水位函數求導得到，用三次多項式擬合曲線，程序如下：用三次多項式擬合曲線，程序如下：c1=polyfit(t(1:10),h(1:10),3); c2=polyfit(t(11:21),h(11:21),3)

48、;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);% (為后面的計算準備數據為后面的計算準備數據)x23=-polyval(a2,tp2); x114=-polyval(a2,10.95:0.01:20.84);wgsysl2=100*trapz(tp2,x114);x24=-polyval(a

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