Matlab数据拟合..ppt
用用Matlab进行进行数据拟合数据拟合1.多项式曲线拟合多项式曲线拟合:polyfit.y0=polyval(p,x0)p=polyfit(x,y,m)其中其中,x,y为已知数据点向量为已知数据点向量,分别表示横分别表示横,纵纵坐标坐标,m为拟合多项式的次数为拟合多项式的次数,结果返回结果返回m次次拟合多项式系数拟合多项式系数,从高次到低次存放在向量从高次到低次存放在向量p中中.可求得多项式在可求得多项式在x0处的值处的值y0.例例1 已知观测数据点如表所示已知观测数据点如表所示xy0-0.4470.11.9780.23.280.36.160.47.080.57.340.67.660.79.560.89.480.99.3111.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(0 1.3-2 16)p3=polyfit(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)gridx=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 用切削机床进行金属品加工时用切削机床进行金属品加工时,为了适当地调整为了适当地调整机床机床,需要测定刀具的磨损速度需要测定刀具的磨损速度.在一定的时间测量在一定的时间测量刀具的厚度刀具的厚度,得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm解解:描出散点图描出散点图,在命令窗口输入在命令窗口输入: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=-0.3012 29.3804hold onplot(t,y1),hold offa=polyfit(t,y,1)y1=-0.3012*t+29.3804例例2 用切削机床进行金属品加工时用切削机床进行金属品加工时,为了适当地调整为了适当地调整机床机床,需要测定刀具的磨损速度需要测定刀具的磨损速度.在一定的时间测量在一定的时间测量刀具的厚度刀具的厚度,得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm拟合曲线为拟合曲线为:y=-0.3012t+29.3804例例3 一个一个15.4cm30.48cm的混凝土柱在加压实验中的应的混凝土柱在加压实验中的应力力-应变关系测试点的数据如表所示应变关系测试点的数据如表所示1.552.472.933.03已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述,即假设即假设式中式中,表示应力表示应力,单位是单位是 N/m2;表示应表示应变变.2.89已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述,即假设即假设式中式中,表示应力表示应力,单位是单位是 N/m2;表示应表示应变变.解解 选取选取指数函数指数函数作拟合时作拟合时,在拟合前需作变量代换在拟合前需作变量代换,化为化为 k1,k2 的线性函的线性函数数.于是于是,令令即即在命令窗口输入在命令窗口输入:x=500*1.0e-6 1000*1.0e-6 1500*1.0e-6 2000*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;表示应表示应变变.拟合曲线为拟合曲线为:令令则则求得求得于是于是在实际应用中常见的拟合曲线有在实际应用中常见的拟合曲线有:直线直线多项式多项式一般一般 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 已知观测数据点如表所示已知观测数据点如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三个参数求三个参数 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 已知观测数据点如表所示已知观测数据点如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.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.0022 4.0304 0.9404resnorm=0.0912例例4 已知观测数据点如表所示已知观测数据点如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.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.已知观测数据点如表所示已知观测数据点如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求用三次多项式进行拟合的曲线方程求用三次多项式进行拟合的曲线方程.2.已知观测数据点如表所示已知观测数据点如表所示xy1.617.72.7491.313.14.1189.43.6110.82.334.50.644.9409.13652.436.9求求a,b,c的值的值,使得曲线使得曲线 f(x)=aex+bsin x+c lnx 与已知数与已知数据点在最小二乘意义上充分接近据点在最小二乘意义上充分接近.插值问题插值问题已知已知 n+1+1个节点个节点其中其中互不相同互不相同,不妨不妨设设求任一插值点求任一插值点处的插值处的插值节点可视为由节点可视为由产生产生,g 表达式复杂表达式复杂,甚至无表达式甚至无表达式1.1.分段线性插值分段线性插值xjxj-1xj+1x0 xn实用插实用插值方法值方法机翼下轮廓机翼下轮廓线线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作插值计算作插值计算例例 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=interp1(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.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个插值点作图个插值点作图,比较结果比较结果.xyy1y2解解 在命令窗口输入在命令窗口输入:例例 7 在一天在一天24h内内,从零点开始每间隔从零点开始每间隔2h测得的环境温度为测得的环境温度为12,9,9,10,18,24,28,27,25,20,18,15,13(单位单位:)推测在每推测在每1s时的温度时的温度.并描绘温度曲线并描绘温度曲线.t=0:2:24T=12 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 133y 0 5.23 8.1 11.97 16.15 17.1 16.34 14.63 12.16 6.69x 152 171 190y 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次观测中次观测中,测得地球与金测得地球与金星之间距离星之间距离(单位单位: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.9468 9.9391 9.9312 9.9232 9.9150 xi=18:1:30yi=interp1(x,y,xi,spline)A=xi;yiA=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等分点作为节点等分点作为节点,分别用三种插值方法分别用三种插值方法:(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.01831.设设 在区间在区间-2,2上用上用10等分点作为节点等分点作为节点,分别用两种插值方法分别用两种插值方法:(2)输出这个函数及两个插值函数的图形输出这个函数及两个插值函数的图形.练习练习:2.已知某型号飞机的机翼断面下缘轮廓线上的部分数已知某型号飞机的机翼断面下缘轮廓线上的部分数据如表所示据如表所示:假设需要得到假设需要得到 x 坐标每改变坐标每改变 0.1 时的时的 y 坐标坐标,分别用两分别用两种插值方法对机翼断面下缘轮廓线上的部分数据加细种插值方法对机翼断面下缘轮廓线上的部分数据加细,并作出插值函数的图形并作出插值函数的图形.xy0031.251.772.092.1112.0121.8131.2141.0151.6例例5 给药方案给药方案。一种新一种新药药用于用于临临床之前床之前,必必须设计给药须设计给药方案方案.在快速静脉注射在快速静脉注射的的给药给药方式下方式下,所所谓给药谓给药方案是指方案是指,每次注射每次注射剂剂量多大量多大,间间隔隔时间时间多多长长.药药物物进进入机体后随血液入机体后随血液输输送到全身送到全身,在在这这个个过过程中不断地被程中不断地被吸收吸收,分布分布,代代谢谢,最最终终排除体外排除体外.药药物在血液中的物在血液中的浓浓度度,即即单单位体位体积积血液中的血液中的药药物含量物含量,称血称血药浓药浓度度.在最在最简单简单的一室模型的一室模型中中,将整个机体看作一个房室将整个机体看作一个房室,称中心室称中心室,室内的血室内的血药浓药浓度是度是均匀的均匀的.快速静脉注射后快速静脉注射后,浓浓度立即上升度立即上升;然后逐然后逐渐渐下降下降.当当浓浓度太低度太低时时,达不到达不到预预期的治期的治疗疗效果效果;血血药浓药浓度太高度太高,又可能又可能导导致致药药物中毒或副作用太物中毒或副作用太强强.临临床上床上,每种每种药药物有一个最小有效物有一个最小有效浓浓度度 c1 和一个最大治和一个最大治疗浓疗浓度度 c2.设计给药设计给药方案方案时时,要使血要使血药浓药浓度保持在度保持在 c1-c2 之之间间.设设本本题题所研究所研究药药物的最小有效物的最小有效浓浓度度c1=10,最大治最大治疗浓疗浓度度 c2=25例例5 给药方案给药方案。显显然然,要要设计给药设计给药方案方案,必必须须知道知道给药给药后血后血药浓药浓度随度随时间时间变变化的化的规规律律.为为此此,从从实验实验和理和理论论两方面着手两方面着手.在在实验实验方面方面,对对某人用快速静脉注射方式一次注入某人用快速静脉注射方式一次注入该药该药物物300mg后后,在一定在一定时时刻刻 t(小小时时)采集血采集血样样,测测得血得血药浓药浓度度c 如表如表:血血药浓药浓度度c(t)的的测试测试数据数据t0.250.511.523468c 19.21 18.1515.3614.1012.899.327.455.243.01例例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 之间之间.x0yc1c22.给药方案设计给药方案设计简单实用的给药方案是简单实用的给药方案是:每隔一定时间每隔一定时间 ,重复注入固定剂量重复注入固定剂量 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,给药方案不妨定为给药方案不妨定为:D0=375 mg,D=225 mg,小时小时.
