建模方法--插值与拟合专题

1、数学建模方法插值与拟合 文理学院数学系：徐艳文理学院数学系：徐艳 插值与拟合的关系 在工程中，常有这样的问题：给定一批数据点（它可以是设计师给定，也可能是从测量与采样中得到），需确定满足特定要求的曲线或曲面。对这个问题有两种方法。q一种是插值法。要求所求曲线（面）通过所给的所有数据点。q另一种方法是数据拟合（曲线拟合与曲面拟合）。人们设法找出某条光滑曲线，它最佳地拟合数据，但不必要经过所有数据点。内容提纲 1、插值问题 2、用Matlab解插值问题 3、数据拟合 4、用Matlab解曲线拟合问题 5、建模案例：水塔流量估计1、插值问题1.1、一维插值 插值问题的一般提法插值问题的一般提法：已知

2、y = f(x)（该函数未知）在互异的n+1个点x0,x1,x2,xn处的函数值y0, y1,y2, yn,构造一个过n+1个点（xk,yk） k=0,1,2,n的次数不超过n的多项式 y = Ln(x)，（称为插值多项式） 使其满足Ln(xk) = yk ，（称为插值条件）然后用y = Ln(x)作为准确函数y = f(x)的近似值。此方法称为插值法插值法。 Theorem：满足插值条件的次数不超过n的多项式是唯一存在的。0 x1xnx0 x1xnx0y1y0y1y两点一次两点一次(线性线性)插值多项式插值多项式: 101001011yxxxxyxxxxxL三点二次三点二次(抛物抛物)插值多

3、项式插值多项式: 2120210121012002010212yxxxxxxxxyxxxxxxxxyxxxxxxxxxL1.1.1 Lagrange插值法就是满足插值条件的n次多项式 Lagrange插值多项式nixxxxxxxxxxxxxxxxxlniiiiiiniii1 ,0,)()()()()()()(110110jijixlji, 0, 1)(满足)()(0 xlyxLiniin则上式称为Lagrange插值基函数例例1 1、已知数据表、已知数据表解：解：基函数为基函数为x x1 12 2f f ( (x x) )0.950.950.820.82写出写出 f(x) 的线性插值函数的线性

4、插值函数 , 并求并求 f(1.5) 的近似值。的近似值。82. 0, 2;95. 0, 11100yxyxxxxxxxxl2212)(10101121)(0101xxxxxxxl线性插值函数为线性插值函数为08. 113. 0) 1(82. 0)2(95. 0)()()(11001xxxxlyxlyxL且且 f(1.5) L1(1.5) = 0.885。Lagrange插值法的缺点 多数情况下，Lagrange插值法效果是不错的，但随着节点数n的增大，Lagrange多项式的次数也会升高，可能造成插值函数的收敛性和稳定性变差。如龙格（Runge）现象。 例：在-1,1上用n+1个等距节点作插

5、值多项式Ln(x),使得它在节点处的值与函数y = 1/(1+25x2)在对应节点的值相等，当n增大时，插值多项式在区间的中间部分趋于y(x)，但对于满足条件0.728|x|1的x， Ln(x)并不趋于y(x)在对应点的值，而是发生突变，产生剧烈震荡，即Runge现象。现象。1.1.2 分段插值法 图中看到，随着节点的增加，Lagrange插值函数次数越高，插值函数在两端容易产生龙格现象，为了改进高次插值的缺陷，就产生了分段插值。 分段插值基本思想：将被插函数逐段多项式化。 处理过程：将区间a,b 划分： 在每个子段 上构造低次多项式，然后将其拼接在一起作为整个区间a,b 上的插值函数，这样构

6、造出的插值函数称为分段多项式，改进了多项式插值整体性太强的缺点，可以进行局部调整而不会影响整体。bxxan0,1iixx分段线性插值 设插值节点 若 ： 为已知。)(,10iiinxfybaxxxx满足)(xLy 上的分段线性插值。为则称是线性函数上在每个子段,)(,)()(,)2()() 1 (11111baxLyxxxyxxxxyxxxxxLxLyxxyxLiiiiiiiiiiiiii分段线性插值ijijjxil,.,n,iixixxixixixxixixxixixixxxil, 0, 1)(,10,01,11,1,11)(则，其它由此构造插值基函数：xjxj-1xj+1x0 xnnixi

7、liyxL0)()(1.1.3 三次样条插值 分段线性插值虽然改善了高次插值的缺陷，但它的光滑性不高（一阶导数一般不存在），这往往不能满足某些工程设计上的要求，如对于飞机的机翼型线以及船体放样型值线等往往要求有二阶光滑度，即有二阶连续导数，这就导致了三次样条插值函数的提出。三次样条插值问题提出 设在区间a,b上，已给n+1个互不相同的节点 a=x0 x1xn=b以及函数y = f(x)在这些节点的值f(xi)=yi,i=0,1,n.如果分段函数S(x)满足下列条件： （1） S(x)在子区间xi,xi+1的表达式Si(x)都是次数为3的多项式； （2）S(xi) = yi; （3） S(x)在

8、区间a,b上有连续的二阶导数。 就称S(x)为f(x)在点x0，x1，xn的三次样条插值函数. 即 Si(x)=aix3+bix2+cix+di i=0,1,n xix xi+1 （4n个变量） 需要4n个方程 S(xi) = yi i=0,1,n (n+1个方程） S(xi-0)= S(xi+0) i=1,n-1 在xi连续 (n-1个方程） S/(xi-0)= S/(xi+0) i=1,n-1 在xi连续(n-1个方程） S/(xi-0)= S/(xi+0) i=1,n-1 在xi连续(n-1个方程） 再加两个条件:可在边界点x0与xn处给出导数的约束条件，称为边界条件。 （1）S/(x0

9、)= y0/ ,S /(xn)= yn/ (2) S/(x0)= y0/ ,S /(xn)= yn/ (3) S/(x0)= S /(xn)= 0 自然边界条件 （2个方程） 可以证明：满足上述满足上述4n个线性方程组有唯一解个线性方程组有唯一解。三次样条插值问题分析 xyO O第一种（网格节点）：第一种（网格节点）：1.2 二维插值 已知已知 m n个节点个节点 ),2 , 1;,.,2 , 1(),(njmizyxijji 其中其中jiyx ,互不相同，不妨设互不相同，不妨设bxxxam 21dyyycn 21 构造一个二元函数构造一个二元函数),(yxfz 通过全部已知节点通过全部已知节

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

11、 y1)(x2, y2)O O 二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。1.2.1网格节点插值法网格节点插值法最邻近插最邻近插值值 将四个插值点（矩形的四个顶点）处的函数值依次简记为： 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)=f41.2.11.2.1网格节点插值法网格节点插值法分片线性插值分片线性插值插值函数为：jii1ij1jy)xx(xxyyy)yy)(ff ()xx)(ff (f)y, x

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

13、 y2)O O1.2.1网格节点插值法网格节点插值法双线性插值双线性插值1.2.2 散乱数据插值法 在T=a,b c,d上散乱分布n个点。一般采用反距离加权平均法。 基本思想：在非给定数据的点处，定义其函数值由已知数据按与该点距离的远近作加权平均决定，记 则二元函数（曲面）定义为：22)()(kkkyyxxrnkkkknkkkkkrryxWzyxWrzyxF1221/11),(,),(0,),(其中否则 如此定义的曲面是全局相关的，对曲面的任一点作数据计算都要涉及到全体数据，这在大量数据中是很慢的，但因为这种做法思想简单，人们对它进行了种种改进。2.1 2.1 一维插值函数一维插值函数yi=i

14、nterp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ： 线性插值；线性插值；spline ： 三次样条插值；三次样条插值；cubic ： 立方插值。立方插值。缺省时：缺省时： 分段线性插值。分段线性插值。 注意：所有的插值方法都要求注意：所有的插值方法都要求x x是单调的，并且是单调的，并且xi不不能够超过能够超过x的范围。的范围。2 2、用、用MATLABMATLAB解插值计算解插值计算解解 在命令窗口输入在命令窗口输入:例例 1 在一天在一天24h内内, 从零点开始每

15、间隔从零点开始每间隔2h测得的环境温度为测得的环境温度为12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13 C(单位单位: )推测在每推测在每1s时的温度时的温度. 并描绘温度曲线并描绘温度曲线.t=0:2:24T=12 9 9 10 18 24 28 27 25 20 18 15 13plot(t,T,*)ti=0:1/3600:24T1i=interp1(t,T,ti)plot(t,T,*,ti,T1i,r-)T2i=interp1(t,T,ti,spline)plot(t,T,*,ti,T1i,r-,ti,T2i,g-)例例 2 在飞机的机翼

16、加工时在飞机的机翼加工时, 由于机翼尺寸很大由于机翼尺寸很大, 通常在图通常在图纸上只能标出部分关键点的数据纸上只能标出部分关键点的数据. 某型号飞机的机翼上缘某型号飞机的机翼上缘轮廓线的部分数据如下轮廓线的部分数据如下:x 0 4.74 9.05 19 38 57 76 95 114 133y 0 5.23 8.1 11.97 16.15 17.1 16.34 14.63 12.16 6.69x 152 171 190y 7.03 3.99 0 x=0 4.74 9.05 19 38 57 76 95 114 133 152 171 190y=0 5.23 8.1 11.97 16.15 1

17、7.1 16.34 14.63 12.16 9.69 7.03 3.99 0 xi=0:0.001:190yi=interp1(x,y,xi,spline)plot(xi,yi) 要求要求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)被插值点插值方法插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插

18、值cubic cubic 双三次插值双三次插值缺省时缺省时, , 双线性插值双线性插值2.2 2.2 用用MATLABMATLAB作网格节点数据的插作网格节点数据的插值值 插值函数插值函数griddata格式为格式为: cz =griddata（x，y，z，cx，cy，method） 要求要求cxcx取行向量，取行向量，cycy取为列向量取为列向量。被插值点插值方法插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值缺省时缺省时, , 双线性插值双线性插值2.3 2.3 用用MA

19、TLABMATLAB作散点数据的插值计算作散点数据的插值计算例3 山区地形地貌图已知某处山区地形选点测量坐标数据为：x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5y=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6海拔高度数据为：z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 8

20、2 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87山区地形地貌图 程序原始地貌图程序：x=0:.5:5;y=0:.

21、5:6;xx,yy=meshgrid(x,y);z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81 82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 8

22、2 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89 98 99 97 96 98 94 92 87;mesh(xx,yy,z)加密后的地貌图x=0:.5:5;y=0:.5:6;z=89 90 87 85 92 91 96 93 90 87 82 92 96 98 99 95 91 89 86 84 82 84 96 98 95 92 90 88 85 84 83 81 85 80 81

23、82 89 95 96 93 92 89 86 86 82 85 87 98 99 96 97 88 85 82 83 82 85 89 94 95 93 92 91 86 84 88 88 92 93 94 95 89 87 86 83 81 92 92 96 97 98 96 93 95 84 82 81 84 85 85 81 82 80 80 81 85 90 93 95 84 86 81 98 99 98 97 96 95 84 87 80 81 85 82 83 84 87 90 95 86 88 80 82 81 84 85 86 83 82 81 80 82 87 88 89

24、98 99 97 96 98 94 92 87;xi=linspace(0,5,50); %加密横坐标数据到50个yi=linspace(0,6,80); %加密纵坐标数据到80个xii,yii=meshgrid(xi,yi); %生成网格数据zii=interp2(x,y,z,xii,yii,cubic); %插值mesh(xii,yii,zii) %加密后的地貌图山区地形地貌图 结果例4 海底曲面图 例：在某海域测得一些点（x,y)处的水深z由下表给出，在矩形区域(75,200) (-50,150)内画出海底曲面的图形.X129140103.5 88185.5195105Y7.5141.5

25、2314722.5137.5 85.5Z 4868688X157.5107.57781162162117.5Y-6.5-813 56.5-66.584-33.5z 9988949海底曲面图 程序 x=129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5; y=7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5; z= 4 8 6 8 6 8 8 9 9 8 8 9 4 9; plot3(x,y,z,o), hold on %原始数据点 % 插值 c

26、x=75:0.5:200; cy=-70:0.5:150; cz=griddata(x,y,z,cx,cy,cubic);% 三次插值 meshz(cx,cy,cz)海底曲面图 结果曲线拟合的一般提法：已知一组（二维）数据，即平面上曲线拟合的一般提法：已知一组（二维）数据，即平面上 n个点个点（xi,yi) i=1,n, 寻求一个函数（曲线）寻求一个函数（曲线）y=f(x), 使使 f(x) 在某种准则下与所有数据点最为接近，即曲线拟合得在某种准则下与所有数据点最为接近，即曲线拟合得最好。常采用的方法是最小二乘拟合法。最好。常采用的方法是最小二乘拟合法。 +xyy=f(x)(xi,yi)i i

27、 为点为点（xi,yi) 与与曲线曲线 y=f(x) 的距离的距离3、数据拟合 3.1 直线拟合 问题1 对于给定数据点（xi,yi) （i=0,1,2n）求作）求作一次式一次式y=a+bx使得总误差使得总误差 为最小。由求极值方法，使Q达到极小值的a,b应满足niiibxayQ12)(方程。，从而得到拟合直线从中可解得亦即即bxaybayxxbxayxbanxbxaybxaybQaQniiininiiiniiniiniiiiniii,)1(0)(20)1)(20000020000例1 给定一组数据用简单式子表示数据间关系。解：描点发现，数据分布可以用一直线来近似设所求拟合直线为y=a+bx,

28、代人（2）式得xybababa5138. 19392.605138. 19392.60108396998647027587025则拟合直线为解得3.2 多项式拟合 有时所给数据点用直线拟合并不合适，此时考虑用多项式拟合。 问题2 对于给定数据点（xi,yi) （i=0,1,2n）求作求作m(mn)次多项式次多项式 使总误差使总误差 为最小。为最小。mjjjxay0nimjjijixayQ020）（）（即）（得令多元函数极值。多元函数，上述问题为关于mkxyxxamkxxayaQmkaQmjaQnikiinimjkijijnimjkijijikkj, 1 , 0, 1 , 00 2, 1 , 0

29、, 0), 1 , 0(00000（1） 称为最小二乘法的正规方程组。可以证明（1）的解存在唯一，即得 表达式。 01(,.)ma aa20120000231012000001201200000(1 )nnnnmiimiiiiiinnnnnmiiimii iiiiiinnnnnmmmm mmiiimii iiiiiia nax axaxyax axaxaxyxaxaxaxaxyx （1）mjjjxay04、用、用Matlab解最小二乘拟合（多项式拟合）方解最小二乘拟合（多项式拟合）方法法 在线性最小二乘拟合中，用的较多的是多项在线性最小二乘拟合中，用的较多的是多项式拟合。则式拟合。则Matla

30、bMatlab中有现成的函数中有现成的函数 a=polyfit(x0,y0,m)a=polyfit(x0,y0,m) 其中输入参数其中输入参数x0,y0 x0,y0为要拟合的数据，为要拟合的数据，m m为拟合多为拟合多项式的次数，输出参数项式的次数，输出参数a a为拟合多项式为拟合多项式 系数系数 多项式在多项式在x x处的值处的值y y可用下面的函数计算可用下面的函数计算 y=polyval(a,xy=polyval(a,x) )10.mmya xa xa10,.,maaa a10.mmya xa xa10,.,maaa a例例2 2、将下列数据拟合成一个二次方程、将下列数据拟合成一个二次方

31、程 x=19 25 31 38 44;x=19 25 31 38 44; y=19.0 32.3 49.0 73.3 97.8; y=19.0 32.3 49.0 73.3 97.8; ab=polyfit(x,y,2)ab=polyfit(x,y,2) x0=19:0.1:44; x0=19:0.1:44; y0=ab(3)+ab(2)y0=ab(3)+ab(2)* *x0+ab(1)x0+ab(1)* *x0.2;x0.2; plot(x,y,o,x0,y0,r) plot(x,y,o,x0,y0,r) 结果：结果： abab = = 0.0497 0.0193 0.6882 0.0497

32、 0.0193 0.6882例例3 3 某乡镇企业某乡镇企业1990-19961990-1996年的生产利润如下表：年的生产利润如下表： 年份年份 1990 1991 1992 1993 1994 1995 19961990 1991 1992 1993 1994 1995 1996利润（万元）利润（万元） 70 122 144 152 174 196 20270 122 144 152 174 196 202试预测试预测19971997年和年和19981998年的利润。年的利润。解解 作已知数据的的散点图，作已知数据的的散点图，x0=1990 1991 1992 1993 1994 1995

33、 1996;x0=1990 1991 1992 1993 1994 1995 1996;y0=70 122 144 152 174 196 202;y0=70 122 144 152 174 196 202;plot(x0,y0,plot(x0,y0,* *) 发现该乡镇企业的年生产利润几乎直线上升。发现该乡镇企业的年生产利润几乎直线上升。因此，我们可以用因此，我们可以用 作为拟合函数来预测作为拟合函数来预测该乡镇企业未来的年利润。编写程序如下：该乡镇企业未来的年利润。编写程序如下：x0=1990 1991 1992 1993 1994 1995 1996;x0=1990 1991 1992

34、1993 1994 1995 1996;y0=70 122 144 152 174 196 202;y0=70 122 144 152 174 196 202;a=polyfit(x0,y0,1)a=polyfit(x0,y0,1)y97=polyval(a,1997)y97=polyval(a,1997)y98=polyval(a,1998)y98=polyval(a,1998) 求得求得19971997年的生产利润年的生产利润y97=233.4286y97=233.4286，19981998年的生产利年的生产利润润y98=253.9286y98=253.9286。10ya xa401100

35、705. 4,20aa10ya xa估计水塔的流量估计水塔的流量2、解题思路解题思路3、算法设计与编程算法设计与编程1、问题问题 美国某州某居民区有一供居民用水的园柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量通常水泵每天供水一两次，每次约两小时.水塔是一个高12.2米，直径17.4米的正园柱按照设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作表1 是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水

36、流量，及一天的总用水量 表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

37、22.96 23.88 24.99 25.91866 843 822 / / 1059 1035 1018流量估计的解题思路流量估计的解题思路拟合水位拟合水位时间函数时间函数确定流量确定流量时间函数时间函数估计一天总用水量估计一天总用水量 拟合水位拟合水位时间函数时间函数 测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），和3个水泵不工作时段（以下称第1时段t=0到t=8.97，第2次时段t=10.95到t=20.84和第3时段t=23以后）对第1、2时段的测量数据直接分别作多项式拟合，得到水位函数为使拟合曲线比较光滑，多项式次数不要太高，一般在36由于第3时段只有3个测量记

38、录，无法对这一时段的水位作出较好的拟合 2、确定流量确定流量时间函数时间函数 对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内3、一天总用水量的估计一天总用水量的估计 总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到。算法设计与编程算法设计与编程1、拟合第拟合第1、2时段的水位，并导出流量时段的水位，并导出流量2、拟合供水时段的流量拟合供水时段的流量3、估计一天总用水量估计一天总用水量4、流量及总用水量的检验、流量及

39、总用水量的检验 1、拟合第拟合第1时段的水位，并导出流量时段的水位，并导出流量 设t，h为已输入的时刻和水位测量记录（水泵启动的4个时刻不输入），第第1时段时段各时刻的流量可如下得：1） c1=polyfit（t（1：10），），h（1：10），），3）；）； %用3次多项式拟合第1时段水位，c1输出3次多项式的系数2）a1=polyder（c1）；）； % a1输出多项式（系数为c1）导数的系数 3）tp1=0：0.1：9； x1=-polyval（a1，tp1）；）； % x1输出多项式（系数为a1）在tp1点的函数值（取负后边为正值），即tp1时刻的流量 4）流量函数为：流量函数为：10

40、79.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）在tp2点的函数值（取负后边为正值），即tp2时刻

41、的流量4）流量函数为：流量函数为：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的流量 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：