《数据拟合.ppt》由会员分享，可在线阅读，更多相关《数据拟合.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:1 y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2 plot(x,y,k.,markersize,25) axis(0 1.3
2、-2 16) p3=polyfit(x,y,3) p6=polyfit(x,y,6),编写Matlab程序如下:,t=0:0.1:1.2 s=polyval(p3,t) s1=polyval(p6,t) hold on plot(t,s,r-,linewidth,2) plot(t,s,b-,linewidth,2) grid,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.3,11.2 plot(x,y,k.,markersize,25) axis(0 1.3 -2 16) p3=polyfit(x,y,3) p6=
3、polyfit(x,y,6),例2 用切削机床进行金属品加工时, 为了适当地调整机床, 需要测定刀具的磨损速度. 在一定的时间测量刀具的厚度, 得数据如表所示:,解: 描出散点图, 在命令窗口输入:,t=0:1:16 y=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.0 plot(t,y,*),解: 描出散点图, 在命令窗口输入:,t=0:1:16 y=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
4、.1 25.7 25.3 24.8 24.0 plot(t,y,*),a = -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
5、.3,12,26.1,13,25.7,14,25.3,15,24.8,16,24.0,刀具厚度 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 的线性函数
6、.,于是,令,即,在命令窗口输入:,x=500*1.0e-6 1000*1.0e-6 1500*1.0e-6 2000*1.0e-6 2375*1.0e-6 y=3.103*1.0e+3 2.465*1.0e+3 1.953*1.0e+3 1.517*1.0e+3 1.219*1.0e+3 z=log(y) a=polyfit(x,z,1) k1=exp(8.3009) w=1.55 2.47 2.93 3.03 2.89 plot(x,w,*),y1=exp(8.3009)*x.*exp( -494.5209*x),plot(x,w,*,x,y1,r-),已知应力-应变关系可以用一条指数曲线
7、来描述, 即假设,式中, 表示应力, 单位是 N/m2; 表示应变.,拟合曲线为:,令,则,求得,于是,在实际应用中常见的拟合曲线有:,直线,多项式,一般 n=2, 3, 不宜过高.,双曲线(一支),指数曲线,2. 非线性曲线拟合: lsqcurvefit.,功能:,x=lsqcurvefit(fun, x0, xdata, ydata),x, resnorm=lsqcurvefit(fun, x0, xdata, ydata),根据给定的数据 xdata, ydata (对应点的横, 纵坐标), 按函数文件 fun 给定的函数, 以x0为初值作最小二乘拟合, 返回函数 fun中的系数向量x和
8、残差的平方和resnorm.,例4 已知观测数据点如表所示,求三个参数 a, b, c的值, 使得曲线 f(x)=aex+bx2+cx3 与已知数据点在最小二乘意义上充分接近.,首先编写存储拟合函数的函数文件.,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.
9、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,
10、13.17; x0=0,0,0; x,resnorm=lsqcurvefit(nihehanshu,x0,xdata,ydata),程序运行后显示,x = 3.0022 4.0304 0.9404,resnorm = 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 与已知数据点在最小二乘意义上充分接近.,说明: 最小二乘意义上的最佳
11、拟合函数为,f(x)= 3ex+ 4.03x2 + 0.94 x3.,此时的残差是: 0.0912.,f(x)= 3ex+ 4.03x2 + 0.94 x3.,拟合函数为:,练习:,1. 已知观测数据点如表所示,求用三次多项式进行拟合的曲线方程.,2. 已知观测数据点如表所示,求a, b, c的值, 使得曲线 f(x)=aex+bsin x+c lnx 与已知数据点在最小二乘意义上充分接近.,插值问题,g 表达式复杂,甚至无表达式,1.分段线性插值,实用插值方法,2. 三次样条插值,细木条: 样条,输入: 节点x0, y0, 插值点x (均为数组,长度自定义); 输出: 插值y (与x同长度数
12、组).,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作插值计算,例 5 对 在-1, 1上, 用n=20的等距分点进行分段线性插值, 绘制 f(x)及插值函数的图形.,解 在命令窗口输入:,x=-1:0.1:1 y=1./(1+9*x.2) xi=-1:0.1:1 yi=interp1(x,y,xi) plot(x,y,r-,xi,yi,*),例 6 对 在-5, 5上, 用n=1
13、1个等距分点作分段线性插值和三次样条插值, 用m=21个插值点作图,比较结果.,解 在命令窗口输入:,n=11, m=21 x=-5:10/(m-1):5 y=1./(1+x.2) z=0*x x0=-5:10/(n-1):5 y0=1./(1+x0.2) y1=interp1(x0,y0,x) y2=interp1(x0,y0,x,spline) x y y1 y2 plot(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.500
14、0 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.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个等距分点
15、作分段线性插值和三次样条插值, 用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:24 T=12 9 9 10 18 24 28 27 25 20 18 15 13 plot(t,T,*),ti=0:1/3600:24 T1i=interp1(t,T,ti) plot(t,T,*,ti,T1i,r-),T2i=interp1(t,T,ti,
16、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.63 12.16 6.69,x 152 171 190,y 7.03 3.99 0,例 8 在飞机的机翼加工时, 由于机翼尺寸很大, 通常在图纸上只能标出部分关键点的数据. 某型号飞机的机翼上缘轮廓线的部分数据如下:,x=0 4.74
17、 9.05 19 38 57 76 95 114 133 152 171 190 y=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:190 yi=interp1(x,y,xi,spline) plot(xi,yi),例9 天文学家在1914年8月份的7次观测中, 测得地球与金星之间距离(单位: m), 并取其常用对数值与日期的一组历史数据如下所示, 试推断何时金星与地球的距离(单位: m)的对数值为 9.9352.,日期,18 20 22 24 26 28 30,距离对数,9.9618 9.9
18、544 9.9468 9.9391 9.9312 9.9232 9.9150,解 由于对数值 9.9352 位于 24 和 26 两天所对应的对数值之间, 所以对上述数据用三次样条插值加细为步长为1的数据:,解 由于对数值 9.9352 位于 24 和 26 两天所对应的对数值之间, 所以对上述数据用三次样条插值加细为步长为1的数据:,x=18:2:30 y=9.9618 9.9544 9.9468 9.9391 9.9312 9.9232 9.9150 xi=18:1:30 yi=interp1(x,y,xi,spline) A=xi;yi,A = 18.0000 19.0000 20.00
19、00 21.0000 22.0000 23.0000 24.0000 25.0000 26.0000 27.0000 28.0000 29.0000 30.0000 9.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等分点作为节点, 分别用三种插值方法:,(1) 计算并输出在该区间的20等分点的函数值.,(2) 输出这个函数及两个插值函数的图形.,(3) 对输出的数据和图形进行分析.,1. 设 在区间-2, 2上用10等分
20、点作为节点, 分别用三种插值方法:,(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等分点作为节点, 分别用两种插值方法:,(2) 输出这个函数及两个插值函数的图形.,练习:,2. 已知某型号飞机的机翼断面下缘轮廓线上的部分数据如表所示:,假设需要得到 x 坐
21、标每改变 0.1 时的 y 坐标, 分别用两种插值方法对机翼断面下缘轮廓线上的部分数据加细, 并作出插值函数的图形.,例5 给药方案,。,一种新药用于临床之前, 必须设计给药方案. 在快速静脉注射的给药方式下, 所谓给药方案是指, 每次注射剂量多大, 间隔时间多长. 药物进入机体后随血液输送到全身, 在这个过程中不断地被吸收, 分布, 代谢, 最终排除体外. 药物在血液中的浓度, 即单位体积血液中的药物含量, 称血药浓度. 在最简单的一室模型中, 将整个机体看作一个房室, 称中心室, 室内的血药浓度是均匀的. 快速静脉注射后, 浓度立即上升; 然后逐渐下降. 当浓度太低时, 达不到预期的治疗效
22、果; 血药浓度太高, 又可能导致药物中毒或副作用太强. 临床上, 每种药物有一个最小有效浓度 c1 和一个最大治疗浓度 c2. 设计给药方案时, 要使血药浓度保持在 c1-c2 之间. 设本题所研究药物的最小有效浓度c1=10, 最大治疗浓度 c2=25,例5 给药方案,。,显然, 要设计给药方案, 必须知道给药后血药浓度随时间变化的规律. 为此, 从实验和理论两方面着手. 在实验方面, 对某人用快速静脉注射方式一次注入该药物300mg后, 在一定时刻 t (小时)采集血样, 测得血药浓度c,如表: 血药浓度c(t) 的测试数据,例5 给药方案,。,近似直线关系, 即 c(t) 有按负指数规律
23、减少的趋势.,例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. 给药方案设计,简单实用的给药方案是:,每隔
24、一定时间 , 重复注入固定剂量 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,给药方案不妨定为:,