教學(xué)項(xiàng)目 數(shù)值法模型

數(shù)學(xué)建模與實(shí)驗(yàn)嚴(yán)可頌主講教學(xué)項(xiàng)目四數(shù)值分析法模型4.1插值法建模

拉格朗日插值分段線(xiàn)性插值三次樣條插值一、插值的定義二、插值的方法三、用Matlab解插值問(wèn)題已知n+1個(gè)節(jié)點(diǎn)其中互不相同，不妨設(shè)求任一插值點(diǎn)處的插值節(jié)點(diǎn)可視為由產(chǎn)生，表達(dá)式復(fù)雜，或無(wú)封閉形式，或未知。4.1.1插值法的定義

構(gòu)造一個(gè)(相對(duì)簡(jiǎn)單的)函數(shù)通過(guò)全部節(jié)點(diǎn),即再用計(jì)算插值，即4.1.2插值的方法

已知函數(shù)f(x)在n+1個(gè)點(diǎn)x0,x1,…,xn處的函數(shù)值為y0,y1,…,yn

。求一n次多項(xiàng)式函數(shù)Pn(x)，使其滿(mǎn)足：拉格朗日(Lagrange)插值稱(chēng)為拉格朗日插值基函數(shù)。解決此問(wèn)題的拉格朗日插值多項(xiàng)式公式如下其中Li(x)為n次多項(xiàng)式：特殊情況兩點(diǎn)一次(線(xiàn)性)插值多項(xiàng)式:三點(diǎn)二次(拋物)插值多項(xiàng)式:采用拉格朗日多項(xiàng)式插值：選取不同插值節(jié)點(diǎn)個(gè)數(shù)n+1，其中n為插值多項(xiàng)式的次數(shù)，當(dāng)n取10時(shí)，繪出插值結(jié)果圖形.例分段線(xiàn)性插值計(jì)算量與n無(wú)關(guān);n越大，誤差越小.xjxj-1xj+1x0xnxoy比分段線(xiàn)線(xiàn)性插值值更光滑滑。xyxi-1xiab在數(shù)學(xué)上上，光滑滑程度的的定量描描述是：：函數(shù)(曲線(xiàn))的k階導(dǎo)數(shù)存存在且連連續(xù)，則則稱(chēng)該曲線(xiàn)具具有k階光滑性性。光滑性的的階次越越高，則則越光滑滑。是否否存在較較低次的的分段多多項(xiàng)式達(dá)達(dá)到較高高階光滑滑性的方方法？三三次樣條條插值就就是一個(gè)個(gè)很好的的例子。。三次樣條條插值三次樣條條插值g(x)為被插值函函數(shù)。用MATLAB作插值計(jì)計(jì)算一維插值值函數(shù)：：yi=interp1(x，y，xi，'method')插值方法被插值點(diǎn)插值節(jié)點(diǎn)xi處的插值結(jié)果‘nearest’：最鄰近近插值‘linear’’：線(xiàn)性插值值；‘spline’’：三次次樣條插插值；‘cubic’：立方方插值。。缺省時(shí)：：分分段線(xiàn)性性插值。。注意：所所有的插插值方法法都要求求x是單調(diào)的的，并且且xi不能夠超超過(guò)x的范圍。。例：在1-12的11小時(shí)內(nèi)，，每隔1小時(shí)測(cè)量量一次溫溫度，測(cè)測(cè)得的溫溫度依次次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計(jì)計(jì)每隔1/10小時(shí)的溫溫度值。。hours=1:12;h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接輸出出數(shù)據(jù)將將是很多多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius’)xy實(shí)驗(yàn)習(xí)題題已知飛機(jī)機(jī)下輪廓廓線(xiàn)上數(shù)數(shù)據(jù)如下下，求x每改變0.1時(shí)的y值。機(jī)翼下輪廓線(xiàn)4.2擬合法簡(jiǎn)介

2.擬合的基基本原理理1.擬合問(wèn)題題引例擬合問(wèn)問(wèn)題題引例例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求600C時(shí)的電阻阻R。設(shè)R=at+ba,b為待定系系數(shù)擬合問(wèn)問(wèn)題題引例例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射300mg)求血藥濃濃度隨時(shí)時(shí)間的變變化規(guī)律律c(t).作半對(duì)數(shù)數(shù)坐標(biāo)系系(semilogy)下的圖形形曲線(xiàn)擬擬合合問(wèn)題題的的提法法已知一組組（二維維）數(shù)據(jù)據(jù)，即平平面上n個(gè)點(diǎn)（xi,yi)i=1,…n,尋求一個(gè)個(gè)函數(shù)（（曲線(xiàn)））y=f(x),使f(x)在某種準(zhǔn)準(zhǔn)則下與與所有數(shù)數(shù)據(jù)點(diǎn)最最為接近近，即曲曲線(xiàn)擬合合得最好好。+++++++++xyy=f(x)(xi,yi)ii為點(diǎn)（xi,yi)與曲線(xiàn)y=f(x)的距離擬合與插值的的關(guān)系函數(shù)插值與曲曲線(xiàn)擬合都是是要根據(jù)一組組數(shù)據(jù)構(gòu)造一一個(gè)函數(shù)作為為近似，由于于近似的要求求不同，二者者的數(shù)學(xué)方法法上是完全不不同的。實(shí)例：下面數(shù)據(jù)是某某次實(shí)驗(yàn)所得得，希望得到到X和f之間的關(guān)系？？問(wèn)題：給定一批數(shù)據(jù)據(jù)點(diǎn)，需確定定滿(mǎn)足特定要要求的曲線(xiàn)或或曲面解決方案：若不要求曲線(xiàn)線(xiàn)（面）通過(guò)過(guò)所有數(shù)據(jù)點(diǎn)點(diǎn)，而是要求求它反映對(duì)象象整體的變化化趨勢(shì)，這就就是數(shù)據(jù)擬合，又稱(chēng)曲線(xiàn)擬擬合或曲面擬擬合。若要求所求曲曲線(xiàn)（面）通通過(guò)所給所有有數(shù)據(jù)點(diǎn)，就就是插值問(wèn)題；最臨近插值、、線(xiàn)性插值、、樣條插值與與曲線(xiàn)擬合結(jié)結(jié)果：曲線(xiàn)擬合問(wèn)題題最常用的解解法——線(xiàn)性最小二乘乘法的基本思思路第一步:先選定一組函函數(shù)r1(x),r2(x),……rm(x),m<n,令f(x)=a1r1(x)+a2r2(x)+……+amrm(x)（1）其中a1,a2,…am為待定系數(shù)。。第二步:確定a1,a2,…am的準(zhǔn)則（最小小二乘準(zhǔn)則））：使n個(gè)點(diǎn)（xi,yi)與曲線(xiàn)y=f(x)的距離i的平方和最小小。記

問(wèn)題歸結(jié)為，，求a1,a2,…am使J(a1,a2,…am)最小。線(xiàn)性最小二乘乘法的求解：：預(yù)備知識(shí)超定方程組：方程個(gè)數(shù)大大于未知量個(gè)個(gè)數(shù)的方程組組即Ra=y其中超定方程一般般是不存在解解的矛盾方程程組。如果有向量a使得達(dá)到最小，則稱(chēng)a為上述超定方程的最小二乘解。線(xiàn)性最小二乘乘法的求解定理：當(dāng)RTR可逆時(shí)，超定定方程組（3）存在最小二二乘解，且即即為方程組RTRa=RTy的解：a=(RTR)-1RTy所以，曲線(xiàn)擬擬合的最小二二乘法要解決決的問(wèn)題，實(shí)實(shí)際上就是求求以下超定方方程組的最小小二乘解的問(wèn)問(wèn)題。其中Ra=y（3）線(xiàn)性最小二乘乘擬合f(x)=a1r1(x)+……+amrm(x)中函數(shù){r1(x),……rm(x)}的選取1.通過(guò)機(jī)理分析析建立數(shù)學(xué)模模型來(lái)確定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx2.將數(shù)據(jù)(xi,yi)i=1,…n作圖，通過(guò)直直觀(guān)判斷確定定f(x)：用MATLAB解擬合問(wèn)題1、線(xiàn)性最小二二乘擬合2、非線(xiàn)性最小小二乘擬合用MATLAB作線(xiàn)性最小二二乘擬合1.作多項(xiàng)式f(x)=a1xm+…+amx+am+1擬合,可利用已有程程序:a=polyfit(x,y,m)2.對(duì)超定方程組可得最小二乘意義下的解。，用3.多項(xiàng)式在x處的值y可用以下命令令計(jì)算：y=polyval（a，x）輸入同長(zhǎng)度的數(shù)組X，Y擬合多項(xiàng)式次數(shù)即要求出二次多項(xiàng)式:中的使得:例對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合1）輸入以下命命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];R=[(x.^2)'x'ones(11,1)]；A=R\y'解法1．用解超定方程程的方法2）計(jì)算結(jié)果：Ａ1）輸入以下命命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出數(shù)據(jù)點(diǎn)和和擬合曲線(xiàn)的的圖形2）計(jì)算結(jié)果：：Ａ解法2．用多項(xiàng)式擬合合的命令1.lsqcurvefit已知數(shù)據(jù)點(diǎn)：xdata=（xdata1，xdata2，…，xdatan），ydata=（ydata1，ydata2，…，ydatan）用MATLAB作非線(xiàn)性最小小二乘擬合Matlab的提供了兩個(gè)個(gè)求非線(xiàn)性最最小二乘擬合合的函數(shù)：lsqcurvefit和lsqnonlin。兩個(gè)命令都都要先建立M-文件fun.m，在其中定義義函數(shù)f(x)，但兩者定義義f(x)的方式是不同同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x，xdata1），…，F(xiàn)（x，xdatan））T中的參變量x(向量),使得輸入格式為:(1)x=lsqcurvefit(‘fun’’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一個(gè)事先建立的定義函數(shù)F(x,xdata)

的M-文件,自變量為x和xdata說(shuō)明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點(diǎn)選項(xiàng)見(jiàn)無(wú)約束優(yōu)化

lsqnonlin用以求含參量x（向量）的向量值函數(shù)

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知數(shù)據(jù)點(diǎn)：：xdata=（xdata1，xdata2，…，xdatan）ydata=（ydata1，ydata2，…，ydatan）輸入格式為：：1）x=lsqnonlin（‘fun’，x0）；2）x=lsqnonlin（‘fun’，x0，options）；3）x=lsqnonlin（‘fun’，x0，options，‘grad’）；4）[x，options]=lsqnonlin（‘fun’，x0，…）；5）[x，options，funval]=lsqnonlin（‘fun’，x0，…）；說(shuō)明：x=lsqnonlin（‘fun’，x0，options）；fun是一個(gè)事先建立的定義函數(shù)f(x)的M-文件，自變量為x迭代初值選項(xiàng)見(jiàn)無(wú)約束優(yōu)化

例2用下面一組數(shù)據(jù)擬合中的參數(shù)a，b，k該問(wèn)題即解最最優(yōu)化問(wèn)題：：1）編寫(xiě)M-文件curvefun1.mfunctionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）輸入命令tdata=100:100:1000cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];x0=[0.2,0.05,0.05];x=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)

F(x，tdata)=，x=(a，b，k)解法1.用命令令lsqcurvefit3）運(yùn)算算結(jié)果果為：4）結(jié)論論：a=0.0063,b=-0.0034,k=0.2542

解法2

用命令lsqnonlin

f(x)=F(x,tdata,ctada)=x=（a，b，k）1）編寫(xiě)M-文件curvefun2.mfunctionf=curvefun2(x)tdata=100:100:1000;cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];f=x(1)+x(2)*exp(-0.02*x(3)*tdata)-cdata2）輸入入命令令:x0=[0.2,0.05,0.05];x=lsqnonlin('curvefun2',x0)f=curvefun2(x)3）運(yùn)算算結(jié)果果為可以看看出,兩個(gè)命命令的的計(jì)算算結(jié)果果是相相同的的.4）結(jié)論論：即擬合合得a=0.0063b=-0.0034k=0.2542MATLAB解應(yīng)用用問(wèn)題題實(shí)例例1、電阻阻問(wèn)題題2、給藥藥方案案問(wèn)題題*3、水塔塔流量量估計(jì)計(jì)問(wèn)題題電阻問(wèn)問(wèn)題溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032例.由數(shù)據(jù)擬合R=a1t+a2方法1.用命令polyfit(x,y,m)得到a1=3.3940,a2=702.4918方法2.直接用結(jié)果相相同。。一室模模型：將整整個(gè)機(jī)機(jī)體看看作一一個(gè)房房室，，稱(chēng)中心室室，室內(nèi)內(nèi)血藥藥濃度度是均均勻的的?？炜焖凫o靜脈注注射后后，濃濃度立立即上上升；；然后后迅速速下降降。當(dāng)當(dāng)濃度度太低低時(shí)，，達(dá)不不到預(yù)預(yù)期的的治療療效果果；當(dāng)當(dāng)濃度度太高高，又又可能能導(dǎo)致致藥物物中毒毒或副副作用用太強(qiáng)強(qiáng)。臨臨床上上，每每種藥藥物有有一個(gè)個(gè)最小小有效效濃度度c1和一個(gè)個(gè)最大大有效效濃度度c2。設(shè)計(jì)計(jì)給藥藥方案案時(shí)，，要使使血藥藥濃度度保保持在在c1~c2之間。。本題題設(shè)c1=10，c2=25(ug/ml).擬合問(wèn)題實(shí)例2給藥方案——一種新新藥用用于臨臨床之之前，，必須須設(shè)計(jì)計(jì)給藥藥方案案.藥物進(jìn)進(jìn)入機(jī)機(jī)體后后血液液輸送送到全全身，，在這這個(gè)過(guò)過(guò)程中中不斷斷地被被吸收收、分分布、、代謝謝，最最終排排出體體外，，藥物物在血血液中中的濃濃度，，即單單位體體積血血液中中的藥藥物含含量，，稱(chēng)為為血藥濃濃度。。

在實(shí)驗(yàn)方面,對(duì)某人用快速靜脈注射方式一次注入該藥物300mg后,在一定時(shí)刻t(小時(shí))采集血藥,測(cè)得血藥濃度c(ug/ml)如下表:

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01要設(shè)計(jì)計(jì)給藥藥方案案,必須知知道給給藥后后血藥藥濃度度隨時(shí)時(shí)間變變化的的規(guī)律律。從從實(shí)驗(yàn)驗(yàn)和理理論兩兩方面面著手手：給藥方方案1.在快速速靜脈脈注射射的給給藥方方式下下，研研究血血藥濃濃度（（單位位體積積血液液中的的藥物物含量量）的的變化化規(guī)律律。tc2cc10問(wèn)題題2.給定定藥藥物物的的最最小小有有效效濃濃度度和和最最大大治治療療濃濃度度，，設(shè)設(shè)計(jì)計(jì)給給藥藥方方案案：：每每次次注注射射劑劑量量多多大大；；間間隔隔時(shí)時(shí)間間多多長(zhǎng)長(zhǎng)。。分析析理論論：：用用一一室室模模型型研研究究血血藥藥濃濃度度變變化化規(guī)規(guī)律律實(shí)驗(yàn)驗(yàn)：：對(duì)對(duì)血血藥藥濃濃度度數(shù)數(shù)據(jù)據(jù)作作擬擬合合，，符符合合負(fù)負(fù)指指數(shù)數(shù)變變化化規(guī)規(guī)律律3.血液液容容積積v,t=0注射射劑劑量量d,血藥藥濃濃度度立立即即為為d/v.2.藥物物排排除除速速率率與與血血藥藥濃濃度度成成正正比比，，比比例例系系數(shù)數(shù)k(>0)模型型假假設(shè)設(shè)1.機(jī)體體看看作作一一個(gè)個(gè)房房室室，，室室內(nèi)內(nèi)血血藥藥濃濃度度均均勻勻———一室室模模型型模型型建建立立在此此，，d=300mg，t及c（t）在在某某些些點(diǎn)點(diǎn)處處的的值值見(jiàn)見(jiàn)前前表表，，需需經(jīng)經(jīng)擬擬合合求求出出參參數(shù)數(shù)k、v用線(xiàn)線(xiàn)性性最最小小二二乘乘擬擬合合c(t)計(jì)算結(jié)果：d=300;t=[0.250.511.523468];c=[19.2118.1515.3614.1012.899.327.455.243.01];y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2))程序：給藥藥方方案案設(shè)設(shè)計(jì)計(jì)cc2c10t設(shè)每每次次注注射射劑劑量量D,間隔隔時(shí)時(shí)間間血藥藥濃濃度度c(t)應(yīng)c1c(t)c2初次次劑劑量量D0應(yīng)加加大大給藥方案記為：2、1、計(jì)算結(jié)果：給藥方案：c1=10,c2=25k=0.2347v=15.02故可可制制定定給給藥藥方方案案：：即:首次次注注射射375mg，其余余每每次次注注射射225mg，注射射的的間間隔隔時(shí)時(shí)間間為為4小時(shí)時(shí)。。估計(jì)計(jì)水水塔塔的的流流量量2、解解題題思思路路3、算算法法設(shè)設(shè)計(jì)計(jì)與與編編程程1、問(wèn)問(wèn)題題某居居民民區(qū)區(qū)有有一一供供居居民民用用水水的的園園柱柱形形水水塔塔，，一一般般可可以以通通過(guò)過(guò)測(cè)測(cè)量量其其水水位位來(lái)來(lái)估估計(jì)計(jì)水水的的流流量量，，但但面面臨臨的的困困難難是是，，當(dāng)當(dāng)水水塔塔水水位位下下降降到到設(shè)設(shè)定定的的最最低低水水位位時(shí)時(shí)，，水水泵泵自自動(dòng)動(dòng)啟啟動(dòng)動(dòng)向向水水塔塔供供水水，，到到設(shè)設(shè)定定的的最最高高水水位位時(shí)時(shí)停停止止供供水水，，這這段段時(shí)時(shí)間間無(wú)無(wú)法法測(cè)測(cè)量量水水塔塔的的水水位位和和水水泵泵的的供供水水量量．．通通常常水水泵泵每每天天供供水水一一兩兩次次，，每每次次約約兩兩小小時(shí)時(shí).水塔塔是是一一個(gè)個(gè)高高12.2米，，直直徑徑17.4米的的正正園園柱柱．．按按照照設(shè)設(shè)計(jì)計(jì)，，水水塔塔水水位位降降至至約約8.2米時(shí)時(shí)，，水水泵泵自自動(dòng)動(dòng)啟啟動(dòng)動(dòng)，，水水位位升升到到約約10.8米時(shí)時(shí)水水泵泵停停止止工工作作．．表1是某某一一天天的的水水位位測(cè)測(cè)量量記記錄錄，，試試估估計(jì)計(jì)任任何何時(shí)時(shí)刻刻（（包包括括水水泵泵正正供供水水時(shí)時(shí)））從從水水塔塔流流出出的的水水流流量量，，及及一一天天的的總總用用水水量量．．流量量估估計(jì)計(jì)的的解解題題思思路路擬合合水水位位~時(shí)間間函函數(shù)數(shù)確定定流流量量~時(shí)間間函函數(shù)數(shù)估計(jì)計(jì)一一天天總總用用水水量量擬合合水水位位~時(shí)間間函函數(shù)數(shù)測(cè)量量記記錄錄看看，，一一天天有有兩兩個(gè)個(gè)供供水水時(shí)時(shí)段段（（以以下下稱(chēng)稱(chēng)第第1供水水時(shí)時(shí)段段和和第第2供水水時(shí)時(shí)段段）），，和和3個(gè)水水泵泵不不工工作作時(shí)時(shí)段段（（以以下下稱(chēng)稱(chēng)第第1時(shí)段段t=0到t=8.97，第第2次時(shí)時(shí)段段t=10.95到t=20.84和第第3時(shí)段段t=23以后后））．．對(duì)對(duì)第第1、2時(shí)段段的的測(cè)測(cè)量量數(shù)數(shù)據(jù)據(jù)直直接接分分別別作作多多項(xiàng)項(xiàng)式式擬擬合合，，得得到到水水位位函函數(shù)數(shù)．．為為使使擬擬合合曲曲線(xiàn)線(xiàn)比比較較光光滑滑，，多多項(xiàng)項(xiàng)式式次次數(shù)數(shù)不不要要太太高高，，一一般般在在3~6．由由于于第第3時(shí)段段只只有有3個(gè)測(cè)測(cè)量量記記錄錄，，無(wú)無(wú)法法對(duì)對(duì)這這一一時(shí)時(shí)段段的的水水位位作作出出較較好好的的擬擬合合．．2、確定定流流量量~時(shí)間間函函數(shù)數(shù)對(duì)于于第第1、2時(shí)段段只只需需將將水水位位函函數(shù)數(shù)求求導(dǎo)導(dǎo)數(shù)數(shù)即即可可，，對(duì)對(duì)于于兩兩個(gè)個(gè)供供水水時(shí)時(shí)段段的的流流量量，，則則用用供供水水時(shí)時(shí)段段前前后后（（水水泵泵不不工工作作時(shí)時(shí)段段））的的流流量量擬擬合合得得到到，，并并且且將將擬擬合合得得到到的的第第2供水時(shí)段段流量外外推，將將第3時(shí)段流量量包含在在第2供水時(shí)段段內(nèi)．3、一天總用用水量的的估計(jì)總用水量量等于兩兩個(gè)水泵泵不工作作時(shí)段和和兩個(gè)供供水時(shí)段段用水量量之和，，它們都都可以由由流量對(duì)對(duì)時(shí)間的的積分得得到。算法設(shè)計(jì)計(jì)與編程程1、擬合第1、2時(shí)段的水水位，并并導(dǎo)出流流量2、擬合供水水時(shí)段的的流量3、估計(jì)一天天總用水水量4、流量及及總用水水量的檢檢驗(yàn)1、擬合第1時(shí)段的水水位，并并導(dǎo)出流流量設(shè)t，h為已輸入入的時(shí)刻刻和水位位測(cè)量記記錄（水水泵啟動(dòng)動(dòng)的4個(gè)時(shí)刻不不輸入）），第1時(shí)段各時(shí)刻的的流量可可如下得得：1）c1=polyfit（t（1：10），h（1：10），3）；%用3次多項(xiàng)式式擬合第第1時(shí)段水位位，c1輸出3次多項(xiàng)式式的系數(shù)數(shù)2）a1=polyder（c1）；%a1輸出多項(xiàng)項(xiàng)式（系系數(shù)為c1）導(dǎo)數(shù)的的系數(shù)3）tp1=0：0.1：9；x1=-polyval（a1，tp1）；%x1輸出多項(xiàng)項(xiàng)式（系系數(shù)為a1）在tp1點(diǎn)的函數(shù)數(shù)值（取取負(fù)后邊邊為正值值），即即tp1時(shí)刻的流流量4）流量函數(shù)數(shù)為：2、擬合第2時(shí)段的水水位，并并導(dǎo)出流流量設(shè)t，h為已輸入入的時(shí)刻刻和水位位測(cè)量記記錄（水水泵啟動(dòng)動(dòng)的4個(gè)時(shí)刻不不輸入）），第2時(shí)段各時(shí)刻的的流量可可如下得得：1）c2=polyfit(t(10.9:21),h(10.9:21),3)；%用3次多項(xiàng)式式擬合第第2時(shí)段水位位，c2輸出3次多項(xiàng)式式的系數(shù)數(shù)2）a2=polyder(c2)；%a2輸出多項(xiàng)項(xiàng)式（系系數(shù)為c2）導(dǎo)數(shù)的的系數(shù)3）tp2=10.9:0.1:21;x2=-polyval(a2,tp2);%x2輸出多項(xiàng)項(xiàng)式（系系數(shù)為a2）在tp2點(diǎn)的函數(shù)數(shù)值（取取負(fù)后邊邊為正值值），即即tp2時(shí)刻的流流量4）流量函數(shù)數(shù)為：3、擬合供水水時(shí)段的的流量在第1供水時(shí)段段（t=9~11）之前（（即第1時(shí)段）和和之后（（即第2時(shí)段）各各取幾點(diǎn)點(diǎn)，其流流量已經(jīng)經(jīng)得到，，用它們們擬合第第1供水時(shí)段段的流量量．為使使流量函函數(shù)在t=9和t=11連續(xù)，我我們簡(jiǎn)單單地只取取4個(gè)點(diǎn)，擬擬合3次多項(xiàng)式式（即曲曲線(xiàn)必過(guò)過(guò)這4個(gè)點(diǎn)），，實(shí)現(xiàn)如如下：xx1=-polyval（a1，[89]）；%取第1時(shí)段在t=8,9的流量xx2=-polyval（a2，[1112]）；%取第2時(shí)段在t=11,12的流量xx12=[xx1xx2]；c12=polyfit（[891112]，xx12，3）；%擬合3次多項(xiàng)式式tp12=9：0.1：11；x12=polyval（c12，tp12）；%x12輸出第1供水時(shí)段段各時(shí)刻的的流量擬合的流流量函數(shù)數(shù)為：在第2供水時(shí)段段之前取取t=20，20.8兩點(diǎn)的流流水量，，在該時(shí)時(shí)刻之后后（第3時(shí)段）僅僅有3個(gè)水位記記錄，我我們用差差分得到