《建模方法--插值与拟合专题优秀PPT.ppt》由会员分享，可在线阅读，更多相关《建模方法--插值与拟合专题优秀PPT.ppt（66页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、数学建模方法插值与拟合 文理学院数学系：徐艳 xuyan555sina 插值与拟合的关系 在工程中，常有这样的问题：给定一批数据点（它可以是设计师给定，也可能是从测量与采样中得到），需确定满足特定要求的曲线或曲面。对这个问题有两种方法。一种是插值法。要求所求曲线（面）通过所给的全部数据点。另一种方法是数据拟合（曲线拟合与曲面拟合）。人们设法找出某条光滑曲线，它最佳地拟合数据，但不必要经过全部数据点。内容提纲1、插值问题2、用Matlab解插值问题3、数据拟合4、用Matlab解曲线拟合问题5、建模案例：水塔流量估计1、插值问题1.1、一维插值 插值问题的一般提法：已知y=f(x)（该函数未知）
2、在互异的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的多项式是唯一存在的。两点一次两点一次(线性线性)插值多项式插值多项式:三点二次三点二次(抛物抛物)插值多项式插值多项式:1.1.1 Lagrange插值法就是满足插值条件的n次多项式 Lagrange插值多项式上式称为Lagrange插值基函数例例1 1、
3、已知数据表、已知数据表解：解：基函数为基函数为x12f(x)0.950.82写出写出 f(x)的线性插值函数的线性插值函数 ,并求并求 f(1.5)的近似值。的近似值。线性插值函数为线性插值函数为且且 f(1.5)L1(1.5)=0.885。Lagrange插值法的缺点多数状况下，Lagrange插值法效果是不错的，但随着节点数n的增大，Lagrange多项式的次数也会上升，可能造成插值函数的收敛性和稳定性变差。如龙格（Runge）现象。例：在-1,1上用n+1个等距节点作插值多项式Ln(x),使得它在节点处的值与函数y=1/(1+25x2)在对应节点的值相等，当n增大时，插值多项式在区间的中
4、间部分趋于y(x)，但对于满足条件0.728|x|1的x，Ln(x)并不趋于y(x)在对应点的值，而是发生突变，产生猛烈震荡，即Runge现象。1.1.2 分段插值法图中看到，随着节点的增加，Lagrange插值函数次数越高，插值函数在两端简洁产生龙格现象，为了改进高次插值的缺陷，就产生了分段插值。分段插值基本思想：将被插函数逐段多项式化。处理过程：将区间a,b 划分：在每个子段 上构造低次多项式，然后将其拼接在一起作为整个区间a,b 上的插值函数，这样构造出的插值函数称为分段多项式，改进了多项式插值整体性太强的缺点，可以进行局部调整而不会影响整体。分段线性插值设插值节点若 ：分段线性插值xj
5、xj-1xj+1x0 xn1.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)在区间a,b上有
6、连续的二阶导数。就称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)=y0/,S/(xn)=yn/(2)S/(x0)=
7、y0/,S/(xn)=yn/(3)S/(x0)=S/(xn)=0 自然边界条件（2个方程）可以证明：满足上述4n个线性方程组有唯一解。三次样条插值问题分析 xyO O第一种（网格节点）：第一种（网格节点）：1.2 二维插值 已知已知 m n个节点个节点 其中其中互不相同，不妨设互不相同，不妨设 构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，即 yx0 0其次种（散乱节点）：其次种（散乱节点）：已知已知n个节点个节点其中其中互不相同，互不相同，构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，
8、即 留意：最邻近插值一般不连续。具有连续性的最简洁留意：最邻近插值一般不连续。具有连续性的最简洁的插值是分片线性插值。的插值是分片线性插值。x y(x1,y1)(x1,y2)(x2,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)=f4网格节点插值法网
9、格节点插值法分片线性插值分片线性插值插值函数为：其次片(上三角形区域)：(x,y)满足插值函数为：留意：留意：(x,y)当然应当是在插值节点所形成的矩形区当然应当是在插值节点所形成的矩形区域内。明显，分片线性插值函数是连续的；域内。明显，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足 双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O网格节点插值法网格节点插值法双线性
10、插值双线性插值1.2.2 散乱数据插值法在T=a,b c,d上散乱分布n个点。一般接受反距离加权平均法。基本思想：在非给定数据的点处，定义其函数值由已知数据按与该点距离的远近作加权平均确定，记则二元函数（曲面）定义为：如此定义的曲面是全局相关的，对曲面的任一点作数据计算都要涉及到全体数据，这在大量数据中是很慢的，但因为这种做法思想简洁，人们对它进行了种种改进。2.1 2.1 一维插值函数一维插值函数yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值：最邻近插值linear ：线性插值；线
11、性插值；spline ：三次样条插三次样条插值；值；cubic ：立方插值。立方插值。缺省时：缺省时：分段线性插值。分段线性插值。留意：全部的插值方法都要求留意：全部的插值方法都要求x x是单调的，并且是单调的，并且xixi不不能够超过能够超过x x的范围。的范围。2、用、用MATLAB解插值计算解插值计算解解 在吩咐窗口输入在吩咐窗口输入:例例 1 在一天在一天24h内内,从零点起先每间隔从零点起先每间隔2h测得的环境温度为测得的环境温度为12,9,9,10,18,24,28,27,25,20,18,15,13(单位单位:)推想在每推想在每1s时的温度时的温度.并描绘温度曲线并描绘温度曲线.
12、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 在飞机的机翼加工时在飞机的机翼加工时,由于机翼尺寸很大由于机翼尺寸很大,通常在图通常在图纸上只能标出部分关键点的数据纸上只能标出部分关键点的数据.某型号飞机的机翼上缘某型号飞机的机翼上缘轮廓线的部分数据如下轮廓线的部分数据如下:x 0 4.74 9.05
13、 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 17.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可取可取为矩阵
14、，或为矩阵，或x x取取行向量，行向量，y y取为列向量，取为列向量，x,yx,y的值分别不能超出的值分别不能超出x0,y0 x0,y0的范围。的范围。z=interp2(x0,y0,z0,x,y,method)被插值点插值方法插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值缺省时缺省时,双线性插值双线性插值2.2 2.2 用用MATLABMATLAB作网格节点数据的插作网格节点数据的插值值 插值函数插值函数griddata格式为格式为:cz=griddata（x，y，z，c
15、x，cy，method）要求要求cxcx取行向量，取行向量，cycy取为列向量取为列向量。被插值点插值方法插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值缺省时缺省时,双线性插值双线性插值2.3 2.3 用用MATLABMATLAB作散点数据的插值计算作散点数据的插值计算例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海拔高度
16、数据为：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 82 80 80 81 85 90 93
17、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:.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
18、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 98
19、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 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
20、 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;xi=linspace(0,5,50);%加密横坐标数据到50个yi=linspace(0,6,80);%加密纵坐标数据到80个xii,yii=meshgrid(xi,yi);%生成网格数据zii=interp2(x,y
21、,z,xii,yii,cubic);%插值mesh(xii,yii,zii)%加密后的地貌图山区地形地貌图 结果例4 海底曲面图例：在某海疆测得一些点（x,y)处的水深z由下表给出，在矩形区域(75,200)(-50,150)内画出海底曲面的图形.X129140103.5 88185.5195105Y7.5141.52314722.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
22、 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%原始数据点%插值cx=75:0.5:200;cy=-70:0.5:150;cz=griddata(x,y,z,cx,cy,cubic);%三次插值meshz(cx,cy,cz)海底曲面图 结果曲线拟合的一般提法：已知一组（二维）数据，即平面上曲线拟合的一般提法：已知一组（二维）数据，即平面上 n n个点（个点（xi,yi)i=
23、1,n,xi,yi)i=1,n,寻求一个函数（曲线）寻求一个函数（曲线）y=f(x),y=f(x),使使 f(x)f(x)在某种准则下与全部数据点最为接近，即曲线拟在某种准则下与全部数据点最为接近，即曲线拟合得最好。常接受的方法是最小二乘拟合法。合得最好。常接受的方法是最小二乘拟合法。+xyy=f(x)(xi,yi)i i 为点为点（xi,yi)与与曲线曲线 y=f(x)的距离的距离3、数据拟合3.1 直线拟合问题1 对于给定数据点（xi,yi)（i=0,1,2n）求作一次式y=a+bx使得总误差 为最小。由求极值方法，使Q达到微小值的a,b应满足例1 给定一组数据用简洁式子表示数据间关系。解
24、：描点发觉，数据分布可以用始终线来近似设所求拟合直线为y=a+bx,代人（2）式得3.2 多项式拟合有时所给数据点用直线拟合并不合适，此时考虑用多项式拟合。问题2 对于给定数据点（xi,yi)（i=0,1,2n）求作求作m(mn)次多项式次多项式 使总误差使总误差 为最小。为最小。（1）称为最小二乘法的正规方程组。可以证明（1）的解存在唯一，即得 表达式。（1）4、用、用Matlab解最小二乘拟合（多项式拟合）方解最小二乘拟合（多项式拟合）方法法 在在线线性最小二乘性最小二乘拟拟合中，用的合中，用的较较多的是多多的是多项项式式拟拟合。合。则则Matlab中有中有现现成的函数成的函数 a=pol
25、yfit(x0,y0,m)其中其中输输入参数入参数x0,y0为为要要拟拟合的数据，合的数据，m为拟为拟合多合多项项式的次数，式的次数，输输出参数出参数a为拟为拟合多合多项项式式 系数系数 多多项项式在式在x处处的的值值y可用下面的函数可用下面的函数计计算算 y=polyval(a,x)例例2、将下列数据、将下列数据拟拟合成一个二次方程合成一个二次方程 x=19 25 31 38 44;y=19.0 32.3 49.0 73.3 97.8;ab=polyfit(x,y,2)x0=19:0.1:44;y0=ab(3)+ab(2)*x0+ab(1)*x0.2;plot(x,y,o,x0,y0,r)结
26、结果：果：ab=0.0497 0.0193 0.6882例例3 某某乡镇乡镇企企业业1990-1996年的生年的生产产利利润润如下表：如下表：年份年份 1990 1991 1992 1993 1994 1995 1996利利润润（万元）（万元）70 122 144 152 174 196 202试预试预料料1997年和年和1998年的利年的利润润。解解 作已知数据的的散点作已知数据的的散点图图，x0=1990 1991 1992 1993 1994 1995 1996;y0=70 122 144 152 174 196 202;plot(x0,y0,*)发觉该乡镇发觉该乡镇企企业业的年生的年生
27、产产利利润润几乎直几乎直线线上升。上升。因此，我因此，我们们可以用可以用 作作为拟为拟合函数来合函数来预预料料该乡镇该乡镇企企业业将来的年利将来的年利润润。编编写程序如下：写程序如下：x0=1990 1991 1992 1993 1994 1995 1996;y0=70 122 144 152 174 196 202;a=polyfit(x0,y0,1)y97=polyval(a,1997)y98=polyval(a,1998)求得求得1997年的生年的生产产利利润润y97=233.4286，1998年的生年的生产产利利润润y98=253.9286。估计水塔的流量估计水塔的流量2、解题思路解题
28、思路3、算法设计与编程算法设计与编程1、问题问题 美国某州某居民区有一供居民用水的园柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量通常水泵每天供水一两次，每次约两小时.水塔是一个高12.2米，直径17.4米的正园柱依据设计，水塔水位降至约8.2米时，水泵自动启动，水位升到约10.8米时水泵停止工作表1 是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量流量估计的解题思路流量估计的解题思路拟合水位拟合水位时间函
29、数时间函数确定流量确定流量时间函数时间函数估计一天总用水量估计一天总用水量 拟合水位拟合水位时间函数时间函数 测测量量记记录录看看，一一天天有有两两个个供供水水时时段段（以以下下称称第第1供供水水时时段段和和第第2供供水水时时段段），和和3个个水水泵泵不不工工作作时时段段（以以下下称称第第1时时段段t=0到到t=8.97，第第2次次时时段段t=10.95到到t=20.84和和第第3时时段段t=23以以后后）对对第第1、2时时段段的的测测量量数数据据干干脆脆分分别别作作多多项项式式拟拟合合，得得到到水水位位函函数数为为使使拟拟合合曲曲线线比比较较光光滑滑，多多项项式式次次数数不不要要太太高高，一
30、一般般在在36由由于于第第3时时段段只只有有3个个测测量量记记录录，无无法法对对这这一一时时段段的的水水位位作作出出较较好的拟合好的拟合 2、确定流量确定流量时间函数时间函数 对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内3、一天总用水量的估计一天总用水量的估计 总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到。算法设计与编程算法设计与编程1、拟合第拟合第1、2时段的水位，并导出流量时段的水位，并导出流量2、拟合供
31、水时段的流量拟合供水时段的流量3、估计一天总用水量估计一天总用水量4、流量及总用水量的检验、流量及总用水量的检验 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点
32、的函数值（取负后边为正值），即tp1时刻的流量 4）流量函数为：流量函数为：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点的函数值（取负后边为正值），即tp
33、2时刻的流量4）流量函数为：流量函数为：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
34、（c12，tp12）；%x12输出第1供水时段 各时刻的流量拟合的流量函数为：拟合的流量函数为：在第2供水时段之前取t=20，20.8两点的流水量，在该时刻之后（第3时段）仅有3个水位记录，我们用差分得到流量，然后用这4个数值拟合第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，x
35、x3，3）；%拟合3次多项式 t3=20.8：0.1：24；x3=polyval（c3，tp3）；%x3输出第2供水时段 （外推至t=24）各时刻的流量拟合的流量函数为：拟合的流量函数为：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
36、)；%第2供水时段用水量 y=（y1+y2+y12+y3）*237.8*0.01；%一天总用水量（）计算结果：计算结果：y1=146.2,y2=266.8,y12=47.4,y3=77.3，y=1250.4 4、流量及总用水量的检验流量及总用水量的检验 计算出的各各时时刻刻的的流流量量可用水位记录的数值微分来检验用水量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计算结果计算结果流量函数为：流量函数为：流量曲线见图流量曲线见图n=(3,4)n=(5,6)