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

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

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

令

可以證明，對(duì)于n+1個(gè)不同結(jié)點(diǎn)，必存在唯一的次數(shù)不超過(guò)n的滿足條件的多項(xiàng)式，這個(gè)多項(xiàng)式稱為插值多項(xiàng)式，這種方法稱為n次多項(xiàng)式插值（或代數(shù)插值）。為了以后使用方便，先編制一個(gè)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

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

12345y-11.52.13.64.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=[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次多項(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ù)，沒(méi)有承接性。下面的牛頓插值法就避免了這個(gè)問(wèn)題。差商具有如下性質(zhì)：（1）m階差商是零階差商的線性組合；（2）差商與插值結(jié)點(diǎn)的次序無(wú)關(guān)；（3）若f(x)是m次多項(xiàng)式，則

f[x0,x1,…,xm]=0由差商公式，逐次代入得稱為牛頓插值公式，最后一項(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ù)。結(jié)點(diǎn)零階差商一階差商二階差商三階差商在做牛頓插值時(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)處滿足二階導(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是矩陣，則對(duì)y的每一列與x配對(duì)進(jìn)行插值；xb表示待求函數(shù)值的插值結(jié)點(diǎn)向量，可以缺?。弧甿ethod’是可選項(xiàng)，說(shuō)明插值使用的方法。缺省時(shí)為線性插值，也可選擇：nearest（最近插值），linear（線性），spline（三次樣條），cubic（三次插值）。命令返回值yb是插值曲線在xb處的縱坐標(biāo)值。2二維插值命令zb=interp2(x,y,z,xb,yb,'method')

根據(jù)同維數(shù)據(jù)向量x，y，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')zb=griddata(x,y,z,xb,yb,’method’)

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

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

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

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

分別用拉格朗日、分段線性、3次樣條進(jìn)行插值，并繪出插值多項(xiàng)式曲線圖。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次多項(xiàng)式曲線，在數(shù)據(jù)點(diǎn)之間產(chǎn)生較大的起伏波動(dòng)，與數(shù)據(jù)點(diǎn)的變化趨勢(shì)有明顯的偏離，這時(shí)曲線并不能很好地反映數(shù)據(jù)點(diǎn)的變化規(guī)律。而且進(jìn)一步實(shí)驗(yàn)發(fā)現(xiàn)，隨著分點(diǎn)的增加，Lagrange插值出現(xiàn)大的起伏波動(dòng)越明顯，這就是插值問(wèn)題中典型的“龍格現(xiàn)象”。針對(duì)高次多項(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=[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é)果分析分段線性插值有效地回避了插值問(wèn)題中的“龍格現(xiàn)象”，結(jié)果連線也大致描述了已知數(shù)據(jù)點(diǎn)的變化規(guī)律。但很明顯，由分段直線連接的插值曲線在節(jié)點(diǎn)處不光滑，不可導(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ù)點(diǎn)，既有效地回避了插值問(wèn)題中的“龍格現(xiàn)象”，又是連續(xù)光滑的，用此曲線來(lái)近似描述已知數(shù)據(jù)點(diǎn)的變化規(guī)律，能比較好地進(jìn)行數(shù)據(jù)點(diǎn)之間的預(yù)測(cè)分析和求值。例4.3某海域上頻繁地有各種噸位的船只經(jīng)過(guò)，為保證船只的航行安全，有關(guān)機(jī)構(gòu)在低潮時(shí)對(duì)水深進(jìn)行了測(cè)量，得到了如下的測(cè)量數(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為高度。請(qǐng)根據(jù)測(cè)量的數(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);三實(shí)驗(yàn)問(wèn)題求解1畫出河床觀測(cè)點(diǎn)的散點(diǎn)圖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ù)可以進(jìn)行分段線性插值,在此基礎(chǔ)上，利用梯形法求積分命令trapz來(lái)計(jì)算河床截面積，同時(shí)計(jì)算每一段連接線長(zhǎng)度之和來(lái)近似河床曲線長(zhǎng)度。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)從運(yùn)行結(jié)果得到：河床截面面積：S=337.15，河床曲線長(zhǎng)度：

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

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

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

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

問(wèn)題的難點(diǎn)在于:如何估計(jì)水泵供水時(shí)段的流量。

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

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

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

三、問(wèn)題求解

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

根據(jù)給出的散點(diǎn)數(shù)據(jù)圖，將數(shù)據(jù)分為三段，然后對(duì)每一段數(shù)據(jù)用數(shù)據(jù)插值或者數(shù)據(jù)擬合的方法得到流量與時(shí)間的近似函數(shù)關(guān)系。具體方法如下：中點(diǎn)流速=(左端點(diǎn)的水位-右端點(diǎn)的水位)/時(shí)間區(qū)間長(zhǎng)度每段數(shù)據(jù)首尾點(diǎn)的流速用下面的公式計(jì)算：經(jīng)過(guò)計(jì)算，得到時(shí)間與流速之間的關(guān)系數(shù)據(jù)表3數(shù)據(jù)表3對(duì)應(yīng)的時(shí)間--流速散點(diǎn)圖如下：下面分別用數(shù)據(jù)插值和數(shù)據(jù)擬合兩種方法來(lái)估計(jì)水塔水流量。(1)數(shù)據(jù)插值法采用拉格朗日、分段線性和三次樣條三種數(shù)據(jù)插值方法。根據(jù)表3的數(shù)據(jù)，對(duì)水泵不工作時(shí)段1,2通過(guò)采取數(shù)據(jù)插值方法可以得到任意時(shí)刻的流速，從而知道任意時(shí)刻的流量。對(duì)于水泵工作時(shí)段1，應(yīng)用前后時(shí)期的流速進(jìn)行插值。由于最后水泵工作時(shí)段2的數(shù)據(jù)太少，我們將它與水泵不工作時(shí)段3合并，一同進(jìn)行數(shù)據(jù)插值處理（簡(jiǎn)稱混合時(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];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')計(jì)算結(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未供水時(shí)段的總用水高度為146(=968-822)。上述三種插值方法計(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中測(cè)量記錄看，一天有兩次供水時(shí)段和三次未供水時(shí)段，分別對(duì)第1,2未供水時(shí)段的測(cè)量數(shù)據(jù)直接作多項(xiàng)式擬合，可得到水位函數(shù)。第3未供水時(shí)段數(shù)據(jù)過(guò)少,擬合效果不是很理想。應(yīng)當(dāng)注意:根據(jù)多項(xiàng)式擬合的特點(diǎn)，此處擬合多項(xiàng)式的次數(shù)不宜過(guò)高，一般以3-6次為宜）。1)擬合水位--時(shí)間函數(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為步長(zhǎng)算出的各個(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í)段的流量可直接對(duì)水位函數(shù)求導(dǎo)得到，用三次多項(xiàng)式擬合曲線，程序如下：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]);%(為后面的計(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,12.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)式重新擬合。對(duì)于第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*t