《Matlab数据拟合实用教程.ppt》由会员分享，可在线阅读，更多相关《Matlab数据拟合实用教程.ppt（44页珍藏版）》请在三一办公上搜索。
1、用Matlab进行数据拟合,1.多项式曲线拟合:polyfit.,y0=polyval(p,x0),p=polyfit(x,y,m),其中,x,y为已知数据点向量,分别表示横,纵坐标,m为拟合多项式的次数,结果返回m次拟合多项式系数,从高次到低次存放在向量p中.,可求得多项式在x0处的值y0.,例1 已知观测数据点如表所示,分别用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(0 1.3-2 16)p3=pol
2、yfit(x,y,3)p6=polyfit(x,y,6),编写Matlab程序如下:,t=0:0.1:1.2s=polyval(p3,t)s1=polyval(p6,t)hold onplot(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(0 1.3-2 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6),例2 用切削机床进行金
3、属品加工时,为了适当地调整机床,需要测定刀具的磨损速度.在一定的时间测量刀具的厚度,得数据如表所示:,解:描出散点图,在命令窗口输入:,t=0:1:16y=30.0 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0plot(t,y,*),解:描出散点图,在命令窗口输入:,t=0:1:16y=30.0 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0plot(t,y,*),a=-
4、0.3012 29.3804,hold on,plot(t,y1),hold off,a=polyfit(t,y,1),y1=-0.3012*t+29.3804,例2 用切削机床进行金属品加工时,为了适当地调整机床,需要测定刀具的磨损速度.在一定的时间测量刀具的厚度,得数据如表所示:,切削时间 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,刀具厚
5、度 y/cm,拟合曲线为:,y=-0.3012t+29.3804,例3 一个15.4cm30.48cm的混凝土柱在加压实验中的应力-应变关系测试点的数据如表所示,1.55,2.47,2.93,3.03,已知应力-应变关系可以用一条指数曲线来描述,即假设,式中,表示应力,单位是 N/m2;表示应变.,2.89,已知应力-应变关系可以用一条指数曲线来描述,即假设,式中,表示应力,单位是 N/m2;表示应变.,解 选取指数函数作拟合时,在拟合前需作变量代换,化为 k1,k2 的线性函数.,于是,令,即,在命令窗口输入:,x=500*1.0e-6 1000*1.0e-6 1500*1.0e-6 200
6、0*1.0e-6 2375*1.0e-6y=3.103*1.0e+3 2.465*1.0e+3 1.953*1.0e+3 1.517*1.0e+3 1.219*1.0e+3z=log(y)a=polyfit(x,z,1)k1=exp(8.3009)w=1.55 2.47 2.93 3.03 2.89plot(x,w,*),y1=exp(8.3009)*x.*exp(-494.5209*x),plot(x,w,*,x,y1,r-),已知应力-应变关系可以用一条指数曲线来描述,即假设,式中,表示应力,单位是 N/m2;表示应变.,拟合曲线为:,令,则,求得,于是,在实际应用中常见的拟合曲线有:,直
7、线,多项式,一般 n=2,3,不宜过高.,双曲线(一支),指数曲线,2.非线性曲线拟合:lsqcurvefit.,功能:,x=lsqcurvefit(fun,x0,xdata,ydata),x,resnorm=lsqcurvefit(fun,x0,xdata,ydata),根据给定的数据 xdata,ydata(对应点的横,纵坐标),按函数文件 fun 给定的函数,以x0为初值作最小二乘拟合,返回函数 fun中的系数向量x和残差的平方和resnorm.,例4 已知观测数据点如表所示,求三个参数 a,b,c的值,使得曲线 f(x)=aex+bx2+cx3 与已知数据点在最小二乘意义上充分接近.,
8、首先编写存储拟合函数的函数文件.,function f=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.2+x(3)*xdata.3,保存为文件 nihehanshu.m,例4 已知观测数据点如表所示,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:
9、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.0022 4.0304 0.9404,re
10、snorm=0.0912,例4 已知观测数据点如表所示,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.94 x3.,此时的残差是:0.0912.,f(x)=3ex+4.03x2+0.94 x3.,拟合函数为:,练习:,1.已知观测数据点如表所示,求用三次多项式进行拟合
11、的曲线方程.,2.已知观测数据点如表所示,求a,b,c的值,使得曲线 f(x)=aex+bsin x+c lnx 与已知数据点在最小二乘意义上充分接近.,插值问题,g 表达式复杂,甚至无表达式,1.分段线性插值,实用插值方法,2.三次样条插值,细木条:样条,输入:节点x0,y0,插值点x(均为数组,长度自定义);输出:插值y(与x同长度数组).,1.分段线性插值:已有程序 y=interp1(x0,y0,x)y=interp1(x0,y0,x,linear),2.三次样条插值:已有程序 y=interp1(x0,y0,x,spline)或 y=spline(x0,y0,x),用Matlab作插
12、值计算,例 5 对 在-1,1上,用n=20的等距分点进行分段线性插值,绘制 f(x)及插值函数的图形.,解 在命令窗口输入:,x=-1:0.1:1y=1./(1+9*x.2)xi=-1:0.1:1yi=interp1(x,y,xi)plot(x,y,r-,xi,yi,*),例 6 对 在-5,5上,用n=11个等距分点作分段线性插值和三次样条插值,用m=21个插值点作图,比较结果.,解 在命令窗口输入:,n=11,m=21x=-5:10/(m-1):5y=1./(1+x.2)z=0*xx0=-5:10/(n-1):5y0=1./(1+x0.2)y1=interp1(x0,y0,x)y2=in
13、terp1(x0,y0,x,spline)x y y1 y2plot(x,z,r,x,y,k:,x,y1,b,x,y2,g)gtext(Piece.-linear.),gtext(Spline),gtext(y=1/(1+x2),0 1.0000 1.0000 1.0000 0.5000 0.8000 0.7500 0.8205 1.0000 0.5000 0.5000 0.5000 1.5000 0.3077 0.3500 0.2973 2.0000 0.2000 0.2000 0.2000 2.5000 0.1379 0.1500 0.1401 3.0000 0.1000 0.1000 0
14、.1000 3.5000 0.0755 0.0794 0.0745 4.0000 0.0588 0.0588 0.0588 4.5000 0.0471 0.0486 0.0484 5.0000 0.0385 0.0385 0.0385,例 6 对 在-5,5上,用n=11个等距分点作分段线性插值和三次样条插值,用m=21个插值点作图,比较结果.,x,y,y1,y2,解 在命令窗口输入:,例 7 在一天24h内,从零点开始每间隔2h测得的环境温度为,12,9,9,10,18,24,28,27,25,20,18,15,13,(单位:),推测在每1s时的温度.并描绘温度曲线.,t=0:2:24T=1
15、2 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-),例 8 在飞机的机翼加工时,由于机翼尺寸很大,通常在图纸上只能标出部分关键点的数据.某型号飞机的机翼上缘轮廓线的部分数据如下:,x 0 4.74 9.05 19 38 57 76 95 114 133,y 0 5.23 8.1 11.97 16.15 17.1 16.34 14.
16、63 12.16 6.69,x 152 171 190,y 7.03 3.99 0,例 8 在飞机的机翼加工时,由于机翼尺寸很大,通常在图纸上只能标出部分关键点的数据.某型号飞机的机翼上缘轮廓线的部分数据如下:,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),例9 天文学家在1914年8月份的7次观测中,测得地球与金星
17、之间距离(单位:m),并取其常用对数值与日期的一组历史数据如下所示,试推断何时金星与地球的距离(单位:m)的对数值为 9.9352.,日期,18 20 22 24 26 28 30,距离对数,9.9618 9.9544 9.9468 9.9391 9.9312 9.9232 9.9150,解 由于对数值 9.9352 位于 24 和 26 两天所对应的对数值之间,所以对上述数据用三次样条插值加细为步长为1的数据:,解 由于对数值 9.9352 位于 24 和 26 两天所对应的对数值之间,所以对上述数据用三次样条插值加细为步长为1的数据:,x=18:2:30y=9.9618 9.9544 9.
18、9468 9.9391 9.9312 9.9232 9.9150 xi=18:1:30yi=interp1(x,y,xi,spline)A=xi;yi,A=18.0000 19.0000 20.0000 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000 27.0000 28.0000 29.0000 30.00009.9618 9.9581 9.9544 9.9506 9.9468 9.9430 9.9391 9.9352 9.9312 9.9272 9.9232 9.9191 9.9150,练习:,1.设 在区间-2,2上用10等分点作为节点,分
19、别用三种插值方法:,(1)计算并输出在该区间的20等分点的函数值.,(2)输出这个函数及两个插值函数的图形.,(3)对输出的数据和图形进行分析.,1.设 在区间-2,2上用10等分点作为节点,分别用三种插值方法:,(1)计算并输出在该区间的20等分点的函数值.,zi=0.0183 0.0387 0.0773 0.1411 0.2369 0.3685 0.5273 0.6980 0.8521 0.9599 1.0000 0.9599 0.8521 0.6980 0.5273 0.3685 0.2369 0.1411 0.0773 0.0387 0.0183,1.设 在区间-2,2上用10等分点作
20、为节点,分别用两种插值方法:,(2)输出这个函数及两个插值函数的图形.,练习:,2.已知某型号飞机的机翼断面下缘轮廓线上的部分数据如表所示:,假设需要得到 x 坐标每改变 0.1 时的 y 坐标,分别用两种插值方法对机翼断面下缘轮廓线上的部分数据加细,并作出插值函数的图形.,例5 给药方案,。,一种新药用于临床之前,必须设计给药方案.在快速静脉注射的给药方式下,所谓给药方案是指,每次注射剂量多大,间隔时间多长.药物进入机体后随血液输送到全身,在这个过程中不断地被吸收,分布,代谢,最终排除体外.药物在血液中的浓度,即单位体积血液中的药物含量,称血药浓度.在最简单的一室模型中,将整个机体看作一个房
21、室,称中心室,室内的血药浓度是均匀的.快速静脉注射后,浓度立即上升;然后逐渐下降.当浓度太低时,达不到预期的治疗效果;血药浓度太高,又可能导致药物中毒或副作用太强.临床上,每种药物有一个最小有效浓度 c1 和一个最大治疗浓度 c2.设计给药方案时,要使血药浓度保持在 c1-c2 之间.设本题所研究药物的最小有效浓度c1=10,最大治疗浓度 c2=25,例5 给药方案,。,显然,要设计给药方案,必须知道给药后血药浓度随时间变化的规律.为此,从实验和理论两方面着手.在实验方面,对某人用快速静脉注射方式一次注入该药物300mg后,在一定时刻 t(小时)采集血样,测得血药浓度c,如表:血药浓度c(t)
22、的测试数据,例5 给药方案,。,近似直线关系,即 c(t)有按负指数规律减少的趋势.,例5 给药方案,1.确定血药浓度的变化规律,假设:,a)药物向体外排除的速率与中心室的血药浓度成正比,比例系数为 k(0),称排除速率.,b)中心室血液容积为常数 V,t=0 瞬时注入药物的剂量为 d,血药浓度立即为,由假设 a),中心室的血药浓度 c(t)应满足微分方程,由假设 b),方程的初始条件为:,求解得:,即血药浓度c(t)按指数规律下降.,2.给药方案设计,简单实用的给药方案是:,每隔一定时间,重复注入固定剂量 D,使血药浓度 c(t)呈周期性变化,并保持在 c1-c2 之间.,c1,c2,2.给
23、药方案设计,简单实用的给药方案是:,每隔一定时间,重复注入固定剂量 D,使血药浓度 c(t)呈周期性变化,并保持在 c1-c2 之间.,为此,初次剂量需加大到 D0.,由式 得到:,显然,当 c1,c2 给定后,要确定给药方案,必须知道参数 V 和 k.,2.由实验数据作曲线拟合以确定参数,问题化为由数据 ti,yi(i=1,8)拟合直线,记,为了用线性最小二乘法拟合 的系数 V 和 k,先取对数得,用Matlab作线性最小二乘法拟合,得到,问题化为由数据 ti,yi(i=1,8)拟合直线,记,为了用线性最小二乘法拟合 的系数 V 和 k,先取对数得,用Matlab作线性最小二乘法拟合,得到,由实验数据 d=300(mg)算出:,拟合曲线为:,3.结论,将 k,V 和给出的 c1=10,c2=25 代入,得:,D0=375.5,D=225.3,给药方案不妨定为:,
