MATLAB应用13.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,Matlab,应用,案例分析,如：,下面数据是某次实验所得，希望得到,X,和,f,之间的关系,问题：,给定一批数据点，需确定满足特定要求的曲线或曲面,插值与拟合,解决方案：,插值与拟合的关系,函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。,若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的,变化趋势,，这就是,数据拟合,，又称曲线拟合或曲面拟合。,若要求所求曲线（面）,通过,所给所有,数据点,，就是,插值,问题,；,插值指令表,interp1,一维插值,interpft,利用,FFT,的一维插值,interp2,二维插值,spline,样条插值,interp3,三维插值,meshgrid,产生经纬矩阵,interpn,高维插值,ndgrid,产生高维经纬矩阵,各种,Matlab,插值指令的用法基本相同，随着维数的增加，插值的具体操作相应复杂，但原理相同。,特别适用于：,由,于实验数据较少无法画出比较好的数据曲线，可以用插值的方法增加数据点。,调用的格式：,yi,=interp1(x,y,xi,method),%,一维插值,(,x,y,),为实验样本数据，其中,x,为数组，必须单调。,y,为对应的函数值。,xi,为拟插入的节点组成的一维数组，,yi,为对应节点的函数值。,method,为指定的插值方法，如果忽略,method,则默认为线性插值。,一维插值,method,四种选择：,nearest,邻近插值,linear,线性插值,spline,立方样条插值,cubic,立方插值,选择方法：,先根据基准数据，调用,matlab,绘图指令（如,plot,，,mesh,等）获得数据的图形表现。然后根据图形，选择插值方法。,例：正弦函数在,0,2,之间取,5,个点。并用各种插值方法增加数据点绘图。,x=0:pi/2:2*pi,%,在正弦曲线上取,5,个数据点,y=,sin(x,),plot(x,y,x,y,o,),xi=0:pi/20:2*pi;,%,指定插入节点数组,xi,线性插值：,yi,=interp1(x,y,xi);,plot(xi,yi,*),邻近插值：,yi,=interp1(x,y,xi,nearest);,plot(xi,yi,*),样条插值：,yi,=interp1(x,y,xi,spline);,plot(xi,yi,*),立方插值：,yi,=interp1(x,y,xi,cubic);,plot(xi,yi,*),四种,method,选择比较,线性,邻近,样条,立方,解：,x=1 2 4 7 9 12 13 15 17;,f=1.5 3.9 6.6 11.7 15.6 18.8 19.6 20.6 21.1;,plot(x,f,x,f,o,);,xi=1:0.5:17;,例：,x,和,f,之间的关系,线性插值：,yi,=interp1(x,f,xi);,plot(xi,yi,*),邻近插值：,yi,=interp1(x,f,xi,nearest);,plot(xi,yi,*),样条插值：,yi,=interp1(x,f,xi,spline);,plot(xi,yi,*),立方插值：,yi,=interp1(x,f,xi,cublic);,plot(xi,yi,*),二维插值,ZI=interp2(X,Y,Z,XI,YI,method),%,二维插值,%,X,Y,Z,是进行插值的基准数据，,X,Y,必须为同维矩阵，矩阵中的元素行列向都必须单调排列。,XI,YI,是待求插补函数值,ZI,的自变量对数组，可以借助,meshgrid,指令形成。,method,为插值方法，四种选项同一维插值。,method,四种选择：,nearest,邻近插值,linear,线性插值,spline,立方样条插值,cubic,立方插值,二维插值,例：有一组海底深度测量数据，采用插值方式绘制海底形状图,（,1,）产生一组海底深度数据,randn(state,2),x=-5:5;y=-5:5;X,Y=meshgrid(x,y);,%,产生经纬矩阵,Z=-500+,1.2*exp(-(X-1).2+(Y-2).2)-0.7*exp(-(X+2).2+(Y+1).2),+randn(size(X)*0.05;,%,产生海底深度矩阵,surf(X,Y,Z);view(-25,25),%,画原始数据图,（,2,）通过插值画更细致的海底图,xi=linspace(-5,5,50);,yi,=linspace(-5,5,50);XI,YI=,meshgrid(xi,yi,);,%,产生插值坐标矩阵,ZI=interp2(X,Y,Z,XI,YI,*cubic);,%,计算插值坐标对应的海底深度,surf(XI,YI,ZI),view(-25,25),%,画插值后的海底深度图,二维插值,原始图,插值后图,样条插值,yy,=,spline(x,y,xx,),%,根据样点,(,x,y,),数据，求,xx,所对应的三次样条插值,yy,pp=,spline(x,y,),%,从样点,(,x,y,),数据，获得逐段多项式样条函数数据,pp,yy,=,ppval(pp,xx,),%,根据逐段多项式样条函数数据,pp,计算,xx,对应的函数值,yy,%,样点数据,(,x,y,),中,x,必须是单调排列的；,若,xx,是行向量，那么,yy,是列维与,y,相同，而行维与,xx,相同的矩阵；,pp,是样条函数的,PP,形式参数，供,ppval,等样条处理指令输入用,样条插值,例：在某一天每隔二小时测量一次气温，具体数据为：,12 12 13 14 16 19 22 26 25 20 18 15 11,，用,spline,函数实现该数据拟合曲线，了解其它时段的气温情况。,程序设计：,调用格式：,yi,=,spline(x,y,xi,),x=0:2:24;t=12 12 13 14 16 19 22 26 25 20 18 15 11;,xi=0:0.5:24;,%,新数据每隔,0.5,小时,yi,=,spline(x,t,xi,);,%,对比,pp=,spline(x,t,);,yy,=,ppval(pp,xi,),plot(x,t,o,xi,yi,),用,MATLAB,解拟合问题,1,、线性最小二乘拟合,2,、非线性最小二乘拟合,线性最小二乘拟合,1.,作多项式,f(x)=a,1,x,m,+,+a,m,x+a,m+1,拟合,:,a,S,=,polyfit(x,y,m,),2.,多项式在,x,处的值,y,计算：,y=,polyval,（,a,，,x,）,ye,delta,=,polyval(a,x,s,),%,利用估计量计算数据带，,ye,原数据的估计，,delta,为离差,输出拟合多项式系数,a=,a,1,a,m,a,m+1,(,数组,),输入同长度,的数组,x,，,y,拟合多项,式次数,S,是一个供,polyval,使用的结构变量，可略,例 对下面一组数据作二次多项式拟合,用多项式拟合：,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),%,计算结果,=-9.8108 20.1293 -0.0317,z=,polyval(A,x,);,plot(x,y,k+,x,z,r,),%,作出数据点和拟合曲线的图形,拟合多项式阶数太低,-,拟合粗糙；,阶次太高,-,噪声放大,非线性最小二乘拟合,Matlab,提供了两个求非线性最小二乘拟合的函数：,lsqcurvefit,和,lsqnonlin,。,两个命令都要先建立,M-,文件,fun.m,，,在其中定义函数,f(x,),。,但两者定义,f(x),的方式不同,1.,lsqcurvefit,格式,:,fun,是一个事先建立的定义函数,F(x,xdata,),的,M-,文件,自变量为,x,和,xdata,x=,lsqcurvefit,(fun,x0,xdata,ydata,options);,迭代初值,已知数据点,选项见无,约束优化,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,);,2.,lsqnonlin,格式,x=,lsqnonlin(fun,x0),；,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=,lsqnonlin,（,fun,，,x0,，,options,）；,fun,是一个事先建立的定义函数,f(x),的,M-,文件，,自变量为,x,迭代初值,选项见无,约束优化,例：用下面一组数据拟合 中的参数,a,，,b,，,k,t,j,100,200,300,400,500,600,700,800,900,1000,c,j,10,3,4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59,解法,1,.,用命令,lsqcurvefit,1,）,编写,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:1000,cdata,=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(1)=a;x(2)=b;x(3)=k,的初值,x=,lsqcurvefit,(curvefun1,x0,tdata,cdata),f=curvefun1(x,tdata),%,运算结果为,：,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.2542,解法,2,用命令,lsqnonlin,1,）,编写,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=x(1)+x(2)*exp(-0.02*x(3)*,tdata,)-,cdata,2,）,输入命令,:,x0=0.2,0.05,0.05;,%,给定,x(1)=a;x(2)=b;x(3)=k,的初值,x=lsqnonlin(curvefun2,x0),f=curvefun2(x),%,运算结果为,f=1.0e-003*(0.2322 -0.1243 -0.2495 -0.2413,-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792);,x=0.0063 -0.0034 0.2542,函数,curvefun2,的自变量是,x,，,cdata,和,tdata,是已知参数，故应将,cdata,tdata,的值写在,curvefun2.m,中,两个命令的拟合结果相同,已知连续,15,天实际日平均气温为：,18 19 17 18 20 22,22,23 21,21,20 22 23,23,22,，求该数据组的拟和方程，并预测以后,5,天的日平均气温。,程序设计：,x=1:15;T=18 19 17 18 20 22,22,23 21,21,20 22 23,23,22;,P=polyfit(x,T,3);,%,拟合三次多项式,Tx,=poly2str(P,x),%,将多项式转换为字符形式,%,Tx,=0.0010474 x3-0.053823 x2+0.96436 x+16.4623,T1=,polyval(P,x,);,%,对应数据点,x,求多项式的值,plot(x,T,x,T1),计算后,5,天的日平均气温：,x=16:20;,Tx,=0.0010474*x.3-0.053823*x.2+0.96436*x+16.4623,Tx,=,22.4035 22.4474 22.4906 22.5392 22.5995,%,得到后,5,天的日平均气温分别为：,22.4,，,22.4,，,22.5,，,22.5,，,22.6,。,应用案例,1,预测日平均气温,城市污水做冷却水的设备腐蚀过程,已知：设备铜管的腐蚀速率,v,与城市二级出水浓缩倍率,k,满足如下的关系式：,v=a,1,k,n,+a,2,k,n-1,+a,3,k,n-2,+,a,n,k+A,（其中,a,1,a,2,a,3,a,n,A,是常数）。,不同浓缩倍率城市二级出水对铜腐蚀的试验数据列表如下：,v/g/(m2,h),0.00362,0.00234,0.00444,0.00571,0.00786,0.00929,k,0,2,4,6,8,10,%,循环水浓缩倍率是循环冷却水运行工况的一个指标。它是指循环水中氯离子浓度与其在补充水中浓度之比，表示循环水在运行中蒸发而使盐类浓缩的倍数。提高循环水的浓缩倍率可减少循环水补充水 用量,求解思路：,用,Matlab,求拟合曲线函数的程序，画图比较分析不同次曲线拟合与实际情况的符合情况，选择适当的拟合函数，给出方程。,应用案例,2,Matlab,求解拟合曲线的程序：,x=0 2 4 6 8 10,y=0.00362 0.00234 0.00444 0.00571 0.00786 0.00929,format long,p3=polyfit(x,y,3),p4=polyfit(x,y,4),p5=polyfit(x,y,5),xcurve,=-1.5:0.5:14,p3curve=polyval(p3,xcurve),p4curve=polyval(p4,xcurve),p5curve=polyval(p5,xcurve),plot(xcurve,p3curve,-,xcurve,p4curve,-.,xcurve,p5curve,-,x,y,*),hold on,lx=0 1.5,ly,=0 0,plot(lx,ly-0.02,-,lx,ly-0.03,-.,lx,ly-0.04,-),text(2,-0.02,degree 3);,%n=3,时的拟合曲线图,text(2,-0.03,degree 4);,%n=4,时的拟合曲线图,text(2,-0.04,degree,5,);,%n=5,时的拟合曲线图,hold off,应用案例,2,%,取,n=4,，拟合方程为：,v=0.00000339843750k,4,-0.00008575810185k,3,+0.00074574652778k,2,-0.00170576719577k+0.00358337301587,某居民区有一个供居民用水的园柱形水塔。,应用案例,3,表,1,是某天的水位测量记录。,试估计任何时刻（包括水泵正供水时）从水塔流出的水的流量以及一天的总用水量,一般可以通过测量其水位来估计水的流量。,但面临的困难是，当水塔水位下降到设定的最低水位时，水泵自动启动向水塔供水，到设定的最高水位时停止供水，这段时间无法测量水塔的水位和水泵的供水量。,通常水泵每天供水一两次，每次约两小时。,水塔是一个高,12.2,米，直径,17.4,米的园柱体。按照设计，水塔水位降至约,8.2,米时，水泵自动启动，水位升到约,10.8,米时水泵停止工作,表,1,、水位测量记录（符号,/,表示水泵启动）,时刻,(h),0,0:55,1:50,2:57,3:52,4:59,5:54,7:01,7:57,8:58,9:59,10:55,10:57,12:02,水位,(cm),968,948,931,913,898,881,869,852,839,822,/,/,1082,1050,时刻,(h),12:57,13:53,14:59,15:54,16:50,17:56,19:02,19:58,20:50,22:01,22:58,23:53,0:59,1:55,水位,(cm),1021,994,965,941,918,892,866,843,822,/,/,1059,1035,1018,应用案例,3,时刻,(h),0,0.92,1.84,2.95,3.87,4.98,5.90,7.01,7.93,8.97,9.98,10.92,10.95,12.03,水位,(cm),968,948,931,913,898,881,869,852,839,822,/,/,1082,1050,时刻,(h),12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,20.84,22.01,22.96,23.88,24.99,25.91,水位,(cm),1021,994,965,941,918,892,866,843,822,/,/,1059,1035,1018,流量估计的解题思路,拟合水位,时间函数,确定流量,时间函数,估计一天总用水量,应用案例,3,1,、拟合水位,时间函数,测量记录看，一天有,2,个供水时段,（以下称第,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,时段只需将,水位函数求导数,即可，对于两个供水时段的流量，则用供水时段前后（水泵不工作时段）的流量拟合得到，并且将拟合得到的第,2,供水时段流量外推，将第,3,时段流量包含在第,2,供水时段内,流量水位关系：,Q=-,h,S,3,、一天总用水量的估计,总用水量等于两个水泵不工作时段和两个供水时段用水量之和，它们都可以由流量对时间的积分得到。,应用案例,3,水位,(cm),968,948,931,913,898,881,869,852,839,822,/,/,1082,1050,水位,(cm),1021,994,965,941,918,892,866,843,822,/,/,1059,1035,1018,算法设计与编程,1,、,拟合第,1,、,2,时段的水位，并导出流量,2,、,拟合供水时段的流量,3,、,估计一天总用水量,4,、,流量及总用水量的检验,应用案例,3,1-1,、拟合第,1,时段的水位，并导出流量,设,t,，,h,为已输入的时刻和水位测量记录（水泵启动的,4,个时刻不输入），,第,1,时段,各时刻的流量可如下得：,1,）,t1=0 0.92 1.84 2.95 3.87 4.98 5.9 7.01 7.93 8.97;,h1=968 948 931 913 898 881 869 852 839 822;,c1=polyfit(t1,h1,3),%,用,3,次多项式拟合第,1,时段水位，,c1,为拟合,h,的多项式系数,2,）,a1=polyder(c1),%a1,为系数,c1,的多项式的导数多项式系数；用于计算第一供水时段不同时刻的水位差,h,3,）,tp1=0:0.1:9;,%0-9,时刻拟合多项式自变量,x1=-polyval(a1,tp1);,%,系数为,a1,多项式在,0-9,函数值，即,x1,为,0-9,时刻供水位差,h,4,）,plot(tp1,x1,r-),%,第,1,时段流量函数为：,应用案例,3,1-2,、拟合第,2,时段的水位，并导出流量,设,t,，,h,为已输入的时刻和水位测量记录（水泵启动的,4,个时刻不输入），,第,2,时段,各时刻的流量可如下得：,1,）,t2=10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.96 20.84;,h2=1082 1050 1021 994 965 941 918 892 866 843 822,c2=polyfit(t2,h2,4);,%,用,4,次多项式拟合第,2,时段水位，,c2,输出,4,次多项式的系数,2,）,a2=polyder(c2),%a2,输出多项式（系数为,c2,）,导数的系数；用于计算第二供水时段不同时刻的水位差,h,3,）,tp2=10.9:0.1:21;,x2=-polyval(a2,tp2);,%x2,为,11-21,时刻供水位差,h,4,）,plot(tp2,x2,g-),%,第,2,时段流量函数为：,f,(t,)=,-0.0284,t,3,+,1.2173,t,2,-,15.9045,t,+34.1994,应用案例,3,2-1,、拟合第,1,供水时段的流量,在第,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);,%,拟合第一次停水时段,h,的,3,次多项式,tp12=9:0.1:11;,x12=polyval(c12,tp12);,%x12,输出第,9-11,时段各时刻的流量,plot(tp12,x12,b-);,%,拟合的流量函数为：,a12=polyder(c12),f(t,)=,-3.5194,t2+,69.6897,t,-336.5049,应用案例,3,2-2,、拟合第,2-3,供水时段的流量,在第,2,供水时段之前取,t=20,，,20.8,两点的流水量，在该时刻之后（第,3,时段）仅有,3,个水位记录，用差分得到流量，然后用这,4,个数值拟合第,2,供水时段的流量。,t3=23.88,24.99,25.91;,dt3=diff(t3);,%,最后,3,个时刻的两两之差,h3=1059,1035,1018;,dh3=diff(h3);,%,最后,3,个水位的两两之差,dht3=-dh3./dt3;,%,计算第,3,供水时段的水位差,h,tt3=20 20.8 t3(1)t3(2);,%,拟合第,2,停水时段供水量用的时间变量,xx3=-polyval(a2,tt3(1:2),dht3,;,%,用第二供水时段拟合公式计算第二停水时段的供水位差,h,，与,第,3,供水时段的水位差,h,组成,第,2-3,供水时段供水,h,c3=polyfit(tt3,xx3,3);,%,拟合第,2-3,供,/,停水时段供水,h,的,3,次多项式,tp3=20.8:0.1:24;,x3=polyval(c3,tp3);,%x3,输出第,21-24,供水时段的供水,h,plot(tp3,x3,y-),a3=polyder(c3),0.2707 -13.2756 159.8758,%,拟合的流量函数为：,f(t,)=,0.2707,t2-13.2756 t+159.8758,应用案例,3,3,、一天总用水量的计算,第,1,、,2,时段和第,1,、,2,供水时段流量的积分之和，就是一天总用水量。诸时段的流量已经表示为多项式函数形式，积分可以解析地算出。这里给出数值积分计算如下：,y1=0.1,*,trapz(x1),%0-9,时段用水,水位差,h,(cm),，,0.1,为积分步长,y2=0.1,*,trapz(x2),%11-21,时段供水水位差,h,cm,y12=0.1,*,trapz(x12),%9-11,时段用水水位差,h,cm,y3=0.1,*,trapz(x3),%21-24,时段用水水位差,h,cm,y=(y1+y2+y12+y3),*,0.01,*17.42*pi/4,%,(s=,d,2,/4,),%y,一天总用水量,m3,计算结果：,y1=146.2,y2=266.4,y12=48.5,y3=78.5,，,y=1283.2,应用案例,3,4,、流量及总用水量的检验,检验方法一：,用水量,y1,可用第,1,时段水位测量记录中下降高度,968-822=146,来检验，类似地，,y2,用,1082-822=260,检验。,检验方法二：,供水时段的用水量加上水位上升值,(,第一供水段为,1082-822=260cm,；,第二供水段为,1059-822=237cm),是该时段泵入的水量，除以时段长度得到水泵的功率（单位时间泵入的水量），而两个供水时段水泵的功率应大致相等第,1,、,2,时段水泵的功率可计算如下：,p1=(y12+260)/(10.92-8.97);,%,第,1,供水时段水泵的功率（水量仍以高度计）,tp4=20.8:0.1:22.96;,xp2=polyval(c3,tp4);,%xp2,输出第,2,供水时段各时刻的流量,p2=(0.1*trapz(xp2)+237)/(22.96-20.8);,%,第,2,供水时段水泵的功率（水量仍以高度计）,计算结果：,p1=158.2,，,p2=134.5,应用案例,3,流量曲线见图,n=(3,4),n=(5,6),应用案例,3,