数学建模里的MATLAB 插值与拟合

1、1数学建模与数学实验数学建模与数学实验后勤工程学院数学教研室 插插 值值2实验目的实验目的实验内容实验内容2、掌握用数学软件包求解插值问题。、掌握用数学软件包求解插值问题。1、了解插值的基本内容。、了解插值的基本内容。11一维插值一维插值22二维插值二维插值33实验作业实验作业3拉格朗日插值拉格朗日插值分段线性插值分段线性插值三次样条插值三次样条插值一一 维维 插插 值值一、一、插值的定义插值的定义二、插值的方法二、插值的方法三、用三、用Matlab解插值问题解插值问题返回返回4返回返回二维插值二维插值一、一、二维插值定义二维插值定义二、网格节点插值法二、网格节点插值法三、用三、用Matlab

2、Matlab解插值问题解插值问题最邻近插值最邻近插值分片线性插值分片线性插值双线性插值双线性插值网格节点数据的插值网格节点数据的插值散点数据的插值散点数据的插值5一维插值的定义一维插值的定义已知已知 n+1个节点个节点, 1 , 0(),(njyxjj其中其中jx互不相同，不妨设互不相同，不妨设),10bxxxan求任一插值点求任一插值点)(*jxx 处的插值处的插值.*y0 x1xnx0y1y节点可视为由节点可视为由)(xgy 产生产生,，g表达式复杂表达式复杂,，或无封闭形式或无封闭形式,，或未知或未知.。*x*y6 构造一个构造一个(相对简单的相对简单的)函数函数),(xfy 通过全部节

3、点通过全部节点, 即即), 1 ,0()(njyxfjj再用再用)(xf计算插值，即计算插值，即).(*xfy 0 x1xnx0y1y*x*y返回返回7 称为拉格朗日插值基函数拉格朗日插值基函数。n0iiiny)x(L)x(P 已知函数f(x)在n+1个点x0,x1,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拉格朗日拉格朗日(

4、Lagrange)插值插值8拉格朗日拉格朗日(Lagrange)插值插值特别地特别地:两点一次两点一次(线性线性)插值多项式插值多项式: 101001011yxxxxyxxxxxL三点二次三点二次(抛物抛物)插值多项式插值多项式: 2120210121012002010212yxxxxxxxxyxxxxxxxxyxxxxxxxxxL .,满足插值条件直接验证可知xLn9 拉格朗日多项式插值的这种振荡现象叫 Runge现象现象55,11)(2xxxg 采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例例返回返回To

5、 MatlabTo Matlablch(larg1)lch(larg1)10分段线性插值分段线性插值其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxL计算量与n无关;n越大，误差越小.nnnxxxxgxL0),()(limxjxj-1xj+1x0 xnxoy11To MATLABxch11，xch12，xch13，xch14返回返回66,11)(2xxxg例例用分段线性插值法求插值用分段线性插值法求插值,并观察插值误差并观察插值误差.1.在在-6,6中平均选取中平均选取5个点作插值个点作插值(xch11)4.在在-6,6中平均选取中平

6、均选取41个点作插值个点作插值(xch14)2.在在-6,6中平均选取中平均选取11个点作插值个点作插值(xch12)3.在在-6,6中平均选取中平均选取21个点作插值个点作插值(xch13)12比分段线性插值更光滑。比分段线性插值更光滑。xyxi-1 xiab 在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。 光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。三次样条插值三次样条插值13 三次样条插值, 1,),()(1nixxxxsxSiii,)()3), 1 ,0()()2), 1()

7、()10223niiiiiiixxCxSniyxSnidxcxbxaxs) 1, 1()()(),()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii自然边界条件）(0)()()40 nxSxS)(,)4)3)2xSdcbaiiii)()(limxgxSng g( (x x) )为被插值函数为被插值函数。14例例66,11)(2xxxg用三次样条插值选取用三次样条插值选取11个基点计算插值个基点计算插值(ych)返回返回To MATLABych(larg1)15用用MATLABMATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi

8、，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ： 线性插值；线性插值；spline ： 三次样条插三次样条插值；值；cubic ： 立方插值。立方插值。缺省时：缺省时： 分段线性插值。分段线性插值。 注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范围。16 例：在例：在1-121-12的的1111小时内，每隔小时内，每隔1 1小时测量一次小时测量一次温度，测得的温度依次为：温度，测得的温度依次为：5 5，8 8，9

9、9，1515，2525，2929，3131，3030，2222，2525，2727，2424。试估计每隔。试估计每隔1/101/10小时的小时的温度值。温度值。To MATLAB(temp)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),ylabel(Degrees Celsius)17xy机翼下轮廓线X035791 11

10、 21 31 41 5Y01 . 21 . 72 . 02 . 12 . 01 . 81 . 21 . 01 . 6例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值。值。To MATLAB(plane)返回返回18二维插值的定义二维插值的定义 xyO O第一种（网格节点）：第一种（网格节点）：19 已知已知 m n个节点个节点 ),2 , 1;,.,2 , 1(),(njmizyxijji 其中其中jiyx ,互不相同，不妨设互不相同，不妨设bxxxam 21dyyycn 21 构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通

11、过全部已知节点,即即再用再用),(yxf计算插值，即计算插值，即).,(*yxfz ),1 ,0;,1 ,0(),(njmizyxfijji 20第二种（散乱节点）：第二种（散乱节点）： yx0 021已知已知n个节点个节点),.,2 , 1(),(nizyxiii 其中其中),(iiyx互不相同，互不相同， 构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节点,即即),1 ,0(),(nizyxfiii 再用再用),(yxf计算插值，即计算插值，即).,(*yxfz 返回返回22 注意：注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值

12、最邻近插值x y(x1, y1)(x1, y2)(x2, y1)(x2, y2)O O 二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。返回返回23 将四个插值点（矩形的四个顶点）处的函数值依次简记为： 分片线性插值分片线性插值xy (xi, yj)(xi, yj+1)(xi+1, yj)(xi+1, yj+1)O Of (xi, yj)=f1，f (xi+1, yj)=f2，f (xi+1, yj+1)=f3，f (xi, yj+1)=f424插值函数为：jii1ij1jy)xx(xxyyy)yy)(ff ()xx)(ff (f)y, x(fj23i121第二片(上三角形

13、区域)：(x, y)满足iii1ij1jy)xx(xxyyy插值函数为：)xx)(ff ()yy)(ff (f)y, x(fi43j141注意注意：(x, y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)： (x, y)满足返回返回25 双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：)dcy)(bax()y, x(f其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值双线性插值x y(x1, y1)(x1, y2)(x2, y1)(x2,

14、y2)O O返回返回26 要求要求x0,y0 x0,y0单调；单调；x x，y y可取可取为矩阵，或为矩阵，或x x取取行向量，行向量，y y取为列向量，取为列向量，x,yx,y的值分别不能超出的值分别不能超出x0,y0 x0,y0的范围。的范围。z=interp2(x0,y0,z0,x,y,method)被插值点插值方法用用MATLAB作网格节点数据的插值作网格节点数据的插值插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值缺省时缺省时, , 双线性插值双线性插值27例：测得平板表

15、面例：测得平板表面3 3* *5 5网格点处的温度分别为：网格点处的温度分别为： 82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 84 82 85 86 84 84 82 85 86 试作出平板表面的温度分布曲面试作出平板表面的温度分布曲面z=f(x,y)z=f(x,y)的图形。的图形。输入以下命令：x=1:5;y=1:3;temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;mesh(x,y,temps)1.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.2以平滑数

16、据,在x、y方向上每隔0.2个单位的地方进行插值.28再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图. To MATLAB(wendu)29例例 山区地貌：山区地貌： 在某山区测得一些地点的高程如下表。平面区域为在某山区测得一些地点的高程如下表。平面区域为 1200=x=4000,1200=y=3600)试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。 X Y12001600200024002

17、800320036004000120011301250128012301040900500700160013201450142014001300700900850200013901500150014009001100106095024001500120011001350145012001150101028001500120011001550160015501380107032001500155016001550160016001600155036001480150015501510143013001200980 通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。To M

18、ATLAB (moutain)返回返回30 插值函数插值函数griddata格式为格式为: cz =griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作散点数据的插值计算作散点数据的插值计算 要求要求cxcx取行向量，取行向量，cycy取为列向量取为列向量。被插值点插值方法插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubic 双三次插值双三次插值v4- Matlab提供的插值方法提供的插值方法缺省时缺省时, , 双线性插值双线性插值31 例例 在某海域测得一些点在某海域测得

19、一些点(x,y)(x,y)处的水深处的水深z z由下由下表给出，船的吃水深度为表给出，船的吃水深度为5 5英尺，在矩形区域（英尺，在矩形区域（7575，200200）* *（-50-50，150150）里的哪些地方船要避免进入。）里的哪些地方船要避免进入。xyz129 140 103.5 88 185.5 195 1057.5 141.5 23 147 22.5 137.5 85.54 8 6 8 6 8 8xyz157.5 107.5 77 81 162 162 117.5-6.5 -81 3 56.5 -66.5 84 -33.59 9 8 8 9 4 932 ) 1( .150,5020

20、0,75. 2hd三次插值法作二维插值在矩形区域. 3 作海底曲面图.1 输入插值基点数据To MATLAB hd1返回返回4.作出水深小于5的海域范围,即z=5的等高线.33实验作业实验作业 山区地貌：山区地貌：在某山区测得一些地点的高程如下表：在某山区测得一些地点的高程如下表：( (平平面区域面区域1200=x=4000,1200=y=3600)1200=x=4000,1200=y=3600)，试作出该山区的试作出该山区的地貌图和等高线图，并对几种插值方法进行比较。地貌图和等高线图，并对几种插值方法进行比较。36003200280024002000160012001480 1500 155

21、0 1510 1430 1300 1200 9801500 1550 1600 1550 1600 1600 1600 15501500 1200 1100 1550 1600 1550 1380 10701500 1200 1100 1350 1450 1200 1150 10101390 1500 1500 1400 900 1100 1060 9501320 1450 1420 1400 1300 700 900 8501130 1250 1280 1230 1040 900 500 700Y/x1200 1600 2000 2400 2800 3200 3600 4000返回返回34数

22、学建模与数学实验数学建模与数学实验后勤工程学院数学教研室 拟拟 合合 35实验目的实验目的实验内容实验内容2 2、掌握用数学软件求解拟合问题。、掌握用数学软件求解拟合问题。1 1、直观了解拟合基本内容。、直观了解拟合基本内容。1 1、拟合问题引例及基本拟合问题引例及基本理论。理论。4 4、实验作业。实验作业。2 2、用数学软件求解拟合问题。用数学软件求解拟合问题。3 3、应用实例应用实例36拟拟 合合2.2.拟合的基本原理拟合的基本原理1. 拟合问题引例拟合问题引例37拟拟 合合 问问 题题 引引 例例 1 1温度温度t(t(0 0C) 20.5 32.7 C) 20.5 32.7 51.0

23、73.0 95.751.0 73.0 95.7电阻电阻R(R( ) 765 826 ) 765 826 873 942 1032873 942 1032已知热敏电阻数已知热敏电阻数据：据：求求60600 0C C时的电阻时的电阻R R。2040608010070080090010001100 设设 R=at+bR=at+ba,ba,b为待定系为待定系数数38拟拟 合合 问问 题题 引引 例例 2 2 t (h) 0.25 0.5 1 t (h) 0.25 0.5 1 1.5 2 3 4 6 81.5 2 3 4 6 8c (c ( g/ml) g/ml) 19.21 18.15 15.36 1

24、4.10 12.89 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.019.32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0(t=0注注射射300300mg)mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).c(t).作半对数坐标系作半对数坐标系(semilogy)下的下的图形图形为待定系数kcectckt,)(002468100101102MATLAB(aa1)MATLAB(aa1)39曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法已知一组（二维）

25、数据，即平面上已知一组（二维）数据，即平面上 n n个点个点（x xi i,y,yi i) ) i=1,n,i=1,n, 寻求一个函数（曲线）寻求一个函数（曲线）y=f(x),y=f(x), 使使 f(x)f(x) 在某种准则下与所有数据点最为接近，即曲线拟合在某种准则下与所有数据点最为接近，即曲线拟合得最好。得最好。 +xyy=f(x)y=f(x)(xi,yi)i i i 为点为点（x xi i,y,yi i) ) 与与曲线曲线 y=f(x) y=f(x) 的距离的距离40拟合与插值的关系拟合与插值的关系 函数插值与曲线拟合都是要根据一组数据构造函数插值与曲线拟合都是要根据一组数据构造一个函

26、数作为近似，由于近似的要求不同，二者的数一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。学方法上是完全不同的。 实例：实例：下面数据是某次实验所得，希望得到X和 f之间的关系？x124791 21 31 51 7f1 .53 .96 .611 .71 5 .61 8 .81 9 .62 0 .62 1 .1MATLAB(cn)MATLAB(cn)问题：问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合数据拟合，又称曲线拟合或曲面拟合。若要求所求曲线（面）通过所给所有数据点

27、，就是插值问题插值问题；41最临近插值、线性插值、样条插值与曲线拟最临近插值、线性插值、样条插值与曲线拟合结果：合结果：0246810121416180510152025已知数据点spline三次多项式插值0246810121416180510152025已知数据点linest三次多项式插值0246810121416180510152025已知数据点nearest三次多项式插值42曲线拟合问题最常用的解法曲线拟合问题最常用的解法线性最小二乘法的线性最小二乘法的基本思路基本思路第一步: :先选定一组函数先选定一组函数 r r1 1(x), r(x), r2 2(x), r(x), rm m(x)

28、, (x), mn,m0)k(0)模型假设模型假设1. 1. 机体看作一个房室，室内血药浓度均匀机体看作一个房室，室内血药浓度均匀一室模型一室模型模型建立模型建立d/c(0) 3得：由假设-kcdtdc 2得：由假设ktevdtc)( 在此，在此，d=300mgd=300mg，t t及及c c（t t）在某些点处的）在某些点处的值见前表，需经拟合求出参数值见前表，需经拟合求出参数k k、v v用线性最小二乘拟合用线性最小二乘拟合c(t)c(t)ktevdtc)()/ln(,ln21vdakacyktvdc)/ln(ln2/,121aedvakatayMATLAB(lihe1)MATLAB(li

29、he1)计算结果：计算结果：)(02.15),/1 (2347. 0lvhkd=300;t=0.25 0.5 1 1.5 2 3 4 6 8;c=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2)程序：程序：用非线性最用非线性最小二乘拟合小二乘拟合c(t)c(t)给药方案给药方案 设设计计cc2c10t 设每次注射剂量D, 间隔时间 血药浓度c(t) 应c1 c(t) c2 初次剂量D0 应加大,0DD给药方案记为：给药方案记为：kecc2112ln1cck

30、2 2、)( ,1220ccDcD1 1、计算结果：计算结果：9 . 3, 3 .225, 5 .3750DD)(4),(225),(3750hmgDmgD给药方案：给药方案：c c1 1=10,c=10,c2 2=25=25k=0.2347k=0.2347v=15.02v=15.0268故可制定给药方案：故可制定给药方案：)(4),(225),(3750hmgDmgD即即: : 首次注射首次注射375mg375mg， 其余每次注射其余每次注射225mg225mg， 注射的间隔时间为注射的间隔时间为4 4小时。小时。69估计水塔的流量估计水塔的流量2 2、解题思路解题思路3 3、算法设计与编程

31、算法设计与编程1 1、问题问题70 某居民区有一供居民用水的园柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量通常水泵每天供水一两次，每次约两小时.水塔是一个高12.2米，直径17.4米的正园柱按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作表1 是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量71 表 1 水位测量记录 （符号/表示水泵启动）时刻(h)水位(cm)0 0.

32、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.91866 843 822 / / 1059 1035 101872流量估计的解题思路流量估计的解题思路拟合水位

33、拟合水位 时间函数时间函数确定流量确定流量 时间函数时间函数估计一天总用水量估计一天总用水量73 拟合水位拟合水位 时间函数时间函数 测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），和3个水泵不工作时段（以下称第1时段t=0到t=8.97，第2次时段t=10.95到t=20.84和第3时段t=23以后）对第1、2时段的测量数据直接分别作多项式拟合，得到水位函数为使拟合曲线比较光滑，多项式次数不要太高，一般在36由于第3时段只有3个测量记录，无法对这一时段的水位作出较好的拟合74 2、确定流量确定流量 时时间函数间函数 对于第1、2时段只需将水位函数求导数即可，对于两个供水时

34、段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内753、一天总用水量的估计一天总用水量的估计 总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到。76算法设计与编程算法设计与编程1 1、拟合第拟合第1 1、2 2时段的水位，并导出流量时段的水位，并导出流量2 2、拟合供水时段的流量拟合供水时段的流量3 3、估计一天总用水量估计一天总用水量4 4、流量及总用水量的检验、流量及总用水量的检验77 1、拟合第拟合第1 1时段的水位，并导出流时段的水位，并导出流量量 设t，h为已输

35、入的时刻和水位测量记录（水泵启动的4个时刻不输入），第第1 1时段时段各时刻的流量可如下得：1） c1=polyfitc1=polyfit（t t（1 1：1010），），h h（1 1：1010），），3 3）；）； %用3次多项式拟合第1时段水位，c1输出3次多项式的系数2）a1=polydera1=polyder（c1c1）；）； % a1输出多项式（系数为c1）导数的系数 3）tp1=0tp1=0：0.10.1：9 9； x1=-polyvalx1=-polyval（a1a1，tp1tp1）；）； % x1输出多项式（系数为a1）在tp1点的函数值（取负后边为正值），即tp1时刻的流量

36、 MATLAB(llgj1)MATLAB(llgj1)4）流量函数为：流量函数为：1079.227173. 22356. 0)(2tttf78 2、拟合第拟合第2 2时段的水位，并导出流时段的水位，并导出流量量 设t，h为已输入的时刻和水位测量记录（水泵启动的4个时刻不输入），第第2 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(a

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

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

39、合第2供水时段的流量如下： dt3=diff（t(22：24)）； %最后3个时刻的两两之差 dh3=diff（h(22：24)）； %最后3个水位的两两之差 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）各时刻的流量MATLAB(llgj4)MATLAB(llgj4)拟合

40、的流量函数为：拟合的流量函数为：8283.913077. 71405. 0)(2tttf81 3、一天总用水量的估计一天总用水量的估计 第1、2时段和第1、2供水时段流量的积分之和，就是一天总用水量虽然诸时段的流量已表为多项式函数，积分可以解析地算出，这里仍用数值积分计算如下： 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.4Lm3310MATLAB(llgjz)MATLAB(llgjz)82 4、流量及总用水量的检验流量及总用水量的检验 计算出的各时刻的流量各时刻的流量可用水位记录的数值微分来检验用水量y1可用第1时段水位测量记录中下降高度968-822=146来检验，类似地，y2用1082-822=260检验供水时段流量供水时段流量的一种检验方法检验方法如下：供水时段的用水量加上水位上升值260是该时段泵入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量），而两个供