插值与拟合-MATLAB学习之四

1、 MATLAB与插值、拟合与插值、拟合主要参考资料：数学建模与数学实验数学建模与数学实验第三版第三版赵静，但琦等编，高等教育出版社，2010年 一、一、 插插 值值1、插值的定义、插值的定义2、插值的方法、插值的方法3、用、用Matlab解插值问题解插值问题问题来源: 当数据量不够，需要补充，且认定已有数据可信时 , 通常利用函数插值方法建立插值模型. 当数据量不够，需要补充，且认定已有数据可信时 , 通常利用函数插值方法建立插值模型.目标：根据一组观测数据寻找函数关系( ,)0,1, .iix yin( )yf x满足( )0,1, .iiyf xin 一、一、 插插 值值0 x1xnx0y

2、1y1 1、插值的定义、插值的定义, 1 , 0(),(njyxjj其中其中jx互不相同，不妨设互不相同，不妨设),10bxxxan求任一插值点求任一插值点)(*jxx 处的插值处的插值.*y0 x1xnx0y1y*x*y已知已知 n+1个数据点个数据点(节点节点): 构造一个构造一个(相对简单的相对简单的)函数函数),(xfy 通过全部节点通过全部节点, 即即), 1 ,0()(njyxfjj再用再用)(xf计算插值，即计算插值，即).(*xfy 0 x1xnx0y1y*x*y 称为拉格朗日插值基函数拉格朗日插值基函数。n0iiiny)x(L)x(P 已知函数f(x)在n+1个点x0,x1,

3、xn处的函数值为 y0,y1,yn 。求一n次多项式函数Pn(x)，使其满足： Pn(xi)=yi,i=0,1,n. 解决此问题的拉格朗日插值多项式公式如下其中Li(x) 为n次多项式：)xx()xx)(xx()xx)(xx()xx()xx)(xx()xx)(xx()x(Lni1ii1ii1i0in1i1i10i(1). (1). 拉格朗日拉格朗日(Lagrange)插值插值2、插值的方法、插值的方法(1).(1).拉格朗日拉格朗日(Lagrange)插值插值特别地特别地:两点一次两点一次(线性线性)插值多项式插值多项式: 101001011yxxxxyxxxxxL三点二次三点二次(抛物抛物)

4、插值多项式插值多项式: 2120210121012002010212yxxxxxxxxyxxxxxxxxyxxxxxxxxxL .,满足插值条件直接验证可知xLn 拉格朗日多项式插值的这种振荡现象叫 Runge现象现象55,11)(2xxxg 采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例例-5-4-3-2-1012345-1.5-1-0.500.511.52y=1/(1+x2)n=2n=4n=6n=8n=10(2).(2).分段线性插值分段线性插值其它,0,)()()(1111110jjjjjjjjjjjn

5、jjjnxxxxxxxxxxxxxxxlxlyxL计算量与n无关;n越大，误差越小.nnnxxxxgxL0),()(limxjxj-1xj+1x0 xnxoy66,11)(2xxxg例例用分段线性插值法求插值用分段线性插值法求插值,并观察插值误差并观察插值误差.1.在在-6,6中平均选取中平均选取5个点作插值个点作插值(xch11)4.在在-6,6中平均选取中平均选取41个点作插值个点作插值(xch14)2.在在-6,6中平均选取中平均选取11个点作插值个点作插值(xch12)3.在在-6,6中平均选取中平均选取21个点作插值个点作插值(xch13)比分段线性插值更光滑。比分段线性插值更光滑。

6、xyxi-1 xiab 在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。 光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值三次样条插值 (3).(3).三次样条插值, 1,),()(1nixxxxsxSiii,)()3), 1 ,0()()2), 1()()10223niiiiiiixxCxSniyxSnidxcxbxaxs) 1, 1()()(),()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii自然边界条件）(0)()()40 nxSxS)(,)4

7、)3)2xSdcbaiiii)()(limxgxSng g( (x x) )为被插值函数为被插值函数。例例66,11)(2xxxg用三次样条插值选取用三次样条插值选取11个基点计算插值个基点计算插值(ych)-5-4-3-2-101234500.10.20.30.40.50.60.70.80.91x0=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);用用M

8、ATLABMATLAB作插值计算小结作插值计算小结一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ： 线性插值；线性插值；spline ： 三次样条插三次样条插值；值；cubic ： 立方插值。立方插值。缺省时：缺省时： 分段线性插值。分段线性插值。 注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范围。 例：在例：在1-121-12的的1111小时内，每隔小

9、时内，每隔1 1小时测量一次小时测量一次温度，测得的温度依次为：温度，测得的温度依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424。试估计每隔。试估计每隔1/101/10小时的小时的温度值。温度值。hours=1:12;temps=5 8 9 15 25 29 31 30 22 25 27 24;h=1:0.1:12;t=interp1(hours,temps,h,spline); (直接输出数据将是很多的)plot(hours,temps,+,h,t,hours,temps,r:) %作图xlabel(Hour),ylab

10、el(Degrees Celsius)X035791 11 21 31 41 5Y01 . 21 . 72 . 02 . 12 . 01 . 81 . 21 . 01 . 6例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值。值。x0=0 3 5 7 9 11 12 13 14 15 ;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ;x=0:0.1:15;y1=lagr1(x0,y0,x);y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,spline); subplot(3,1,1

11、)plot(x0,y0,k+,x,y1,r)gridtitle(lagrange)subplot(3,1,2)plot(x0,y0,k+,x,y2,r)Grid;title(piecewise linear)subplot(3,1,3)plot(x0,y0,k+,x,y3,r)Grid;title(spline)机翼下轮廓线xy051015-20-10010lagrange0510150123piecewise linear0510150123spline二、拟合二、拟合2.2.拟合的基本原理拟合的基本原理1. 拟合问题引例拟合问题引例拟拟 合合 问问 题题 引引 例例 1 1温度温度t(0C

12、) 20.5 32.7 51.0 73.0 95.7电阻电阻R( ) 765 826 873 942 1032已知热敏电阻数据：已知热敏电阻数据：求求60600C时的电阻时的电阻R。2040608010070080090010001100 设设 R=at+ba,b为待定系数为待定系数拟拟 合合 问问 题题 引引 例例 2 2 t (h) 0.25 0.5 1 1.5 2 3 4 6 8c ( g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注

13、射300mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).作半对数坐标系作半对数坐标系(semilogy)下的图形下的图形为待定系数kcectckt,)(002468100101102曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上 n个点个点（xi,yi) i=1,n, 寻求一个函数（曲线）寻求一个函数（曲线）y=f(x), 使使 f(x) 在某种准则下与所在某种准则下与所有数据点最为接近，即曲线拟合得最好。有数据点最为接近，即曲线拟合得最好。 +xyy=f(x)(xi,yi)i i 为点为点（xi,yi) 与

14、与曲线曲线 y=f(x) 的距离的距离拟合与插值的关系拟合与插值的关系 函数插值与曲线拟合都是要根据一组数据构造一个函数作函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同为近似，由于近似的要求不同，二者的数学方法上是完全不同的。的。 实例：实例：下面数据是某次实验所得，希望得到X和 f之间的关系？x124791 21 31 51 7f1 .53 .96 .611 .71 5 .61 8 .81 9 .62 0 .62 1 .1问题：问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：解决方案：若不要求曲线（面）通过所有数据点，而

15、是要求它反映对象整体的变化趋势，这就是数据拟合数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所有数据点，就是插值问题插值问题；x=1 2 4 7 9 12 13 15 17;f=1.5 3.9 6 11.7 12.6 18.8 20.3 20.6 21.1;axis(0 18 1 22)xlabel(x); ylabel(f)y=1:0.1:17; figure(1);plot(x,f,o); gtext(已知数据点)hold on;pausebb1=interp1(x,f,y,nearest)plot(y,bb1); gtext(nearest)hold on;pausea=

16、polyfit(x,f,3)aa=polyval(a,y)plot(y,aa)gtext(曲线拟合); hold off;pausefigure(2);plot(x,f,o)gtext(已知数据点)hold on;pausebb1=interp1(x,f,y,linest)plot(y,bb1)gtext(linest)hold on;pausea=polyfit(x,f,3)aa=polyval(a,y)plot(y,aa)gtext(曲线拟合)hold off;pause%最临近插值、线性插值、样条插值与曲线拟合结果：最临近插值、线性插值、样条插值与曲线拟合结果：0246810121416

17、180510152025已知数据点spline三次多项式插值0246810121416180510152025已知数据点linest三次多项式插值0246810121416180510152025已知数据点nearest三次多项式插值曲线拟合问题最常用的解法曲线拟合问题最常用的解法线性最小二乘法的基本思路线性最小二乘法的基本思路第一步: :先选定一组函数先选定一组函数 r1(x), r2(x), rm(x), mn, 令令 f(x)=a1r1(x)+a2r2(x)+ +amrm(x) （1）其中其中 a1,a2, am 为待定系数。为待定系数。 第二步: 确定确定a1,a2, am 的准则（最

18、小二乘准则）：的准则（最小二乘准则）：使使n个点个点（xi,yi) 与与曲线曲线 y=f(x) 的距离的距离 i 的平方和最小的平方和最小 。记记 )2()()(),(211211221iiknimkkininiiimyxrayxfaaaJ 问题归结为，求问题归结为，求 a1,a2, am 使使 J(a1,a2, am) 最小。最小。线性最小二乘法的求解：预备知识线性最小二乘法的求解：预备知识超定方程组超定方程组：方程个数大于未知量个数的方程组：方程个数大于未知量个数的方程组)( 221111212111mnyarararyarararnmnmnnmm即即 Ra=ynmnmnnmyyyaaar

19、rrrrrR112111211,其中其中超定方程一般是不存在解的矛盾方程组。超定方程一般是不存在解的矛盾方程组。 如果有向量如果有向量a使得使得 达到最小，达到最小，则称则称a为上述为上述超定方程的最小二乘解超定方程的最小二乘解。 212211)(imniimiiyararar线性最小二乘法的求解线性最小二乘法的求解 定理：定理：当当R RT TR R可逆时，超定方程组（可逆时，超定方程组（3 3）存在最小二乘解，）存在最小二乘解，且即为方程组且即为方程组 R RT TRa=RRa=RT Ty y的解：的解：a=(Ra=(RT TR)R)-1-1R RT Ty y 所以，曲线拟合的最小二乘法要

20、解决的问题，实际上就是所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。求以下超定方程组的最小二乘解的问题。nmnmnmyyyaaaxrxrxrxrR111111,)()()()(其中其中Ra=y （3）线性最小二乘拟合线性最小二乘拟合 f(x)=a1r1(x)+ +amrm(x)中中函数函数rr1 1(x), (x), r rm m(x)(x)的选取的选取 1. 1. 通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 f(x)f(x)；+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-

21、bx 2. 2. 将数据将数据 (xi,yi) i=1, n 作图，通过直观判断确定作图，通过直观判断确定 f(x)：三三 用用MATLABMATLAB解拟合问题解拟合问题1 1、线性最小二乘拟合、线性最小二乘拟合2 2、非线性最小二乘拟合、非线性最小二乘拟合用用MATLAB作线性最小二乘拟合作线性最小二乘拟合1. 1. 作多项式作多项式f(x)=a1xm+ +amx+am+1拟合拟合, ,可利用已有程序可利用已有程序:a=polyfit(x,y,m)2. 2. 对超定方程组对超定方程组)(11nmyaRnmmn可得最小二乘意义下的解。可得最小二乘意义下的解。，用，用yRa3.3.多项式在多项

22、式在x x处的值处的值y y可用以下命令计算：可用以下命令计算： y=polyvaly=polyval（a a，x x）输出拟合多项式系数输出拟合多项式系数a=a1, am , am+1 (数组数组) )）输入同长度输入同长度的数组的数组X，Y拟合多项拟合多项式次数式次数即要求即要求 出二次多项式出二次多项式:3221)(axaxaxf中中 的的),(321aaaA 使得使得:最小 )(1112iiiyxf例例 对下面一组数据作二次多项式拟合对下面一组数据作二次多项式拟合xi0.10.20.40.50.60.70.80.91yi1.9783.286.167.347.669.589.489.30

23、11.21）输入以下命令）输入以下命令：x=0:0.1:1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2; R=(x.2) x ones(11,1)； A=Ry11 11211121xxxxR此时解法解法1 1用解超定方程的方法用解超定方程的方法2）计算结果）计算结果： = -9.8108 20.1293 -0.03170317.01293.208108.9)(2xxxf1）输入以下命令：）输入以下命令： x=0:0.1:1; y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9

24、.48 9.30 11.2; A=polyfit(x,y,2) z=polyval(A,x); plot(x,y,k+,x,z,r) %作出数据点和拟合曲线的图形作出数据点和拟合曲线的图形2）计算结果：）计算结果： = -9.8108 20.1293 -0.0317解法解法2用多项式拟合的命令用多项式拟合的命令00.20.40.60.81-20246810120317.01293.208108.9)(2xxxf1. lsqcurvefit1. lsqcurvefit已知数据点数据点： xdata=xdata=（xdata1，xdata2，xdataxdatan n），）， ydata=ydat

25、a=（ydataydata1 1，ydataydata2 2，ydataydatan n） 3.3.用用MATLAB作非线性最小二乘拟合作非线性最小二乘拟合 Matlab Matlab的提供了两个求非线性最小二乘拟合的函数：的提供了两个求非线性最小二乘拟合的函数：lsqcurvefitlsqcurvefit和lsqnonlinlsqnonlin。两个命令都要先建立。两个命令都要先建立M-M-文件文件fun.mfun.m，在其中定义函数在其中定义函数f(x)f(x)，但两者定义，但两者定义f(x)f(x)的方式是不同的的方式是不同的, ,可参可参考例题考例题.最小 ),(21niiiydatax

26、dataxF lsqcurvefitlsqcurvefit用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数F(x,xdata)=F(x,xdata)=（F F（x x，xdataxdata1 1），），F F（x x，xdataxdatan n）T T中的参变量中的参变量x(x(向量向量),),使得使得 输入格式为输入格式为: : (1) x = lsqcurvefit (fun,x0,xdata,ydata); (2) x =lsqcurvefit (fun,x0,xdata,ydata,options); (3) x = lsqcurvefit (fun,x0,xda

27、ta,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是一个事先建立的是一个事先建立的定义函数定义函数F(x,xdata) 的的M-文件文件, 自变量为自变量为x和和xdata说明：x = lsqcurvefit (fun,x0,xdata,yda

28、ta,options);迭代初值迭代初值已知数据点已知数据点选项见无选项见无约束优化约束优化 lsqnonlin用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数 f(x)f(x)=(f=(f1 1(x),f(x),f2 2(x),(x),f,fn n(x)(x)T T 中的参量中的参量x x，使得，使得 最小。最小。 其中其中 fi（x）=f（x，xdatai，ydatai） =F(x,xdatai)-ydatai 22221)()()()()(xfxfxfxfxfnT2. lsqnonlin已知数据点：已知数据点： xdata=xdata=（xdata1，xdata2，

29、xdataxdatan n） ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）输入格式为：输入格式为： 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，）；说明：x= lsqnonlinlsqnonlin （fun，x0，options）；）；

30、fun是一个事先建立的是一个事先建立的定义函数定义函数f(x)的的M-文件，文件，自变量为自变量为x迭代初值迭代初值选项见无选项见无约束优化约束优化100200 30040050060070080090010004.54 4.99 5.35 5.65 5.90 6.10 6.26 6.39 6.50 6.59jt310jc210102. 0),(minjjktcbeakbaFj 例例2 用下面一组数据拟合用下面一组数据拟合 中的参数中的参数a，b，kktbeatc2 . 0 . 0)(该问题即解最优化问题：该问题即解最优化问题： 1 1）编写）编写M-M-文件文件 curvefun1.mcur

31、vefun1.m function f=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:1000tdata=100:100:1000cdata=cdata=1e-03* *4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;6.50,6.59; x0=0.2,0.05,0.05; x0=0.2,0.05,0

32、.05; x=lsqcurvefit (curvefun1,x0,tdata,cdata)x=lsqcurvefit (curvefun1,x0,tdata,cdata) f= f= curvefun1(x,tdata) F(x，tdata)= ，x=(a，b，k)Tktktbeabea),(10102. 002. 0解法解法1 1. 用命令用命令lsqcurvefitlsqcurvefit3 3）运算结果为）运算结果为：f =0.0043 0.0051 0.0056 0.0059 0.0061 f =0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.00

33、62 0.0063 0.0063 0.0063 0.0062 0.0062 0.0063 0.0063 0.0063 x = 0.0063 -0.0034 0.2542 x = 0.0063 -0.0034 0.25424）结论）结论：a=0.0063, b=-0.0034, k=0.2542Tktktcbeacbea),(102. 0102. 0101 解法解法 2 用命令用命令lsqnonlin f(x)=F(x,tdata,ctada)= x=（a，b，k）1）编写编写M-M-文件文件 curvefun2.mcurvefun2.m function f=curvefun2(x) tdat

34、a=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)函数函数curvefun2的自变量是的自变量是x，cdata和和tdata是已知参数，故应是已知参数，故应将将cdata tdata的值写在的值写在curvefun2.m中中3 3）运算结果为）运算结果为 f =1.0e

35、-003 f =1.0e-003 * *(0.2322 -0.1243 -0.2495 -0.2413 (0.2322 -0.1243 -0.2495 -0.2413 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792 x =0.0063 -0.0034 0.2542 x =0.0063 -0.0034 0.2542可以看出可以看出,两个命令的计算结果是相同的两个命令的计算结果是相同的.4）结论）结论：即拟合得即拟合得a=0.0063 b=-0.0034 k=0.25420

36、.0063 b=-0.0034 k=0.2542估计水塔的流量估计水塔的流量(AMCM 91)四 MATLAB解应用问题实例 某居民区有一供居民用水的园柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量通常水泵每天供水一两次，每次约两小时.一、问题重述一、问题重述 水塔是一个高12.2米，直径17.4米的正园柱按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作表1 是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）

37、从水塔流出的水流量，及一天的总用水量 表 1 水位测量记录 （符号/表示水泵启动）时刻(h)水位(cm)0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97968 948 931 913 898 881 869 852 839 822时刻(h)水位(cm)9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93/ / 1082 1050 1021 994 965 941 918 892时刻(h)水位(cm)19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.

38、91866 843 822 / / 1059 1035 1018二、问题分析二、问题分析我们很容易想到应通过对所给的数据进行数值拟合来建模在讨论具体的建模方法以前，我们应先给出一些合理的假设(1) 影响水从水塔中流出的流量的唯一因素是公众对水的传统要求因为表只给出了某一天（近26小时）水塔的水位数据，并没有对这些数据的产生有影响的因素作出具体的说明，我们只能假定所给数据反映了有代表性的一天，而不包括任何特殊情况，如自然灾害、火灾、水塔溢水、水塔漏水等对水的特殊要求(2) 水塔中的水位不影响水流量的大小，气候条件、温度变化等也不影响水流量(3) 水泵工作起止时间有它的水位决定，每次充水时间大约为

39、两个小时(4) 水泵充水速度恒定，且水泵充水的水流量远大于水塔的水流量，以保证人们对水的需求水泵工作时不需要维修，也不中途停止工作(5) 水塔的水流量与水泵状态独立，并不因水泵工作而增加或减少水流量的大小(6) 水塔的水流量曲线可以用一条光滑的曲线了逼近这时，在每一个数据点，水流量的两阶导数是连续的因为水的消耗量是基于社区公众一天的活动，如洗澡、做饭、洗衣服等，每一个使用者的要求与整个社区的要求相比是微不足道的，而整个社区的需求是不可能同时增加或减少的，由于水的消耗的自然性，可以假设水流量曲线是一条连续光滑的曲线(7) 表的数据是准确的.对所给的问题，其建模方法是经典对所给的问题，其建模方法是

40、经典的，基本上是分成三步的，基本上是分成三步: 1. 拟合水位拟合水位时间函数时间函数分析：分析： 从 测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），和3个水泵不工作时段:第1时段t=0到t=8.97，第2次时段t=10.95到t=20.84第3时段t=23以后对第1、2时段的测量数据直接分别作多项式拟合，得到水位函数为使拟合曲线比较光滑，多项式次数不要太高，一般在36由于第3时段只有3个测量记录，无法对这一时段的水位作出较好的拟合 2. 确定流量确定流量时间函数时间函数 由流量与水位的关系可知，对于第1、2时段只需将水位函数求导数即可得流量; 对于两个供水时段的流量，则

41、用供水时段前后（水泵不工作时段）的流量拟合得到;A. 将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内3、一天总用水量的估计一天总用水量的估计 总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分积分得到。三、模型的求解三、模型的求解1、拟合第、拟合第1、2时段的水位，得基本模型时段的水位，得基本模型2、导出流量关系、导出流量关系, 拟合供水时段的流量拟合供水时段的流量3、估计一天总用水量、估计一天总用水量4、流量及总用水量的检验、流量及总用水量的检验 1、拟合第拟合第1时段的水位，并导出流量时段的水位，并导出流量 设t，h为已输入的时刻和水

42、位测量记录（水泵启动的4个时刻不输入），第第1时段时段各时刻的流量可如下得：1） c1=polyfit（t（1：10），），h（1：10），），3）；）； %用3次多项式拟合第1时段水位，c1输出3次多项式的系数32( )-0.0785t +1.3586 -22.1079967.7356F ttt第第1时段时间和水位的函数关系：时段时间和水位的函数关系：2）a1=polyder（c1）；）； % a1输出多项式（系数为c1）导数的系数 3）tp1=0：0.1：9； x1=-polyval（a1，tp1）；）； % x1输出多项式（系数为a1）在tp1点的函数值（取负后边为正值），即tp1时刻的

43、流量 4）第第1时段流量函数为：时段流量函数为：1079.227173. 22356. 0)(2tttf 2、拟合第拟合第2时段的水位，并导出流量时段的水位，并导出流量 设t，h为已输入的时刻和水位测量记录（水泵启动的4个时刻不输入），第第2时段时段各时刻的流量可如下得：1） c2=polyfit(t(10.9:21),h(10.9:21),3)； %用3次多项式拟合第2时段水位，c2输出3次多项式的系数2） a2=polyder(c2)； % a2输出多项式（系数为c2）导数的系数 3）tp2=10.9:0.1:21; x2=-polyval(a2,tp2); % x2输出多项式（系数为a2

44、）在tp2点的函数值（取负后边为正值），即tp2时刻的流量.4）第第2时段流量函数为：时段流量函数为：8313. 17512. 87529. 00186. 0)(23ttttf 3、拟合供水时段的流量拟合供水时段的流量 在第1供水时段（t=911）之前（即第1时段）和之后（即第2时段）各取几点，其流量已经得到，用它们拟合第1供水时段的流量为使流量函数在t=9和t=11连续，我们简单地只取4个点，拟合3次多项式（即曲线必过这4个点），实现如下： xx1=-polyval（a1，8 9）； %取第1时段在t=8,9的流量 xx2=-polyval（a2，11 12）； %取第2时段在t=11,12

45、的流量 xx12=xx1 xx2； c12=polyfit（8 9 11 12，xx12，3）； %拟合3次多项式 tp12=9：0.1：11； x12=polyval（c12，tp12 % x12输出第1供水时段各时刻的流量拟合的流量函数为：拟合的流量函数为：078.3555879.737207. 3)(2tttf 在第2供水时段之前取t=20，20.8两点的流水量，在该时刻之后（第3时段）仅有3个水位记录，我们用差分得到流量，然后用这4个数值拟合第2供水时段的流量如下： dt3=diff（t(22：24)）； %最后3个时刻的两两之差 dh3=diff（h(22：24)）； %最后3个水位

46、的两两之差 dht3=-dh3./dt3； %t(22)和t(23)的流量 t3=20 20.8 t(22) t(23)； xx3=-polyval(a2，t3(1：2)，dht3)； %取t3各时刻的流量 c3=polyfit（t3，xx3，3）； %拟合3次多项式 t3=20.8：0.1：24； x3=polyval（c3，tp3）；% x3输出第2供水时段 （外推至t=24）各时刻的流量拟合的流量函数为：拟合的流量函数为：8283.913077. 71405. 0)(2tttf 3、一天总用水量的估计一天总用水量的估计 第1、2时段和第1、2供水时段流量的积分之和，就是一天总用水量虽然诸

47、时段的流量已表为多项式函数，积分可以解析地算出，这里仍用数值积分计算如下： y1=0.1*trapz(x1)； %第1时段用水量（仍按高 度计），0.1为积分步长 y2=0.1*trapz(x2)； %第2时段用水量 y12=0.1*trapz(x12)； %第1供水时段用水量 y3=0.1*trapz(x3)； %第2供水时段用水量 y=（y1+y2+y12+y3）*237.8*0.01； %一天总用水量（ ） 计算结果：计算结果：y1=146.2, y2=266.8, y12=47.4, y3=77.3，y=1250.4Lm3310 4、流量及总用水量的检验流量及总用水量的检验 计算出的各

48、时刻的流量各时刻的流量可用水位记录的数值微分来检验用水量y1可用第1时段水位测量记录中下降高度968-822=146来检验，类似地，y2用1082-822=260检验供水时段流量供水时段流量的一种检验方法检验方法如下：供水时段的用水量加上水位上升值260是该时段泵入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量），而两个供水时段水泵的功率应大致相等第1、2时段水泵的功率可计算如下： p1=(y12+260)/2； %第1供水时段水泵的功率 （水量仍以高度计） tp4=20.8：0.1：23； xp2=polyval（c3，tp4）； % xp2输出第2供水时段 各时刻的流量 p2=(0.1*trapz(xp2)+260)/2.2； %第2供水时段水泵的功率 （水量仍以高度计）计算结果计算结果：p1=154.5 ，p2=140.1计算结果计算结果(n1,n2) y1, y2 , y12 , y3 y p1 p2(3,4)146.2, 266.8, 47.4, 77.3 1250.4 154.5 140.1(5,6)146.5, 257.8, 46.1, 76.3 1282.4 15