Matlab数据拟合实用教程学习教案

1、会计学1Matlab数据数据(shj)拟合实用教程拟合实用教程第一页，共44页。例例1 已知观测已知观测(gunc)数据点如数据点如表所示表所示xy0-0.4470.11.9780.23.280.36.160.47.080.57.340.67.660.79.560.89.480.99.3111.2分别用分别用3次和次和6次多项式曲线拟合这些次多项式曲线拟合这些(zhxi)数据数据点点.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

2、 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6)编写编写(binxi)Matlab程程序如下序如下:第1页/共44页第二页，共44页。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=polyf

3、it(x,y,3)p6=polyfit(x,y,6)第2页/共44页第三页，共44页。例例2 用切削用切削(qixio)机床进行金属品加工时机床进行金属品加工时, 为了适当地调整机为了适当地调整机床床, 需要测定刀具的磨损速度需要测定刀具的磨损速度. 在一定的时间测量刀具的厚度在一定的时间测量刀具的厚度, 得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具

4、厚度刀具厚度 y/cm第3页/共44页第四页，共44页。解解: 描出散点图描出散点图, 在命令在命令(mng lng)窗口输入窗口输入: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,*)第4页/共44页第五页，共44页。解解: 描出散点图描出散点图, 在命令窗口在命令窗口(chungku)输入输入: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

5、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第5页/共44页第六页，共44页。例例2 用切削机床进行用切削机床进行(jnxng)金属品加工时金属品加工时, 为了适当地为了适当地调整机床调整机床, 需要测定刀具的磨损速度需要测定刀具的磨损速度. 在一定的时间测量在一定的时间测量刀具的厚度刀具的厚度, 得数据如表所示得数据如表所示:切削切削(qixio)时间时间 t/h030.0129.1228.4328.1

6、428.0527.7627.5727.2827.0刀具刀具(doj)厚厚度度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm拟合曲线为拟合曲线为:y=-0.3012t+29.3804第6页/共44页第七页，共44页。例例3 一个一个(y )15.4cm30.48cm的混凝土柱在加压实验中的应力的混凝土柱在加压实验中的应力-应变关系测试点的数据如表所示应变关系测试点的数据如表所示1.55 2/ N/m 2/ / N/m 6500 10 33.103 10 2.4761000 10 32.46

7、5 10 2. 9361500 10 31.953 10 3. 0362000 10 31.517 10 已知应力已知应力-应变关系可以用一条指数应变关系可以用一条指数(zhsh)曲线来描述曲线来描述, 即假设即假设21 kke式中式中, 表示应力表示应力, 单位是单位是 N/m2; 表示应变表示应变. 2.8962375 10 31.219 10 第7页/共44页第八页，共44页。已知应力已知应力-应变关系可以用一条指数曲线应变关系可以用一条指数曲线(qxin)来描述来描述, 即假设即假设21 kke式中式中, 表示应力表示应力, 单位是单位是 N/m2; 表示应变表示应变. 解解 选取指数

8、函数选取指数函数(zh sh hn sh)作拟合时作拟合时, 在拟合在拟合前需作变量代换前需作变量代换, 化为化为 k1, k2 的线性函数的线性函数(hnsh).于是于是,12lnln kk令令0211ln,ln zakak即即01 zaa第8页/共44页第九页，共44页。在命令窗口在命令窗口(chungku)输入输入: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=po

9、lyfit(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-)第9页/共44页第十页，共44页。已知应力已知应力(yngl)-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述, 即即假设假设21 kke式中式中, 表示应力表示应力, 单位是单位是 N/m2; 表示应变表示应变. 拟合拟合(n h)曲曲线为线为:3 -494.52094.0275 10 e0211ln,ln, zakak01 zaa0211-4

10、94.5209,ln8.3009, akak3124.0275 10 ,494.5209kk令令则则求得求得于是于是(ysh)第10页/共44页第十一页，共44页。在实际应用中常见在实际应用中常见(chn jin)的的拟合曲线有拟合曲线有:01ya xa直线直线(zhxin)101 nnnya xa xa多项式多项式一般一般(ybn) n=2, 3, 不宜不宜过高过高.01ayax双曲线双曲线(一支一支) bxyae指数曲线指数曲线第11页/共44页第十二页，共44页。2. 非线性曲线拟合非线性曲线拟合: lsqcurvefit.功能功能(gngnng):x=lsqcurvefit(fun,

11、x0, xdata, ydata)x, resnorm=lsqcurvefit(fun, x0, xdata, ydata)根据给定的数据根据给定的数据 xdata, ydata (对应点的横对应点的横, 纵坐标纵坐标), 按函数文件按函数文件 fun 给定的函数给定的函数, 以以x0为初值作最小二乘拟合为初值作最小二乘拟合, 返回返回(fnhu)函数函数 fun中的系数向量中的系数向量x和残差的平方和和残差的平方和resnorm.第12页/共44页第十三页，共44页。例例4 已知观测已知观测(gunc)数据点如数据点如表所示表所示xy03.10.13.270.23.810.34.50.45.

12、180.560.67.050.78.560.89.690.911.25113.17求三个参数求三个参数(cnsh) a, b, c的值的值, 使得曲线使得曲线 f(x)=aex+bx2+cx3 与已知数据点在最小二乘意义上充分与已知数据点在最小二乘意义上充分接近接近.首先编写存储拟合函数首先编写存储拟合函数(hnsh)的函的函数数(hnsh)文件文件.function f=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.2+x(3)*xdata.3保存为文件保存为文件 nihehanshu.m第13页/共44页第十四页，共44页。例例4 已知观测

13、已知观测(gunc)数据点如数据点如表所示表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三个参数求三个参数 a, b, c的值的值, 使得曲线使得曲线(qxin) f(x)=aex+bx2+cx3 与已知数据点在最小二乘意义上充分接近与已知数据点在最小二乘意义上充分接近.编写下面编写下面(xi mian)的程序调用的程序调用拟合函数拟合函数.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=

14、0,0,0;x,resnorm=lsqcurvefit(nihehanshu,x0,xdata,ydata)第14页/共44页第十五页，共44页。编写下面的程序调用拟合编写下面的程序调用拟合(n h)函数函数.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)程序运行后显示程序运行后显示(xinsh)x = 3.0022 4.0304 0.9404resnorm = 0.0912第15页/

15、共44页第十六页，共44页。例例4 已知观测已知观测(gunc)数据点如数据点如表所示表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三个参数求三个参数(cnsh) a, b, c的值的值, 使得曲线使得曲线 f(x)=aex+bx2+cx3 与已知数据点在最小二乘意义上充分接近与已知数据点在最小二乘意义上充分接近.说明说明: 最小二乘意义上的最佳拟合最小二乘意义上的最佳拟合(n h)函数函数为为f(x)= 3ex+ 4.03x2 + 0.94 x3.此时的残差是此时的残差是: 0.0912

16、.第16页/共44页第十七页，共44页。f(x)= 3ex+ 4.03x2 + 0.94 x3.拟合拟合(n h)函数为函数为:第17页/共44页第十八页，共44页。练习练习(linx):1. 已知观测已知观测(gunc)数据点数据点如表所示如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求用三次多项式进行拟合求用三次多项式进行拟合(n h)的曲的曲线方程线方程.2. 已知观测数据点如表所示已知观测数据点如表所示xy1.617.72.7491.313.14.1189.43.6110.82.3

17、34.50.644.9409.13652.436.9求求a, b, c的值的值, 使得曲线使得曲线 f(x)=aex+bsin x+c lnx 与已知数据点与已知数据点在最小二乘意义上充分接近在最小二乘意义上充分接近.第18页/共44页第十九页，共44页。插值问题插值问题(wnt)(wnt)已知已知 n+1+1个节点个节点, 1 , 0(),(njyxjj其中其中jx互不相同互不相同, 不妨设不妨设01), na xxxb求任一插值点求任一插值点)(*jxx 处的插值处的插值.*y节点可视为节点可视为由由)(xgy 产生产生,g 表达式复杂表达式复杂(fz),甚至无甚至无表达式表达式0 x1x

18、nx0y1y*x*y(,) (0,1,) jjxyjn第19页/共44页第二十页，共44页。1.1.分段分段(fn (fn dun)dun)线性插值线性插值xjxj-1xj+1x0 xn实用实用(shyn(shyng)g)插值插值方法方法机翼下轮机翼下轮廓线廓线2. 三次三次(sn c)样条插值样条插值细木条细木条: 样条样条第20页/共44页第二十一页，共44页。输入输入: 节点节点(ji din)x0, y0, 插值点插值点x (均为数组均为数组,长长度自定义度自定义);输出输出: 插值插值y (与与x同长度数组同长度数组).1. 分段分段(fn dun)线性插值线性插值: 已有程序已有程

19、序 y=interp1(x0,y0,x) y=interp1(x0,y0,x,linear)2. 三次三次(sn c)样条插值样条插值: 已有程序已有程序 y=interp1(x0,y0,x,spline) 或或 y=spline(x0,y0,x)用用Matlab作插值计算作插值计算第21页/共44页第二十二页，共44页。例例 5 对对 在在-1, 1上上, 用用n=20的等距分的等距分点进行分段线性插值点进行分段线性插值, 绘制绘制 f(x)及插值函数的图形及插值函数的图形. 2119 fxx解解 在命令在命令(mng lng)窗窗口输入口输入:x=-1:0.1:1y=1./(1+9*x.2

20、)xi=-1:0.1:1yi=interp1(x,y,xi)plot(x,y,r-,xi,yi,*)第22页/共44页第二十三页，共44页。例例 6 对对 在在-5, 5上上, 用用n=11个等距分点作分段线个等距分点作分段线性插值和三次样条插值性插值和三次样条插值, 用用m=21个插值点作图个插值点作图,比较结果比较结果.211 yx解解 在命令在命令(mng lng)窗口输入窗口输入: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

21、,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)第23页/共44页第二十四页，共44页。 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.1

22、000 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个插值点作图个插值点作图,比较结果比较结果.211 yxxyy1y2第24页/共44页第二十五页，共44页。解解 在命令窗口在命令窗口(chungku)输入输入:例例 7 在一天在一天24h内内, 从

23、零点开始从零点开始(kish)每间隔每间隔2h测得的环测得的环境温度为境温度为12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13 C(单位单位(dnwi): )推测在每推测在每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,T

24、2i,g-)第25页/共44页第二十六页，共44页。例例 8 在飞机的机翼加工时在飞机的机翼加工时, 由于机翼尺寸很大由于机翼尺寸很大, 通常在图纸上只通常在图纸上只能标出部分能标出部分(b fen)关键点的数据关键点的数据. 某型号飞机的机翼上缘轮廓某型号飞机的机翼上缘轮廓线的部分线的部分(b fen)数据如下数据如下: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第26页/共44页第二十七页，共44页。例例 8

25、 在飞机的机翼加工时在飞机的机翼加工时, 由于机翼尺寸很大由于机翼尺寸很大, 通常在图纸通常在图纸上只能标出部分上只能标出部分(b fen)关键点的数据关键点的数据. 某型号飞机的机翼上某型号飞机的机翼上缘轮廓线的部分缘轮廓线的部分(b fen)数据如下数据如下: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)第27页/

26、共44页第二十八页，共44页。例例9 天文学家在天文学家在1914年年8月份月份(yufn)的的7次观测中次观测中, 测得地球与金星测得地球与金星之间距离之间距离(单位单位: m), 并取其常用对数值与日期的一组历史数据如下并取其常用对数值与日期的一组历史数据如下所示所示, 试推断何时金星与地球的距离试推断何时金星与地球的距离(单位单位: m)的对数值为的对数值为 9.9352.日期日期(rq)18 20 22 24 26 28 30距离距离(jl)对数对数9.9618 9.9544 9.9468 9.9391 9.9312 9.9232 9.9150解解 由于对数值由于对数值 9.9352

27、位于位于 24 和和 26 两天所对应的对数值两天所对应的对数值之间之间, 所以对上述数据用三次样条插值加细为步长为所以对上述数据用三次样条插值加细为步长为1的的数据数据:第28页/共44页第二十九页，共44页。解解 由于对数值由于对数值 9.9352 位于位于 24 和和 26 两天所对应的对数值之两天所对应的对数值之间间, 所以对上述数据所以对上述数据(shj)用三次样条插值加细为步长为用三次样条插值加细为步长为1的数据的数据(shj):x=18:2:30y=9.9618 9.9544 9.9468 9.9391 9.9312 9.9232 9.9150 xi=18:1:30yi=inte

28、rp1(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第29页/共44页第三十页，共44页。练习练习(linx):1. 设设 在区间在区间-2, 2上用上用10等分点作为节点等分点作为节点, 分分别用三种别用三种

29、(sn zhn)插值方法插值方法: 2, xf xe(1) 计算并输出计算并输出(shch)在该区间的在该区间的20等分点的等分点的函数值函数值.(2) 输出这个函数及两个插值函数的图形输出这个函数及两个插值函数的图形.(3) 对输出的数据和图形进行分析对输出的数据和图形进行分析.第30页/共44页第三十一页，共44页。1. 设设 在区间在区间-2, 2上用上用10等分点作为等分点作为(zuwi)节节点点, 分别用三种插值方法分别用三种插值方法: 2, xf xe(1) 计算并输出在该区间计算并输出在该区间(q jin)的的20等分点的等分点的函数值函数值.zi = 0.0183 0.0387

30、 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第31页/共44页第三十二页，共44页。1. 设设 在区间在区间-2, 2上用上用10等分点作为节点等分点作为节点, 分别分别(fnbi)用两种插值方法用两种插值方法: 2, xf xe(2) 输出这个输出这个(zh ge)函数及两个插值函数函数及两个插值函数的图形的图形.第32页/共44页第三十三页，共44页。练习练习(linx)

31、:2. 已知某型号飞机已知某型号飞机(fij)的机翼断面下缘轮廓线上的部分数的机翼断面下缘轮廓线上的部分数据如表所示据如表所示:假设需要得到假设需要得到 x 坐标每改变坐标每改变 0.1 时的时的 y 坐标坐标, 分别用两种插值分别用两种插值方法对机翼断面下缘轮廓线上的部分方法对机翼断面下缘轮廓线上的部分(b fen)数据加细数据加细, 并作并作出插值函数的图形出插值函数的图形.xy0031.251.772.092.1112.0121.8131.2141.0151.6第33页/共44页第三十四页，共44页。例例5 给药方案给药方案(fng n) /g ml 一种新药用于临床之前一种新药用于临床

32、之前, 必须设计给药方案必须设计给药方案. 在快速静脉注射在快速静脉注射的给药方式下的给药方式下, 所谓给药方案是指所谓给药方案是指, 每次注射剂量多大每次注射剂量多大, 间隔间隔时间多长时间多长.药物进入机体后随血液输送到全身药物进入机体后随血液输送到全身, 在这个过程中不断地被在这个过程中不断地被吸收吸收, 分布分布, 代谢代谢, 最终排除体外最终排除体外. 药物在血液中的浓度药物在血液中的浓度, 即单即单位体积血液中的药物含量位体积血液中的药物含量, 称血药浓度称血药浓度. 在最简单的一室模型在最简单的一室模型中中, 将整个机体看作一个房室将整个机体看作一个房室, 称中心室称中心室, 室

33、内的血药浓度是室内的血药浓度是均匀的均匀的. 快速静脉注射后快速静脉注射后, 浓度立即上升浓度立即上升; 然后逐渐下降然后逐渐下降. 当浓当浓度太低时度太低时, 达不到预期的治疗效果达不到预期的治疗效果; 血药浓度太高血药浓度太高, 又可能导又可能导致药物中毒或副作用太强致药物中毒或副作用太强. 临床上临床上, 每种药物有一个最小有效每种药物有一个最小有效浓度浓度 c1 和一个最大治疗浓度和一个最大治疗浓度 c2. 设计给药方案时设计给药方案时, 要使血药要使血药浓度保持在浓度保持在 c1-c2 之间之间. 设本题所研究药物的最小有效浓度设本题所研究药物的最小有效浓度c1=10, 最大治疗浓度

34、最大治疗浓度 c2=25第34页/共44页第三十五页，共44页。例例5 给药方案给药方案(fng n) 显然显然, 要设计给药方案要设计给药方案, 必须知道给药后血药浓度随时必须知道给药后血药浓度随时间变化的规律间变化的规律. 为此为此(wi c), 从实验和理论两方面着手从实验和理论两方面着手.在实验方面在实验方面, 对某人用快速静脉注射方式一次注入该药对某人用快速静脉注射方式一次注入该药物物300mg后后, 在一定时刻在一定时刻 t (小时小时)采集血样采集血样, 测得血药浓度测得血药浓度c 如表如表: 血药浓度血药浓度c(t) 的测试数据的测试数据t0.250.511.523468c 1

35、9.21 18.1515.3614.1012.899.327.455.243.01第35页/共44页第三十六页，共44页。例例5 给药方案给药方案(fng n)近似直线关系近似直线关系, 即即 c(t) 有按负指数规律减少的趋势有按负指数规律减少的趋势. logtc第36页/共44页第三十七页，共44页。例例5 给药方案给药方案(fng n)1. 确定确定(qudng)血药浓度的血药浓度的变化规律变化规律假设假设(jish):a) 药物向体外排除的速率与中心室的血药浓度成正比药物向体外排除的速率与中心室的血药浓度成正比, 比例系数为比例系数为 k(0), 称排除速率称排除速率.b) 中心室血液容积为常数中心室血液容积为常数 V, t=0 瞬时注入药物的瞬时注入药物的剂量为剂量为 d, 血药浓度立即为血药浓度立即为.dV由假设由假设 a), 中心室的血药浓度中心室的血药浓度 c(t)应满足微分方程应满足微分方程dckcdt 由假设由假设 b), 方程的初始条件为方程的初始条件为: 0.dcV 求解得求解得: .ktdc teV 即血药浓度即血药浓度c(t)按指数规律下降按指数规律下降.第37页/共44页第三十八页，共4