插值和擬合講義

插值和擬合講義第1頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月[1]一維插值[2]二維插值第2頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第3頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第4頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第5頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月注：Hermite插值（略）第6頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月Runge現(xiàn)象第7頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第8頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第9頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月用MATLAB作插值計(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的范圍．第10頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月例1.用三次樣條插值選取11個(gè)基點(diǎn)計(jì)算插值x0=linspace(-5,5,11);y0=1./(1+x0.^2);x=linspace(-5,5,100);y=interp1(x0,y0,x,'spline');x1=linspace(-5,5,100);y1=1./(1+x1.^2);plot(x1,y1,'k',x0,y0,'+',x,y,'r')此例，可以看出插值函數(shù)得到的函數(shù)圖像與原函數(shù)很接近。第11頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

例2：從1點(diǎn)12點(diǎn)的11小時(shí)內(nèi)，每隔1小時(shí)測(cè)量一次溫度，測(cè)得的溫度的數(shù)值依次為：5，8，9，15，25，29，31，30，22，25，27，24．試估計(jì)每隔1/10小時(shí)的溫度值．hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');%(直接輸出數(shù)據(jù)將是很多的)plot(hours,temps,'+',h,t,'b',hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius')第12頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

xy機(jī)翼下輪廓線(xiàn)例3.已知飛機(jī)下輪廓線(xiàn)上數(shù)據(jù)如下，求x每改變0.1時(shí)的y值．第13頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月x0=[035791112131415];y0=[01.21.72.02.12.01.81.21.01.6];x=0:0.1:15;y1=interp1(x0,y0,x);y2=interp1(x0,y0,x,'spline');subplot(2,1,1)plot(x0,y0,'k+',x,y1,'r')gridtitle('piecewiselinear')subplot(2,1,2)plot(x0,y0,'k+',x,y2,'r')gridtitle('spline')第14頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第15頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第16頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第17頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第18頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

要求x0,y0單調(diào)；x，y可取為矩陣，或x取行向量，y取為列向量，x,y的值分別不能超出x0,y0的范圍．z=interp2(x0,y0,z0,x,y,’method’)被插值點(diǎn)插值方法用MATLAB作網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)的插值插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值‘nearest’

最鄰近插值；‘linear’

雙線(xiàn)性插值（缺省值）；‘cubic’

雙三次插值；缺省時(shí)雙線(xiàn)性插值.第19頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月例4.測(cè)得平板表面3×5網(wǎng)格點(diǎn)處的溫度分別為：828180828479636165818484828586試作出平板表面的溫度分布曲面z=f(x,y)的圖形．輸入以下命令：x=1:5;y=1:3;temps=[8281808284;7963616581;8484828586];mesh(x,y,temps)1.先在三維坐標(biāo)畫(huà)出原始數(shù)據(jù)，畫(huà)出粗糙的溫度分布曲線(xiàn)圖.2．以平滑數(shù)據(jù),在x、y方向上每隔0.2個(gè)單位的地方進(jìn)行插值.第20頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫(huà)出插值后的溫度分布曲面圖.第21頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

通過(guò)此例對(duì)最近鄰點(diǎn)插值、雙線(xiàn)性插值方法和雙三次插值方法的插值效果進(jìn)行比較．第22頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月figure(1);meshz(x,y,z)xlabel('X')ylabel('Y')zlabel('Z')xi=0:50:5600;yi=0:50:4800;figure(2)z1i=interp2(x,y,z,xi,yi','nearest');surfc(xi,yi,z1i)xlabel('X'),ylabel('Y'),zlabel('Z')figure(3)z2i=interp2(x,y,z,xi,yi');surfc(xi,yi,z2i)xlabel('X'),ylabel('Y'),zlabel('Z')

figure(4)z3i=interp2(x,y,z,xi,yi','cubic');surfc(xi,yi,z3i)xlabel('X'),ylabel('Y'),zlabel('Z')figure(5)subplot(1,3,1),contour(xi,yi,z1i,10,'r');subplot(1,3,2),contour(xi,yi,z2i,10,'r');subplot(1,3,3),contour(xi,yi,z3i,10,'r');第23頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

插值函數(shù)griddata格式為:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散點(diǎn)數(shù)據(jù)的插值計(jì)算

要求cx取行向量，cy取為列向量．被插值點(diǎn)插值方法插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值‘nearest’最鄰近插值‘linear’

雙線(xiàn)性插值‘cubic’

雙三次插值'v4'-MATLAB提供的插值方法缺省時(shí),雙線(xiàn)性插值第24頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

例6.在某海域測(cè)得一些點(diǎn)(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）×（-50，150）里的哪些地方船要避免進(jìn)入．第25頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月x=[129140103.588185.5195105.5157.5107.57781162162117.5];y=[7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5];z=[-4-8-6-8-6-8-8-9-9-8-8-9-4-9];cx=75:0.5:200;cy=-50:0.5:150;cz=griddata(x,y,z,cx,cy','cubic');第26頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月4.作出水深小于5的海域范圍,即z=5的等高線(xiàn)...3作海底曲面圖meshz(cx,cy,cz),rotate3dxlabel('X'),ylabel('Y'),zlabel('Z')figure(2),contour(cx,cy,cz,[-5-5]);gridholdonplot(x,y,'+')xlabel('X'),ylabel('Y')第27頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月用MATLAB解擬合問(wèn)題1.線(xiàn)性最小二乘擬合2.非線(xiàn)性最小二乘擬合第28頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第29頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第30頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第31頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第32頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第33頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第34頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月用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）輸出擬合多項(xiàng)式系數(shù)a=[a1,…,am,

am+1](數(shù)組)）輸入同長(zhǎng)度的數(shù)組x，y擬合多項(xiàng)式次數(shù)第35頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月即要求出二次多項(xiàng)式:中的使得:例7.對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合第36頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月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é)果：Ａ

=-9.810820.1293-0.0317第37頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月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é)果：Ａ=-9.810820.1293-0.0317解法2．用多項(xiàng)式擬合的命令第38頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月1.lsqcurvefit已知數(shù)據(jù)點(diǎn)：xdata=（xdata1,xdata2,…,xdatan）,ydata=（ydata1,ydata2,…,ydatan）

用MATLAB作非線(xiàn)性最小二乘擬合

MATLAB提供了兩個(gè)求非線(xiàn)性最小二乘擬合的函數(shù)：lsqcurvefit和lsqnonlin．兩個(gè)命令都要先建立M文件fun.m，在其中定義函數(shù)f(x)，但兩者定義f(x)的方式是不同的,可參考例題.

用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x,xdata1）,…,F（x，xdatan））T中的參變量x(向量),使得

第39頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

輸入格式為:

(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)化第40頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月lsqnonlin用以求含參量x（向量）的向量值函數(shù)

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

中的參量x，使得

最?。渲衒i（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知數(shù)據(jù)點(diǎn)：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）第41頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月輸入格式為：

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)化第42頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

例8.用下面一組數(shù)據(jù)擬合中的參數(shù)a，b，k該問(wèn)題即解最優(yōu)化問(wèn)題：第43頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

1）編寫(xiě)M文件curvefun1.m

functionf=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.用命令lsqcurvefit第44頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月3）運(yùn)算結(jié)果為：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結(jié)論：a=0.0063,b=-0.0034,k=0.2542第45頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

解法2

用命令lsqnonlin

f(x)=F(x,tdata,ctada)=

x=（a，b，k）1）編寫(xiě)M文件curvefun2.m

functionf=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=cdata-x(1)-x(2)*exp(-0.02*x(3)*tdata)

2）輸入命令:

x0=[0.2,0.05,0.05];x=lsqnonlin('curvefun2',x0)f=curvefun2(x)函數(shù)curvefun2的自變量是x，cdata和tdata是已知參數(shù)，故應(yīng)將cdatatdata的值寫(xiě)在curvefun2.m中第46頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月3）運(yùn)算結(jié)果為f=1.0e-003*(0.2322-0.1243-0.2495-0.2413-0.1668-0.07240.02410.11590.20300.2792x=0.0063-0.00340.2542可以看出,兩個(gè)命令的計(jì)算結(jié)果是相同的.4）結(jié)論：即擬合得a=0.0063b=-0.0034k=0.2542第47頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月MATLAB解應(yīng)用問(wèn)題實(shí)例1.電阻問(wèn)題2.給藥方案問(wèn)題*3.水塔流量估計(jì)問(wèn)題第48頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月電阻問(wèn)題溫度t(oC)20.532.751.073.095.7電阻R()7658268739421032例.由數(shù)據(jù)擬合R=a1t+b方法1.用命令polyfit(x,y,m)t=[20.532.5517395.7];r=[7658268739421032];aa=polyfit(t,r,1);a=aa(1)b=aa(2)y=polyval(aa,t);plot(t,r,'k+',t,y,'r')得到a=3.3940,b=702.4918第49頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月方法2.直接用結(jié)果相同．t=[20.532.5517395.7];r=[7658268739421032];R=[t'ones(5,1)];aa=R\r';a=aa(1)b=aa(2)y=polyval(aa,t);plot(t,r,'k+',t,y,'r')第50頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月一室模型：將整個(gè)機(jī)體看作一個(gè)房室，稱(chēng)中心室，室內(nèi)血藥濃度是均勻的．快速靜脈注射后，濃度立即上升；然后迅速下降．當(dāng)濃度太低時(shí)，達(dá)不到預(yù)期的治療效果；當(dāng)濃度太高，又可能導(dǎo)致藥物中毒或副作用太強(qiáng)．臨床上，每種藥物有一個(gè)最小有效濃度c1和一個(gè)最大有效濃度c2．設(shè)計(jì)給藥方案時(shí)，要使血藥濃度保持在c1~c2之間．本題設(shè)c1=10ug/ml，c2=25ug/ml.擬合問(wèn)題實(shí)例2給藥方案——一種新藥用于臨床之前，必須設(shè)計(jì)給藥方案.

藥物進(jìn)入機(jī)體后通過(guò)血液輸送到全身，在這個(gè)過(guò)程中不斷地被吸收、分布、代謝，最終排出體外，藥物在血液中的濃度，即單位體積血液中的藥物含量，稱(chēng)為血藥濃度．第51頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

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

要設(shè)計(jì)給藥方案,必須知道給藥后血藥濃度隨時(shí)間變化的規(guī)律．從實(shí)驗(yàn)和理論兩方面著手：第52頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月給藥方案1.在快速靜脈注射的給藥方式下，研究血藥濃度（單位體積血液中的藥物含量）的變化規(guī)律．tc2cc1O

問(wèn)題2.給定藥物的最小有效濃度和最大治療濃度，設(shè)計(jì)給藥方案：每次注射劑量多大；間隔時(shí)間多長(zhǎng)．分析

理論：用一室模型研究血藥濃度變化規(guī)律

實(shí)驗(yàn)：對(duì)血藥濃度數(shù)據(jù)作擬合，符合負(fù)指數(shù)變化規(guī)律第53頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月3.血液容積v,t=0注射劑量d,血藥濃度立即為d/v.2.藥物排除速率與血藥濃度成正比，比例系數(shù)k(>0)模型假設(shè)1.機(jī)體看作一個(gè)房室，室內(nèi)血藥濃度均勻——一室模型模型建立

在此，d=300mg，t及c（t）在某些點(diǎn)處的值見(jiàn)前表，需經(jīng)擬合求出參數(shù)k、v.第54頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月用線(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))程序：第55頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月用非線(xiàn)性最小二乘擬合c(t)-用lsqcurvefit2.主程序如下cleartdata=[0.250.511.523468];cdata=[19.2118.1515.3614.1012.899.327.455.243.01];x0=[10,0.5];x=lsqcurvefit('curvefun3',x0,tdata,cdata);f=curvefun3(x,tdata)x1.用M文件curvefun3.m定義函數(shù)functionf=curvefun3(x,tdata)d=300f=(x(1)\d)*exp(-x(2)*tdata)%x(1)=v;x(2)=k第56頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月給藥方案設(shè)計(jì)cc2c1O

t

設(shè)每次注射劑量D,間隔時(shí)間

血藥濃度c(t)

應(yīng)c1

c(t)

c2

初次劑量D0應(yīng)加大給藥方案記為：2.1.計(jì)算結(jié)果：給藥方案：c1=10,c2=25k=0.2347v=15.02第57頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月故可制定給藥方案：即:

首次注射375mg，其余每次注射225mg，注射的間隔時(shí)間為4h．第58頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月估計(jì)水塔的流量2.解題思路3.算法設(shè)計(jì)與編程1.問(wèn)題第59頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

某居民區(qū)有一供居民用水的圓柱形水塔，一般可以通過(guò)測(cè)量其水位來(lái)估計(jì)水的流量，但面臨的困難是，當(dāng)水塔水位下降到設(shè)定的最低水位時(shí)，水泵自動(dòng)啟動(dòng)向水塔供水，到設(shè)定的最高水位時(shí)停止供水，這段時(shí)間無(wú)法測(cè)量水塔的水位和水泵的供水量．通常水泵每天供水一兩次，每次約兩小時(shí).水塔是一個(gè)高12.2m，直徑17.4m的正圓柱．按照設(shè)計(jì)，水塔水位降至約8.2m時(shí)，水泵自動(dòng)啟動(dòng)，水位升到約10.8m時(shí)水泵停止工作．表1是某一天的水位測(cè)量記錄，試估計(jì)任何時(shí)刻（包括水泵正供水時(shí)）從水塔流出的水流量，及一天的總用水量．第60頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月第61頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月流量估計(jì)的解題思路擬合水位~時(shí)間函數(shù)確定流量~時(shí)間函數(shù)估計(jì)一天總用水量第62頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

擬合水位~時(shí)間函數(shù)從測(cè)量記錄看，一天有兩個(gè)供水時(shí)段（以下稱(chēng)第1供水時(shí)段和第2供水時(shí)段），和3個(gè)水泵不工作時(shí)段（以下稱(chēng)第1時(shí)段t=0到t=8.97，第2次時(shí)段t=10.95到t=20.84和第3時(shí)段t=23以后）．對(duì)第1、2時(shí)段的測(cè)量數(shù)據(jù)直接分別作多項(xiàng)式擬合，得到水位函數(shù)．為使擬合曲線(xiàn)比較光滑，多項(xiàng)式次數(shù)不要太高，一般在3~6．由于第3時(shí)段只有3個(gè)測(cè)量記錄，無(wú)法對(duì)這一時(shí)段的水位作出較好的擬合．第63頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

確定流量~時(shí)間函數(shù)

對(duì)于第1、2時(shí)段只需將水位函數(shù)求導(dǎo)數(shù)即可，對(duì)于兩個(gè)供水時(shí)段的流量，則用供水時(shí)段前后（水泵不工作時(shí)段）的流量擬合得到，并且將擬合得到的第2供水時(shí)段流量外推，將第3時(shí)段流量包含在第2供水時(shí)段內(nèi)．第64頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

一天總用水量的估計(jì)

總用水量等于兩個(gè)水泵不工作時(shí)段和兩個(gè)供水時(shí)段用水量之和，它們都可以由流量對(duì)時(shí)間的積分得到．第65頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月算法設(shè)計(jì)與編程1.擬合第1、2時(shí)段的水位，并導(dǎo)出流量2.擬合供水時(shí)段的流量3.估計(jì)一天總用水量4.流量及總用水量的檢驗(yàn)第66頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月1.擬合第1時(shí)段的水位，并導(dǎo)出流量設(shè)t，h為已輸入的時(shí)刻和水位測(cè)量記錄（水泵啟動(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ù)2）a1=polyder（c1）；

%a1輸出多項(xiàng)式（系數(shù)為c1）導(dǎo)數(shù)的系數(shù)

3）tp1=0：0.1：9；

x1=-polyval（a1，tp1）；%x1輸出多項(xiàng)式（系數(shù)a1）在tp1點(diǎn)的函數(shù)值（取負(fù)后邊為正值），即tp1時(shí)刻的流量

4）流量函數(shù)為：第67頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

擬合第2時(shí)段的水位，并導(dǎo)出流量設(shè)t，h為已輸入的時(shí)刻和水位測(cè)量記錄（水泵啟動(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ù)2）a2=polyder(c2)；

%a2輸出多項(xiàng)式（系數(shù)為c2）導(dǎo)數(shù)的系數(shù)

3）tp2=10.9:0.1:21;x2=-polyval(a2,tp2);%x2輸出多項(xiàng)式（系數(shù)為a2）在tp2點(diǎn)的函數(shù)值（取負(fù)后邊為正值），即tp2時(shí)刻的流量4）流量函數(shù)為：第68頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

2.擬合供水時(shí)段的流量在第1供水時(shí)段（t=9~11）之前（即第1時(shí)段）和之后（即第2時(shí)段）各取幾點(diǎn)，其流量已經(jīng)得到，用它們擬合第1供水時(shí)段的流量．為使流量函數(shù)在t=9和t=11連續(xù)，我們簡(jiǎn)單地只取4個(gè)點(diǎn)，擬合3次多項(xiàng)式（即曲線(xiàn)必過(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ù)為：第69頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月

在第2供水時(shí)段之前取t=20，20.8兩點(diǎn)的流水量，在該時(shí)刻之后（第3時(shí)段）僅有3個(gè)水位記錄，我們用差分得到流量，然后用這4個(gè)數(shù)值擬合第2供水時(shí)段的流量如下：

dt3=diff（t(22：24)）；

%最后3個(gè)時(shí)刻的兩兩之差

dh3=diff（h(22：24)）；

%最后3個(gè)水位的兩兩之差dht3=-dh3./dt3；

%t(22)和t(23)的流量t3=[2020.8t(22)t(23)]；

xx3=[-polyval(a2，t3(1：2)，dht3)]；

%取t3各時(shí)刻的流量

c3=polyfit（t3，xx3，3）；%擬合3次多項(xiàng)式

t3=20.8：0.1：24；

x3=polyval（c3，tp3）；%x3輸出第2供水時(shí)段（外推至t=24）各時(shí)刻的流量擬合的流量函數(shù)為：第70頁(yè)，課件共77頁(yè)，創(chuàng)作于2023年2月3.一天總用水量的估計(jì)第1、2時(shí)段和第1、2供水時(shí)段流量的積分之和，就是一天總用水量