數(shù)學(xué)實驗：實驗14 水塔水流量估計模型與數(shù)據(jù)插值

MathematicsLaboratory阮小娥博士ExperimentsinMathematics數(shù)學(xué)實驗實驗14水塔水流量估計模型與數(shù)據(jù)插值2、學(xué)會用matlab軟件進行數(shù)據(jù)插值計算1、掌握四種經(jīng)典的插值方法：拉格朗日插值法、牛頓插值法、分段插值法、三次樣條插值法3、學(xué)會用數(shù)據(jù)插值、數(shù)據(jù)擬合方法建立數(shù)學(xué)模型并求解實驗問題一條100米寬的河道見圖1所示,為了測量其流量需要知道河道的截面積。為此從一端開始每隔5米測量出河床的深度如表1

試根據(jù)以上數(shù)據(jù)，估計出河道的截面積，進而在已知河流速度（設(shè)為1m/s）的情況下計算出流量。若沿河床鋪設(shè)一條光纜，試估計光纜的長度。一數(shù)據(jù)插值給定n個數(shù)據(jù)點稱P(x)為插值函數(shù)，也即求解一條嚴(yán)格通過各數(shù)據(jù)點的曲線，用它來進行分析研究和預(yù)測，這種方法常稱為數(shù)據(jù)插值法。對于這類問題，選取一條何種類型的曲線作為插值函數(shù)是求解的關(guān)鍵。P(x)不同，就得到不同的插值函數(shù)。1拉格朗日插值

令

可以證明，對于n+1個不同結(jié)點，必存在唯一的次數(shù)不超過n的滿足條件的多項式，這個多項式稱為插值多項式，這種方法稱為n次多項式插值（或代數(shù)插值）。為了以后使用方便，先編制一個Lagrange插值函數(shù)程序：functionp=lagrange(x,y)L=length(x);A=ones(L);%產(chǎn)生元素全是1的L階方陣forj=2:L

A(:,j)=A(:,j-1).*x';endX=inv(A)*y';fori=1:L

p(i)=X(L-i+1);end例4.1

已知觀測數(shù)據(jù)x

12345y-11.52.13.64.9求其插值多項式曲線。利用Matlab提供的計算以向量p為系數(shù)的多項式值的命令polyval

,可以求得多項式函數(shù)在任意一點的值。其使用格式為

y0=polyval(p,x0)其中p為多項式系數(shù)（從高次到底次排列）向量，x0為所求值的點或向量，y0為返回x0處的函數(shù)值，從而可以繪制多項式曲線。程序如下：x=[12345];y=[-11.52.13.64.9];plot(x,y,'b.','markersize',30)axis([15-15])gridholdonp=lagrange(x,y);t=1:0.1:5;u=polyval(p,t);plot(t,u,'r-','linewidth',3)從結(jié)果可以看到，所插值的4次多項式曲線較好地連接了5個數(shù)據(jù)點，從而可以用此多項式曲線作為這5個數(shù)據(jù)的一個近似變化。2牛頓插值法

拉格朗日插值法的最大缺點在于當(dāng)增加差值點時，需要重新計算多項式的系數(shù)，沒有承接性。下面的牛頓插值法就避免了這個問題。差商具有如下性質(zhì)：（1）m階差商是零階差商的線性組合；（2）差商與插值結(jié)點的次序無關(guān)；（3）若f(x)是m次多項式，則

f[x0,x1,…,xm]=0由差商公式，逐次代入得稱為牛頓插值公式，最后一項稱為牛頓插值余項，記為Rn(x)，余項前的多項式稱為插值多項式，記為Pn(x)。牛頓插值多項式具有以下特點：（1）在插值結(jié)點處與拉格朗日插值一樣，誤差為零；（2）多項式k次項的系數(shù)是f(x)的k階差商；（3）增加插值節(jié)點時，只增加最后一項，不必像拉格朗日插值公式那樣需要重新計算系數(shù)。結(jié)點零階差商一階差商二階差商三階差商在做牛頓插值時，一般先做出差商表，然后套用公式。3樣條插值法

特別地，如果m=1，則稱為分段線性插值，即在每一個子區(qū)間上S(x)是線性函數(shù)，在整個區(qū)間上是分段線性函數(shù)。這種情況一般很難滿足實際要求，通常使用最多的是3次樣條插值函數(shù)，即S(x)在每一個子區(qū)間上是三次多項式函數(shù)，而且在每個結(jié)點處滿足二階導(dǎo)數(shù)連續(xù)和相關(guān)的邊界條件。二MATLAB軟件實現(xiàn)數(shù)據(jù)插值1一維插值命令

yb=interp1(x,y,xb,'method')x,y是同維數(shù)據(jù)向量，表示插值結(jié)點的橫坐標(biāo)和縱坐標(biāo)。若x是向量，y是矩陣，則對y的每一列與x配對進行插值；xb表示待求函數(shù)值的插值結(jié)點向量，可以缺??；‘method’是可選項，說明插值使用的方法。缺省時為線性插值，也可選擇：nearest（最近插值），linear（線性），spline（三次樣條），cubic（三次插值）。命令返回值yb是插值曲線在xb處的縱坐標(biāo)值。2二維插值命令zb=interp2(x,y,z,xb,yb,'method')

根據(jù)同維數(shù)據(jù)向量x，y，z，按照指定的方法做插值運算，然后返回插值函數(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')zb=griddata(x,y,z,xb,yb,’method’)

x,y,z是同維數(shù)據(jù)向量，表示插值結(jié)點的橫坐標(biāo)、縱坐標(biāo)和豎坐標(biāo)。xb，yb表示待求函數(shù)值的插值結(jié)點向量，通常由meshgrid生成?！甿ethod’是可選項，說明插值使用的方法。缺省時為線性插值，也可選擇：nearest（最近插值），linear（線性），cubic（三次插值），v4。

v4是一種插值算法，沒有具體的名字，原文稱為“MATLAB4griddatamethod”，是一種很圓滑的差值算法，效果不錯。

命令返回值zb是插值曲面在(xb,yb)處的豎坐標(biāo)值。5用griddata命令插值生成三維曲面

interp2和griddata都是二維插值，兩者的區(qū)別是：interp2的插值數(shù)據(jù)必須是矩形域，即已知數(shù)據(jù)點(x,y)組成規(guī)則的矩陣，或稱之為柵格，可使用meshgid生成。而griddata函數(shù)的已知數(shù)據(jù)點（X，Y）不要求規(guī)則排列，特別是對試驗中隨機沒有規(guī)律采取的數(shù)據(jù)進行插值具有很好的效果。例4.2已知觀測數(shù)據(jù)

分別用拉格朗日、分段線性、3次樣條進行插值，并繪出插值多項式曲線圖。x00.10.20.30.40.50.60.70.80.91y-0.4471.9783.286.167.087.347.669.569.489.311.2(1)拉格朗日插值法x=0:0.1:1；y=[-0.4471.9783.286.167.087.347.669.569.489.311.2];plot(x,y,'b.','markersize',30)axis([01-216])gridholdonp=lagrange(x,y);t=0:0.01:1;u=polyval(p,t);plot(t,u,'r-','linewidth',3)

結(jié)果顯示，所插值出的10次多項式曲線，在數(shù)據(jù)點之間產(chǎn)生較大的起伏波動，與數(shù)據(jù)點的變化趨勢有明顯的偏離，這時曲線并不能很好地反映數(shù)據(jù)點的變化規(guī)律。而且進一步實驗發(fā)現(xiàn)，隨著分點的增加，Lagrange插值出現(xiàn)大的起伏波動越明顯，這就是插值問題中典型的“龍格現(xiàn)象”。針對高次多項式插值時容易發(fā)生“龍格現(xiàn)象”，在實際插值時，常常采用分段低次插值方法，即在相鄰兩個數(shù)據(jù)點構(gòu)成的子區(qū)間上分別進行低次插值,整個區(qū)間上的插值函數(shù)將是一個分段的多項式函數(shù)。(2)分段線性插值x=[00.10.20.30.40.50.60.70.80.91];y=[-0.4471.9783.286.167.087.347.669.569.489.311.2];plot(x,y,‘b.’,‘markersize’,30)；axis([01-115])；gridholdont=0:0.01:1;u=interp1(x,y,t);plot(t,u,'r-','linewidth',3)結(jié)果分析分段線性插值有效地回避了插值問題中的“龍格現(xiàn)象”，結(jié)果連線也大致描述了已知數(shù)據(jù)點的變化規(guī)律。但很明顯，由分段直線連接的插值曲線在節(jié)點處不光滑，不可導(dǎo)。(3)3次樣條插值x=[00.10.20.30.40.50.60.70.80.91];y=[-0.4471.9783.286.167.087.347.669.569.489.311.2];plot(x,y,'b.','markersize',30)axis([01-116])gridholdont=0:0.01:1;u=spline(x,y,t);plot(t,u,'r-','linewidth',3)結(jié)果分析從圖可以看出，樣條插值所得的曲線比較好地連接了已知的數(shù)據(jù)點，既有效地回避了插值問題中的“龍格現(xiàn)象”，又是連續(xù)光滑的，用此曲線來近似描述已知數(shù)據(jù)點的變化規(guī)律，能比較好地進行數(shù)據(jù)點之間的預(yù)測分析和求值。例4.3某海域上頻繁地有各種噸位的船只經(jīng)過，為保證船只的航行安全，有關(guān)機構(gòu)在低潮時對水深進行了測量，得到了如下的測量數(shù)據(jù)：x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];

y=[7.5141.523.0147.022.5137.585.560.5121.03.056.5116.584.043.5];

Z=[48686889988949];

x,y為平面坐標(biāo)，z為高度。請根據(jù)測量的數(shù)據(jù)描述該海域的地貌。程序如下：clc;x=[129.0140.0103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5];y=[7.5141.523.0147.022.5137.585.560.5121.03.056.5116.584.043.5];z=[48686889988949];h=-z;xi=linspace(min(x),max(x),100);yi=linspace(min(y),max(y),100);[X,Y]=meshgrid(xi,yi);H=griddata(x,y,h,X,Y,‘v4');mesh(X,Y,H);三實驗問題求解1畫出河床觀測點的散點圖x=0:5:100;y=[02.412.962.152.653.124.235.126.215.684.22…3.913.262.852.353.023.634.123.462.080];axis([0100010])gridy1=10-y;plot(x,y1,'b.','markersize',30)2利用分段線性插值繪制河床曲線根據(jù)已給的數(shù)據(jù)可以進行分段線性插值,在此基礎(chǔ)上，利用梯形法求積分命令trapz來計算河床截面積，同時計算每一段連接線長度之和來近似河床曲線長度。x=0:5:100;y=[02.412.962.152.653.124.235.126.215.684.22…3.913.262.852.353.023.634.123.462.080];y1=10-y;plot(x,y1,'b.','markersize',30);axis([0100010])gridholdont=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=%.2f\n',S,L)從運行結(jié)果得到：河床截面面積：S=337.15，河床曲線長度：

L=102.093利用樣條插值計算截面積與曲線長度另一方面,為了提高河床曲線的模擬精度，可以根據(jù)已給的數(shù)據(jù)進行3次樣條插值,在此基礎(chǔ)上，利用梯形法求積分命令trapz來計算河床截面積，同時對樣條曲線加密，計算每一段連接線長度之和來近似河床曲線長度。x=0:5:100;y=[02.412.962.152.653.124.235.126.215.684.22...3.913.262.852.353.023.634.123.462.080];y1=10-y;plot(x,y1,'b.','markersize',30)axis([0100210])gridholdont=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=%.2f\n',S,L)河床截面面積：S=339.50，河床曲線長度：L=102.23.數(shù)據(jù)插值模型試驗：水塔水流量估計一．實驗問題某一地區(qū)的用水管理機構(gòu)要求各社區(qū)提供每小時以加侖計的用水量以及每天所用水的總量。由于許多社區(qū)沒有測量流入或流出水塔的水量裝置，所以他們只能通過測量水塔每小時的水位高度來代替每小時的用水量(誤差不超過0.05)。當(dāng)水塔中的水位下降到最低水位L(27.00ft)時水泵就自動啟動向水塔輸水,當(dāng)水塔水位達(dá)到限定最高水位H(35.00ft)時，水泵停止供水。水泵供水期間無法測量水泵的供水量。水泵每天輸水一次或兩次,每次約二小時。

已知某一小鎮(zhèn)的水塔是一高為40ft，直徑為57ft的正圓柱，下表記錄的是某一天水塔水位的真實數(shù)據(jù)(1ft=0.3024米(m))。試估計任何時刻(包括水泵正在輸水時間)從水塔流出的水流量和一天的用水總量。二、問題分析

流量是單位時間內(nèi)流出水的體積，由于水塔是正圓柱形，截面面積是常數(shù)，所以，流量很容易根據(jù)水位相對時間的變化率算出。

水泵不工作時段的用水量可以由測量記錄直接得到:由表1中記錄的水位下降高度乘以水塔的截面面積就是這一時段的用水量。

問題的難點在于:如何估計水泵供水時段的流量。

水泵供水時段的流量只能靠供水時段前后的流量經(jīng)插值或擬合得到。而作為用于插值或擬合的原始數(shù)據(jù)，我們希望水泵不工作時段的流量數(shù)據(jù)越準(zhǔn)確越好。這個數(shù)值可以用來檢驗數(shù)據(jù)插值或擬合結(jié)果的準(zhǔn)確性。二、問題分析流量大體上可由兩種方法計算:

直接對表1中的用水量用數(shù)值插值方法算出各時段的流量，用它們擬合其它時刻或連續(xù)時間的流量;

先用表中數(shù)據(jù)擬合水位--時間函數(shù),求導(dǎo)數(shù)即可得到連續(xù)時間的流量。有了任何時刻的流量，就不難計算一天的總用水量。

三、問題求解

為了表示方便，將所給表1中的數(shù)據(jù)全部化為國際標(biāo)準(zhǔn)單位，時間用小時(h),高度用米(m)。3.1模型假設(shè)(1)流量只取決于水位差，與水位本身無關(guān)。(2)Torricelli定律：從小孔流出液體的流速正比于液面高度的平方根。在假設(shè)出口的水位為零的前提下，題目給出水塔的最低和最高水位分別是8.1648m(27×0.3024)和10.7352m(35.50×0.3024)。因為：sqrt(10.7352/8.1648)=1.1467≈1所以，可忽略水位對流速的影響。(3)流量可以看作是時間的連續(xù)函數(shù)。因此，為了計算簡單，不妨將流量定義成單位時間流出水的高度，即水位對時間變化率的絕對值。(4)水塔截面面積為3.2流量估計方法由表2給出的數(shù)據(jù)，用MTALAB軟件做出時間—水位散點圖1：為了計算水箱水流量與時間的關(guān)系，最簡單的方法就是：

根據(jù)給出的散點數(shù)據(jù)圖，將數(shù)據(jù)分為三段，然后對每一段數(shù)據(jù)用數(shù)據(jù)插值或者數(shù)據(jù)擬合的方法得到流量與時間的近似函數(shù)關(guān)系。具體方法如下：中點流速=(左端點的水位-右端點的水位)/時間區(qū)間長度每段數(shù)據(jù)首尾點的流速用下面的公式計算：經(jīng)過計算，得到時間與流速之間的關(guān)系數(shù)據(jù)表3數(shù)據(jù)表3對應(yīng)的時間--流速散點圖如下：下面分別用數(shù)據(jù)插值和數(shù)據(jù)擬合兩種方法來估計水塔水流量。(1)數(shù)據(jù)插值法采用拉格朗日、分段線性和三次樣條三種數(shù)據(jù)插值方法。根據(jù)表3的數(shù)據(jù)，對水泵不工作時段1,2通過采取數(shù)據(jù)插值方法可以得到任意時刻的流速，從而知道任意時刻的流量。對于水泵工作時段1，應(yīng)用前后時期的流速進行插值。由于最后水泵工作時段2的數(shù)據(jù)太少，我們將它與水泵不工作時段3合并，一同進行數(shù)據(jù)插值處理（簡稱混合時段）。以第1未供水時段數(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];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,'b')gtext('lglr')gtext('fdxx')gtext('scyt')計算結(jié)果：lglrjf=145.6231，fdxxjf=147.1430，sancytjf=145.6870functiony=lglrcz(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;forj=1:nifj~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;end

y(i)=s;end由表2知，第1未供水時段的總用水高度為146(=968-822)。上述三種插值方法計算的結(jié)果與實際值146相比都比較接近。計算結(jié)果：lglrjf=145.6231，fdxxjf=147.1430，sancytjf=145.6870(2)數(shù)據(jù)擬合法1)擬合水位--時間函數(shù)2)確定流量—時間函數(shù)3)一天總用水量的估計4)流量及總用水量的檢驗從表2中測量記錄看，一天有兩次供水時段和三次未供水時段，分別對第1,2未供水時段的測量數(shù)據(jù)直接作多項式擬合，可得到水位函數(shù)。第3未供水時段數(shù)據(jù)過少,擬合效果不是很理想。應(yīng)當(dāng)注意:根據(jù)多項式擬合的特點，此處擬合多項式的次數(shù)不宜過高，一般以3-6次為宜）。1)擬合水位--時間函數(shù)t=[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.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.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];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=-polyval(c2,tp2);plot(tp2,x2)第1、2未供水時段時間與水位圖程序：2)確定流量—時間函數(shù)第1,2未供水時段第1供水時段第2供水時段＋第3未供水時段＝混合時段第1,2未供水時段的流量可直接對水位函數(shù)求導(dǎo)得到，用三次多項式擬合曲線，程序如下：c1=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);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);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)subplot(1,2,2)plot(tp2,x23*100)上下曲線比較發(fā)現(xiàn)，用三次多項式擬合效果不是很好，故改用五次多項式重新擬合。對于第1供水時段的流量,則用前后時期的流量進行擬合得到。為使流量函數(shù)在t=8.97和t=10.93連續(xù)，只取4個點，用3次多項式擬合得到第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,nhsj);gssj1=8.97:0.01:10.95;gs1=polyval(nhjg,[8.97:0.01:10.95]);gsysl1=100*t