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

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

令

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

已知觀測數據x

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

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

y0=polyval(p,x0)其中p為多項式系數（從高次到底次排列）向量，x0為所求值的點或向量，y0為返回x0處的函數值，從而可以繪制多項式曲線。程序如下：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)從結果可以看到，所插值的4次多項式曲線較好地連接了5個數據點，從而可以用此多項式曲線作為這5個數據的一個近似變化。2牛頓插值法

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

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

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

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

根據同維數據向量x，y，z，按照指定的方法做插值運算，然后返回插值函數的豎坐標。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是同維數據向量，表示插值結點的橫坐標、縱坐標和豎坐標。xb，yb表示待求函數值的插值結點向量，通常由meshgrid生成?！甿ethod’是可選項，說明插值使用的方法。缺省時為線性插值，也可選擇：nearest（最近插值），linear（線性），cubic（三次插值），v4。

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

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

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

分別用拉格朗日、分段線性、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)

結果顯示，所插值出的10次多項式曲線，在數據點之間產生較大的起伏波動，與數據點的變化趨勢有明顯的偏離，這時曲線并不能很好地反映數據點的變化規(guī)律。而且進一步實驗發(fā)現(xiàn)，隨著分點的增加，Lagrange插值出現(xiàn)大的起伏波動越明顯，這就是插值問題中典型的“龍格現(xiàn)象”。針對高次多項式插值時容易發(fā)生“龍格現(xiàn)象”，在實際插值時，常常采用分段低次插值方法，即在相鄰兩個數據點構成的子區(qū)間上分別進行低次插值,整個區(qū)間上的插值函數將是一個分段的多項式函數。(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)結果分析分段線性插值有效地回避了插值問題中的“龍格現(xiàn)象”，結果連線也大致描述了已知數據點的變化規(guī)律。但很明顯，由分段直線連接的插值曲線在節(jié)點處不光滑，不可導。(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)結果分析從圖可以看出，樣條插值所得的曲線比較好地連接了已知的數據點，既有效地回避了插值問題中的“龍格現(xiàn)象”，又是連續(xù)光滑的，用此曲線來近似描述已知數據點的變化規(guī)律，能比較好地進行數據點之間的預測分析和求值。例4.3某海域上頻繁地有各種噸位的船只經過，為保證船只的航行安全，有關機構在低潮時對水深進行了測量，得到了如下的測量數據：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為平面坐標，z為高度。請根據測量的數據描述該海域的地貌。程序如下：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利用分段線性插值繪制河床曲線根據已給的數據可以進行分段線性插值,在此基礎上，利用梯形法求積分命令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)從運行結果得到：河床截面面積：S=337.15，河床曲線長度：

L=102.093利用樣條插值計算截面積與曲線長度另一方面,為了提高河床曲線的模擬精度，可以根據已給的數據進行3次樣條插值,在此基礎上，利用梯形法求積分命令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.數據插值模型試驗：水塔水流量估計一．實驗問題某一地區(qū)的用水管理機構要求各社區(qū)提供每小時以加侖計的用水量以及每天所用水的總量。由于許多社區(qū)沒有測量流入或流出水塔的水量裝置，所以他們只能通過測量水塔每小時的水位高度來代替每小時的用水量(誤差不超過0.05)。當水塔中的水位下降到最低水位L(27.00ft)時水泵就自動啟動向水塔輸水,當水塔水位達到限定最高水位H(35.00ft)時，水泵停止供水。水泵供水期間無法測量水泵的供水量。水泵每天輸水一次或兩次,每次約二小時。

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

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

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

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

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

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

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

三、問題求解

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

根據給出的散點數據圖，將數據分為三段，然后對每一段數據用數據插值或者數據擬合的方法得到流量與時間的近似函數關系。具體方法如下：中點流速=(左端點的水位-右端點的水位)/時間區(qū)間長度每段數據首尾點的流速用下面的公式計算：經過計算，得到時間與流速之間的關系數據表3數據表3對應的時間--流速散點圖如下：下面分別用數據插值和數據擬合兩種方法來估計水塔水流量。(1)數據插值法采用拉格朗日、分段線性和三次樣條三種數據插值方法。根據表3的數據，對水泵不工作時段1,2通過采取數據插值方法可以得到任意時刻的流速，從而知道任意時刻的流量。對于水泵工作時段1，應用前后時期的流速進行插值。由于最后水泵工作時段2的數據太少，我們將它與水泵不工作時段3合并，一同進行數據插值處理（簡稱混合時段）。以第1未供水時段數據為例分別用三種方法算出流量函數和用水量（用水高度）。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)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.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)。上述三種插值方法計算的結果與實際值146相比都比較接近。計算結果：lglrjf=145.6231，fdxxjf=147.1430，sancytjf=145.6870(2)數據擬合法1)擬合水位--時間函數2)確定流量—時間函數3)一天總用水量的估計4)流量及總用水量的檢驗從表2中測量記錄看，一天有兩次供水時段和三次未供水時段，分別對第1,2未供水時段的測量數據直接作多項式擬合，可得到水位函數。第3未供水時段數據過少,擬合效果不是很理想。應當注意:根據多項式擬合的特點，此處擬合多項式的次數不宜過高，一般以3-6次為宜）。1)擬合水位--時間函數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)確定流量—時間函數第1,2未供水時段第1供水時段第2供水時段＋第3未供水時段＝混合時段第1,2未供水時段的流量可直接對水位函數求導得到，用三次多項式擬合曲線，程序如下：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]);%(為后面的計算準備數據)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]);%(為后面的計算準備數據)x25=-polyval(a2,[19.96,20.84]);%(為后面的計算準備數據)subplot(1,2,1)plot(tp1,x13*100)subplot(1,2,2)plot(tp2,x23*100)上下曲線比較發(fā)現(xiàn)，用三次多項式擬合效果不是很好，故改用五次多項式重新擬合。對于第1供水時段的流量,則用前后時期的流量進行擬合得到。為使流量函數在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