用MATLAB进行数据拟合.ppt

拟合,对于情况较复杂的实际问题（因素不易化简，作用机理不详）可直接使用数据组建模，寻找简单的因果变量之间的数量关系，从而对未知的情形作预报。这样组建的模型为拟合模型。拟合模型的组建主要是处理好观测数据的误差，使用数学表达式从数量上近似因果变量之间的关系。拟合模型的组建是通过对有关变量的观测数据的观察、分析和选择恰当的数学表达方式得到的。,拟合,2.拟合的基本原理,1.拟合问题引例,拟合问题引例1,求600C时的电阻R。,设R=at+ba,b为待定系数,拟合问题引例2,求血药浓度随时间的变化规律c(t).,作半对数坐标系(semilogy)下的图形,MATLAB(aa1),曲线拟合问题的提法,已知一组（二维）数据，即平面上n个点（xi,yi)i=1,n,寻求一个函数（曲线）y=f(x),使f(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。,y=f(x),i为点（xi,yi)与曲线y=f(x)的距离,拟合与插值的关系,函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。,实例：下面数据是某次实验所得，希望得到X和f之间的关系？,MATLAB(cn),问题：给定一批数据点，需确定满足特定要求的曲线或曲面,解决方案：,若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。,若要求所求曲线（面）通过所给所有数据点，就是插值问题；,最临近插值、线性插值、样条插值与曲线拟合结果：,曲线拟合问题最常用的解法线性最小二乘法的基本思路,第一步:先选定一组函数r1(x),r2(x),rm(x),mn,令f(x)=a1r1(x)+a2r2(x)+amrm(x)（1）其中a1,a2,am为待定系数。,第二步:确定a1,a2,am的准则（最小二乘准则）：使n个点（xi,yi)与曲线y=f(x)的距离i的平方和最小。,记,问题归结为，求a1,a2,am使J(a1,a2,am)最小。,线性最小二乘法的求解：预备知识,超定方程组：方程个数大于未知量个数的方程组,超定方程一般是不存在解的矛盾方程组。,如果有向量a使得达到最小，则称a为上述超定方程的最小二乘解。,线性最小二乘法的求解,定理：当RTR可逆时，超定方程组（3）存在最小二乘解，且即为方程组RTRa=RTy的解：a=(RTR)-1RTy,所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。,线性最小二乘拟合f(x)=a1r1(x)+amrm(x)中函数r1(x),rm(x)的选取,1.通过机理分析建立数学模型来确定f(x)；,2.将数据(xi,yi)i=1,n作图，通过直观判断确定f(x)：,用MATLAB解拟合问题,1、线性最小二乘拟合,2、非线性最小二乘拟合,用MATLAB作线性最小二乘拟合,1.作多项式f(x)=a1xm+amx+am+1拟合,可利用已有程序:,a=polyfit(x,y,m),2.对超定方程组,3.多项式在x处的值y可用以下命令计算：y=polyval（a，x）,例1对下面一组数据作二次多项式拟合,1）输入以下命令：x=0:0.1:1;y=-0.4471.9783.286.167.087.347.669.569.489.3011.2;R=(x.2)xones(11,1)；A=Ry,MATLAB(zxec1),解法1用解超定方程的方法,2）计算结果：=-9.810820.1293-0.0317,1）输入以下命令：x=0:0.1:1;y=-0.4471.9783.286.167.087.347.669.569.489.3011.2;A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,k+,x,z,r)%作出数据点和拟合曲线的图形,2）计算结果：=-9.810820.1293-0.0317,解法2用多项式拟合的命令,MATLAB(zxec2),例2已知观测数据点如表所示,分别用3次和6次多项式曲线拟合这些数据点.,x=0:0.1:1y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2plot(x,y,k.,markersize,25)axis(01.3-216)p3=polyfit(x,y,3)p6=polyfit(x,y,6),编写Matlab程序如下:,t=0:0.1:1.2s=polyval(p3,t)s1=polyval(p6,t)holdonplot(t,s,r-,linewidth,2)plot(t,s,b-,linewidth,2)grid,x=0:0.1:1y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2plot(x,y,k.,markersize,25)axis(01.3-216)p3=polyfit(x,y,3)p6=polyfit(x,y,6),例3用切削机床进行金属品加工时,为了适当地调整机床,需要测定刀具的磨损速度.在一定的时间测量刀具的厚度,得数据如表所示:,解:描出散点图,在命令窗口输入:,t=0:1:16y=30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0plot(t,y,*),解:描出散点图,在命令窗口输入:,t=0:1:16y=30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0plot(t,y,*),a=-0.301229.3804,holdon,plot(t,y1),holdoff,a=polyfit(t,y,1),y1=-0.3012*t+29.3804,例3用切削机床进行金属品加工时,为了适当地调整机床,需要测定刀具的磨损速度.在一定的时间测量刀具的厚度,得数据如表所示:,切削时间t/h,0,30.0,1,29.1,2,28.4,3,28.1,4,28.0,5,27.7,6,27.5,7,27.2,8,27.0,刀具厚度y/cm,切削时间t/h,9,26.8,10,26.5,11,26.3,12,26.1,13,25.7,14,25.3,15,24.8,16,24.0,刀具厚度y/cm,拟合曲线为:,y=-0.3012t+29.3804,例4一个15.4cm30.48cm的混凝土柱在加压实验中的应力-应变关系测试点的数据如表所示,1.55,2.47,2.93,3.03,已知应力-应变关系可以用一条指数曲线来描述,即假设,式中,表示应力,单位是N/m2;表示应变.,2.89,已知应力-应变关系可以用一条指数曲线来描述,即假设,式中,表示应力,单位是N/m2;表示应变.,解选取指数函数作拟合时,在拟合前需作变量代换,化为k1,k2的线性函数.,于是,令,即,在命令窗口输入:,x=500*1.0e-61000*1.0e-61500*1.0e-62000*1.0e-62375*1.0e-6y=3.103*1.0e+32.465*1.0e+31.953*1.0e+31.517*1.0e+31.219*1.0e+3z=log(y)a=polyfit(x,z,1)k1=exp(8.3009)w=1.552.472.933.032.89plot(x,w,*),y1=exp(8.3009)*x.*exp(-494.5209*x),plot(x,w,*,x,y1,r-),已知应力-应变关系可以用一条指数曲线来描述,即假设,式中,表示应力,单位是N/m2;表示应变.,拟合曲线为:,令,则,求得,于是,在实际应用中常见的拟合曲线有:,直线,多项式,一般n=2,3,不宜过高.,双曲线(一支),指数曲线,二、人口预测线性模型,对于开始提出的实验问题,代如数据，计算得,从而得到人口数与年份的函数关系为,把x=1999代如，估算出1999年的人口数为y=1252.1（百万）12.52亿,1999年实际人口数量为.亿。,线性预测模型,英国人口学家Malthus根据百余年的人口统计资料，于1798年提出了著名的人口自然增长的指数增长模型。,三、人口预测的Malthus模型,基本假设:人口(相对)增长率r是常数,设x(t)为t时刻的人口数量,t=0时人口数为x0,指数增长模型,实际中，常用,1.由前100年的数据求出美国的人口增长Malthus模型。,2.预测后100年（每隔10年）的人口状况。,3.根据预测的人口状况和实际的人口数量,讨论人口模型的改进情况。,例,解：,取得最小值.,其中,表示人口数量。,表示年份,解方程组:,prog41.m%Thisprogramistopredictthenumberofpopulation%formatlongt1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;x2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;lnx1=log(x1);lnx2=log(x2);,a12=sum(t1);a11=10;a21=a12;a22=sum(t1.2);d1=sum(lnx1);d2=sum(lnx1.*t1);A=a11,a12;a21,a22;D=d1;d2;ab=inv(A)*D;disp(a=);disp(ab(1);disp(b=);disp(ab(2);fori=1:10 xx1(i)=exp(ab(1)+ab(2)*t1(i);endfori=1:10 xx2(i)=exp(ab(1)+ab(2)*t2(i);endplot(t1,x1,r*-,t1,xx1,b+-,t2,x2,g*-,t2,xx2,m+-);,a=-49.79535457790735b=0.02859807120038,仿真结果表明：人口增加的指数模型在短期内基本上能比较准确地反映人口自然增长的规律，但长期预测误差很大，需要修正预测模型。,拟合曲线,原始数据曲线,四、人口预测的Logistic模型,人口增长到一定数量后，增长率下降的原因：,资源、环境等因素对人口增长的阻滞作用,且阻滞作用随人口数量增加而变大,假设,r固有增长率(x很小时),k人口容量（资源、环境能容纳的最大数量）,例的Logistic模型留给同学们练习,用多项式拟合人口模型,%Thisprogramistopredictthemodelofpopulationby4-degreepolynomial%prog42.m%formatlongt1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;t=t1;t2;P1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;P2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;P=P1;P2;n=4;%Thedegreeofthefittingpolynomial%a,s=polyfit(t1,P1,n);y=polyval(a,t);%aisthecoefficientsvectorfromn-degreeto0-degree%plot(t,P,r*-,t,y,b+-);,23,a=1.0e+006*-0.000000000000140.00000000107892-0.000003048785950.00381927346813-1.79012132225427,仿真结果表明,人口增加的模型用多项式拟合能比较准确地反映人口自然增长的规律，对长期预测具有指导意义。,例2:海底光缆线长度预测模型,某一通信公司在一次施工中,需要在水面宽为20m的河沟底沿直线走向铺设一条沟底光缆.在铺设光缆之前需要对沟底的地形做初,探测到一组等分点位置的深度数据如下表所示.,25,步探测,从而估计所需光缆的长度,为工程预算提供依据.基本情况如图所示.,(1)预测通过这条河沟所需光缆长度的近似值.,(2)作出铺设沟底光缆的曲线图.,解：用12次多项式函数拟合光缆走势的曲线图如下,仿真结果表明,拟合曲线能较准确地反映光缆的走势图.,ThelengthofthelabelisL=26.3809(m),假设所铺设的光缆足够柔软,在铺设过程中光缆触地走势光滑,紧贴地面,并且忽略水流对光缆的冲击.,%prog45.mThisprogramistofitthedatabypolynomial%formatlongt=linspace(0,20,21);x=linspace(0,20,100);P=9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93;a,s=polyfit(t,P,12);yy=polyval(a,x);disp(yy=);disp(yy);plot(x,yy,r*-,t,P,b+-);L=0;fori=2:100L=L+sqrt(x(i)-x(i-1)2+(yy(i)-yy(i-1)2);enddisp(ThelengthofthelabelisL=);disp(L);,1.lsqcurvefit已知数据点：xdata=（xdata1，xdata2，xdatan），ydata=（ydata1，ydata2，ydatan）,用MATLAB作非线性最小二乘拟合,Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的,可参考例题.,lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），F（x，xdatan）T中的参变量x(向量),使得,输入格式为:(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,);,说明：x=lsqcurvefit(fun,x0,xdata,ydata,options);,lsqnonlin用以求含参量x（向量）的向量值函数f(x)=(f1(x),f2(x),fn(x)T中的参量x，使得最小。其中fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydatai,2.lsqnonlin,已知数据点：xdata=（xdata1，xdata2，xdatan）ydata=（ydata1，ydata2，ydatan）,输入格式为：1）x=lsqnonlin（fun，x0）；2）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）；,例6已知观测数据点如表所示,求三个参数a,b,c的值,使得曲线f(x)=aex+bx2+cx3与已知数据点在最小二乘意义上充分接近.,首先编写存储拟合函数的函数文件.,functionf=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.2+x(3)*xdata.3,保存为文件nihehanshu.m,例6已知观测数据点如表所示,x,y,0,3.1,0.1,3.27,0.2,3.81,0.3,4.5,0.4,5.18,0.5,6,0.6,7.05,0.7,8.56,0.8,9.69,0.9,11.25,1,13.17,求三个参数a,b,c的值,使得曲线f(x)=aex+bx2+cx3与已知数据点在最小二乘意义上充分接近.,编写下面的程序调用拟合函数.,xdata=0:0.1:1;ydata=3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17;x0=0,0,0;x,resnorm=lsqcurvefit(nihehanshu,x0,xdata,ydata),编写下面的程序调用拟合函数.,xdata=0:0.1:1;ydata=3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17;x0=0,0,0;x,resnorm=lsqcurvefit(nihehanshu,x0,xdata,ydata),程序运行后显示,x=3.00224.03040.9404,resnorm=0.0912,例6已知观测数据点如表所示,x,y,0,3.1,0.1,3.27,0.2,3.81,0.3,4.5,0.4,5.18,0.5,6,0.6,7.05,0.7,8.56,0.8,9.69,0.9,11.25,1,13.17,求三个参数a,b,c的值,使得曲线f(x)=aex+bx2+cx3与已知数据点在最小二乘意义上充分接近.,说明:最小二乘意义上的最佳拟合函数为,f(x)=3ex+4.03x2+0.94x3.,此时的残差是:0.0912.,f(x)=3ex+4.03x2+0.94x3.,拟合函数为:,最小二乘拟合：lsqcurvefit()例7假设有一组实测数据xi0.10.20.30.40.50.60.70.80.91.0yi2.32012.6472.97073.28853.60083.9094.21474.51914.82325.1275,假设已知该数据可能满足的原型函数为,试求出满足下面数据的最小二乘解的a,b,c,d的值。,例8用下面一组数据拟合中的参数a，b，k,该问题即解最优化问题：,MATLAB(fzxec1),1）编写M-文件curvefun1.mfunctionf=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.用命令lsqcurvefit,3）运算结果为：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.2542,4）结论：a=0.0063,b=-0.0034,k=0.2542,MATLAB(fzxec2),解法2用命令lsqnonlinf(x)=F(x,tdata,ctada)=x=（a，b，k）,1）编写M-文件curvefun2.mfunctionf=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=lsqnonlin(curvefun2,x0)f=curvefun2(x),函数curvefun2的自变量是x，cdata和tdata是已知参数，故应将cdatatdata的值写在curvefun2.m中,3）