插值和拟合讲义精.ppt

插值和拟合讲义第1页，本讲稿共77页11一维插值一维插值22二维插值二维插值第2页，本讲稿共77页第3页，本讲稿共77页第4页，本讲稿共77页第5页，本讲稿共77页注：Hermite插值（略）第6页，本讲稿共77页Runge现象第7页，本讲稿共77页第8页，本讲稿共77页第9页，本讲稿共77页用用MATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xi处的插处的插值结果值结果nearest 最邻近插值；最邻近插值；linear 线性插值线性插值(缺省）；缺省）；spline 三次样条插值；三次样条插值；cubic 立方插值；立方插值；缺省时缺省时 分段线性插值分段线性插值 注注意意：所所有有的的插插值值方方法法都都要要求求x是是单单调调，并并且且xi不能够超过不能够超过x的范围的范围第10页，本讲稿共77页例例1.用三次样条插值选取用三次样条插值选取11个基点计算插值个基点计算插值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)此例，可以看出插值函数得到的函数图像与原函数很接近。此例，可以看出插值函数得到的函数图像与原函数很接近。第11页，本讲稿共77页 例例2 2：从：从1 1点点1212点点的的1111小时内，每隔小时内，每隔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,b,hours,temps,r:)%作图xlabel(Hour),ylabel(Degrees Celsius)第12页，本讲稿共77页xy机翼下轮廓线例例3.3.已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值值第13页，本讲稿共77页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=interp1(x0,y0,x);y2=interp1(x0,y0,x,spline);subplot(2,1,1)plot(x0,y0,k+,x,y1,r)gridtitle(piecewise linear)subplot(2,1,2)plot(x0,y0,k+,x,y2,r)gridtitle(spline)第14页，本讲稿共77页第15页，本讲稿共77页第16页，本讲稿共77页第17页，本讲稿共77页第18页，本讲稿共77页 要求要求x0,y0单调；单调；x，y可取为矩阵，或可取为矩阵，或x取行向量，取行向量，y取取为列向量，为列向量，x,y的值分别不能超出的值分别不能超出x0,y0 0的范围的范围z=interp2(x0,y0,z0,x,y,method)被插值点插值方法用用MATLAB作网格节点数据的插值作网格节点数据的插值插值节点被插值点的函数值nearest 最邻近插值；最邻近插值；linear 双线性插值（缺省值）；双线性插值（缺省值）；cubic 双三次插值；双三次插值；缺省时缺省时 双线性插值双线性插值.第19页，本讲稿共77页例例4.4.测得平板表面测得平板表面3 35 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)的图形的图形输入以下命令：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以平滑数据,在 x、y方向上每隔0.2个单位的地方进行插值.第20页，本讲稿共77页再输入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图.第21页，本讲稿共77页 通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较第22页，本讲稿共77页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页，本讲稿共77页 插值函数插值函数griddata格式为格式为:cz=griddata（x，y，z，cx，cy，method）用用MATLAB作散点数据的插值计算作散点数据的插值计算 要求要求cx取行向量，取行向量，cy取为列向量取为列向量被插值点插值方法插值节点被插值点的函数值nearest最邻近插值最邻近插值linear 双线性插值双线性插值cubic 双三次插值双三次插值v4-MATLAB提供的插值方法提供的插值方法缺省时缺省时,双线性插值双线性插值第24页，本讲稿共77页 例例6.6.在某海域测得一些点在某海域测得一些点(x,y)处的水深处的水深z由下表给由下表给出，船的吃水深度为出，船的吃水深度为5 5英尺，在矩形区域（英尺，在矩形区域（7575，200200）（-50-50，150150）里的哪些地方船要避免进入）里的哪些地方船要避免进入第25页，本讲稿共77页x=129 140 103.5 88 185.5 195 105.5 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;cx=75:0.5:200;cy=-50:0.5:150;cz=griddata(x,y,z,cx,cy,cubic);第26页，本讲稿共77页4.作出水深小于5的海域范围,即z=5的等高线.3 作海底曲面图meshz(cx,cy,cz),rotate3dxlabel(X),ylabel(Y),zlabel(Z)figure(2),contour(cx,cy,cz,-5-5);gridhold onplot(x,y,+)xlabel(X),ylabel(Y)第27页，本讲稿共77页用用MATLAB解拟合问题解拟合问题1.1.线性最小二乘拟合线性最小二乘拟合2.2.非线性最小二乘拟合非线性最小二乘拟合第28页，本讲稿共77页第29页，本讲稿共77页第30页，本讲稿共77页第31页，本讲稿共77页第32页，本讲稿共77页第33页，本讲稿共77页第34页，本讲稿共77页用用MATLAB作线性最小二乘拟合作线性最小二乘拟合1.1.作多项式作多项式f(x)=a1xm+amx+am+1拟合拟合,可利用已有程序可利用已有程序:a=polyfit(x,y,m)2.2.对超定方程组对超定方程组可得最小二乘意义下的解可得最小二乘意义下的解，用，用3.3.多项式在多项式在x处的值处的值y可用以下命令计算：可用以下命令计算：y=polyval（a，x）输出拟合多项式系数输出拟合多项式系数a=a1,am,am+1(数组数组)）输入同长度输入同长度的数组的数组x，y拟合多项拟合多项式次数式次数第35页，本讲稿共77页即要求即要求 出二次多项式出二次多项式:中中 的的使得使得:例例7.对下面一组数据作二次多项式拟合对下面一组数据作二次多项式拟合第36页，本讲稿共77页1）输入以下命令）输入以下命令：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=Ry解法解法1 1用解超定方程的方法用解超定方程的方法2）计算结果）计算结果：=-9.8108 20.1293 -0.0317第37页，本讲稿共77页1）输入以下命令：）输入以下命令：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;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,k+,x,z,r)%作出数据点和拟合曲线的图形作出数据点和拟合曲线的图形2）计算结果：）计算结果：=-9.8108 20.1293 -0.0317解法解法2用多项式拟合的命令用多项式拟合的命令第38页，本讲稿共77页1.1.lsqcurvefit已知数据点数据点：xdata=（xdata1,xdata2,xdatan）,ydata=（ydata1,ydata2,ydatan）用用MATLAB作非线性最小二乘拟合作非线性最小二乘拟合 MATLAB提供了两个求非线性最小二乘拟合的函数：提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin两个命令都要先建立两个命令都要先建立M M文件文件fun.m，在其中定义函数，在其中定义函数f(x)，但两者定义，但两者定义f(x)的方式是不同的的方式是不同的,可参考例题可参考例题.用以求含参量用以求含参量x（向量）的向量值函数（向量）的向量值函数F(x,xdata)=（F（x,xdata1）,F（x，xdatan）T中的参变量中的参变量x(x(向量向量),),使得使得 第39页，本讲稿共77页 输入格式为输入格式为:(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是一个事先建立的是一个事先建立的定义函数定义函数F(x,xdata)的的M文件文件,自变量为自变量为x和和xdata说明:x=lsqcurvefit(fun,x0,xdata,ydata,options);迭代初值迭代初值已知数据点已知数据点选项见无选项见无约束优化约束优化第40页，本讲稿共77页 lsqnonlin用以求含参量用以求含参量x（向量）的向量值函数（向量）的向量值函数 f(x)=(=(f1 1(x),),f2 2(x),),fn(x)T T 中的参量中的参量x，使得，使得 最小最小 其中其中 fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydatai 2.lsqnonlin已知数据点：已知数据点：xdata=（xdata1，xdata2，xdatan）ydata=（ydata1，ydata2，ydatan）第41页，本讲稿共77页输入格式为：输入格式为：1)x=lsqnonlin（fun，x0）；2)x=lsqnonlin（fun，x0，options）；3)x=lsqnonlin（fun，x0，optionsgrad）；4)x，options=lsqnonlin（fun，x0，）；5)x，options，funval=lsqnonlin（funx0，）；说明：x=lsqnonlin（fun，x0，options）；fun是一个事先建立的是一个事先建立的定义函数定义函数f(x)的的M文件，文件，自变量为自变量为x迭代初值迭代初值选项见无选项见无约束优化约束优化第42页，本讲稿共77页 例例8.用下面一组数据拟合用下面一组数据拟合 中的参数中的参数a，b，k该问题即解最优化问题：该问题即解最优化问题：第43页，本讲稿共77页 1 1）编写）编写M M文件文件 curvefun1.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: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 1.用命令用命令lsqcurvefit第44页，本讲稿共77页3 3）运算结果为）运算结果为：f=0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 x=0.0063 -0.0034 0.25424）结论）结论：a=0.0063,b=-0.0034,k=0.2542第45页，本讲稿共77页 解法解法 2 用命令用命令lsqnonlin f(x)=F(x,tdata,ctada)=x=（a，b，k）1）编写编写M M文件文件 curvefun2.m function f=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)函数函数curvefun2的自变量是的自变量是x，cdata和和tdata是已知参数，故应将是已知参数，故应将cdata tdata的值写在的值写在curvefun2.m中中第46页，本讲稿共77页3 3）运算结果为）运算结果为f=1.0e-003*(0.2322 -0.1243 -0.2495 -0.2413-0.1668 -=1.0e-003*(0.2322 -0.1243 -0.2495 -0.2413-0.1668 -0.0724 0.0241 0.1159 0.2030 0.27920.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 0.0063 b=-0.0034 =-0.0034 k=0.2542=0.2542第47页，本讲稿共77页MATLAB解应用问题实例解应用问题实例1.电阻问题电阻问题 2.给药方案问题给药方案问题*3.水塔流量估计问题水塔流量估计问题第48页，本讲稿共77页电阻问题电阻问题 温度温度t(C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032例例.由数据由数据拟合拟合R=a1t+b方法方法1.1.用命令用命令polyfit(x,y,m)t=20.5 32.5 51 73 95.7;r=765 826 873 942 1032;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页，本讲稿共77页方法方法2.直接用直接用结果相同结果相同t=20.5 32.5 51 73 95.7;r=765 826 873 942 1032;R=t ones(5,1);aa=Rr;a=aa(1)b=aa(2)y=polyval(aa,t);plot(t,r,k+,t,y,r)第50页，本讲稿共77页一室模型一室模型：将整个机体看作一个房室，称：将整个机体看作一个房室，称中心室中心室，室内血药浓度是，室内血药浓度是均匀的快速静脉注射后，浓度立即上升；然后迅速下降当均匀的快速静脉注射后，浓度立即上升；然后迅速下降当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导致药物中毒或副作用太强临床上，每种药物有一个最小有效致药物中毒或副作用太强临床上，每种药物有一个最小有效浓度浓度c1和一个最大有效浓度和一个最大有效浓度c2设计给药方案时，要使血药浓度设计给药方案时，要使血药浓度 保保持在持在c1c2之间本题设之间本题设c1=10ug/ml，c2=25ug/ml.拟拟 合合 问问 题题 实实 例例 2 2给药方案给药方案 一种新药用于临床之前，必须设计给药方案一种新药用于临床之前，必须设计给药方案.药物进入机体后通过血液输送到全身，在这个过程中不断地被吸收、药物进入机体后通过血液输送到全身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，称为的药物含量，称为血药浓度血药浓度第51页，本讲稿共77页 在实验方面在实验方面,对某人用快速静脉注射方式一次注对某人用快速静脉注射方式一次注入该药物入该药物300mg后后,在一定时刻在一定时刻t(h)采集血药采集血药,测得血测得血药浓度药浓度c(ug/ml)如下表如下表: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 要设计给药方案要设计给药方案,必须知道给药后血药浓度随时间变必须知道给药后血药浓度随时间变化的规律从实验和理论两方面着手：化的规律从实验和理论两方面着手：第52页，本讲稿共77页给药方案给药方案 1.在快速静脉注射的给药方式下，研究血药浓度（单位在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律体积血液中的药物含量）的变化规律tc2cc1O问问题题2.给定药物的最小有效浓度和最大治疗浓度，设计给药方给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长案：每次注射剂量多大；间隔时间多长分分析析 理论：用一室模型研究理论：用一室模型研究血药浓度变化规律血药浓度变化规律 实验：对血药浓度实验：对血药浓度数据作拟合，符合负数据作拟合，符合负指数变化规律指数变化规律第53页，本讲稿共77页3.3.血液容积血液容积v,t=0=0注射剂量注射剂量d,血药浓度立即为血药浓度立即为d/v.2.2.药物排除速率与血药浓度成正比，比例系数药物排除速率与血药浓度成正比，比例系数 k(0)(0)模型假设模型假设1.1.机体看作一个房室，室内血药浓度均匀机体看作一个房室，室内血药浓度均匀一室模型一室模型模型建立模型建立 在此，在此，d=300mg，t及及c（t）在某些点处的值见前表，需经拟）在某些点处的值见前表，需经拟合求出参数合求出参数k、v.第54页，本讲稿共77页用线性最小二乘拟合用线性最小二乘拟合c(t)计算结果：计算结果：d=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)程序：程序：第55页，本讲稿共77页用非线性最小二乘拟合用非线性最小二乘拟合c(t)-用用lsqcurvefit2.主程序如下主程序如下cleartdata=0.25 0.5 1 1.5 2 3 4 6 8;cdata=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;x0=10,0.5;x=lsqcurvefit(curvefun3,x0,tdata,cdata);f=curvefun3(x,tdata)x1.1.用用M M文件文件curvefun3.m定义函数定义函数function f=curvefun3(x,tdata)d=300f=(x(1)d)*exp(-x(2)*tdata)%x(1)=v;x(2)=k第56页，本讲稿共77页给药方案给药方案 设计设计cc2c1Ot 设每次注射剂量D,间隔时间 血药浓度c(t)应c1 c(t)c2 初次剂量D0 应加大给药方案记为：给药方案记为：2.1.计算结果：计算结果：给药方案：给药方案：c1=10,c2=25k=0.2347v=15.02第57页，本讲稿共77页故可制定给药方案：故可制定给药方案：即即:首次注射首次注射375mg，其余每次注射其余每次注射225mg，注射的间隔时间为注射的间隔时间为4h第58页，本讲稿共77页估计水塔的流量估计水塔的流量2.解题思路解题思路3.算法设计与编程算法设计与编程1.问题问题第59页，本讲稿共77页 某居民区有一供居民用水的圆柱形水塔，一般可以通过测量其水位来估计水的流量，但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量通常水泵每天供水一两次，每次约两小时.水塔是一个高12.2m，直径17.4m的正圆柱按照设计，水塔水位降至约8.2m时，水泵自动启动，水位升到约10.8m时水泵停止工作表1 是某一天的水位测量记录，试估计任何时刻（包括水泵正供水时）从水塔流出的水流量，及一天的总用水量第60页，本讲稿共77页第61页，本讲稿共77页流量估计的解题思路流量估计的解题思路拟合水位拟合水位时间函数时间函数确定流量确定流量时间函数时间函数估计一天总用水量估计一天总用水量第62页，本讲稿共77页 拟合水位拟合水位时间函数时间函数 从测量记录看，一天有两个供水时段（以下称第1供水时段和第2供水时段），和3个水泵不工作时段（以下称第1时段t=0到t=8.97，第2次时段t=10.95到t=20.84和第3时段t=23以后）对第1、2时段的测量数据直接分别作多项式拟合，得到水位函数为使拟合曲线比较光滑，多项式次数不要太高，一般在36由于第3时段只有3个测量记录，无法对这一时段的水位作出较好的拟合第63页，本讲稿共77页 确定流量确定流量时间函数时间函数 对于第1、2时段只需将水位函数求导数即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第2供水时段流量外推，将第3时段流量包含在第2供水时段内第64页，本讲稿共77页 一天总用水量的估计一天总用水量的估计 总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到第65页，本讲稿共77页算法设计与编程算法设计与编程1.拟合第拟合第1、2时段的水位，并导出流量时段的水位，并导出流量2.拟合供水时段的流量拟合供水时段的流量3.估计一天总用水量估计一天总用水量4.流量及总用水量的检验流量及总用水量的检验第66页，本讲稿共77页 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）流量函数为：流量函数为：第67页，本讲稿共77页 拟合第拟合第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时刻的流量4）流量函数为：流量函数为：第68页，本讲稿共77页 2.拟合供水时段的流量拟合供水时段的流量 在第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供水时段 各时刻的流量拟合的流量函数为：拟合的流量函数为：第69页，本讲稿共77页 在第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，xx3，3）；）；%拟合3次多项式 t3=20.8：0.1：24；x3=polyval（c3，tp3）；）；%x3输出第2供水时段 （外推至t=24）各时刻的流量拟合的流量函数为：拟合的流量函数为：第70页，本讲稿共77页 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.4第71页，本讲稿共77页 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第72页，本讲稿共77页计算结果计算结果流量函数为：流量函数为：第73页，本讲稿共77页流量曲线见图流量曲线见图n=(3,4)n=(5,6)第74页，本讲稿共77页练习练习1.用给定的多项式，如y=x3-6x2+5x-3，产生一组数据(xi,yi，i=1,2,n),再在yi上添加随机干扰(可用rand产生(0,1)均匀分布随机数,或用rands产生N(0,1)分布随机数)，然后用xi和添加了随机干扰的yi作的3次多项式拟合，与原系数比较 如果作2或4次多项式拟合，结果如何？实验作业实验作业第75页，本讲稿共77页 练习练习2.用电压V=10V的电池给电容器充电，电容器上t时刻的电压为 ，其中V0是电容器的初始电压,是充电常数.试由下面一组t，v数据确定V0,.t(s)0.51234579V(伏）伏）6.366.487.268.228.668.999.439.63实验作业实验作业第76页，本讲稿共77页实验作业实验作业第77页，本讲稿共77页