数学建模中Matlab数据拟合应用

数学建模中Matlab数据拟合应用第一页，共44页。例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.2]plot(x,y,'k.','markersize',25)axis([01.3-216])p3=polyfit(x,y,3)p6=polyfit(x,y,6)编写Matlab程序如下:第二页，共44页。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)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.2]plot(x,y,'k.','markersize',25)axis([01.3-216])p3=polyfit(x,y,3)p6=polyfit(x,y,6)第三页，共44页。例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第四页，共44页。解:描出散点图,在命令窗口输入:t=[0:1:16]y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0]plot(t,y,'*')第五页，共44页。解:描出散点图,在命令窗口输入:t=[0:1:16]y=[30.029.128.428.128.027.727.527.227.026.826.526.326.125.725.324.824.0]plot(t,y,'*')a=-0.301229.3804holdonplot(t,y1),holdoffa=polyfit(t,y,1)y1=-0.3012*t+29.3804第六页，共44页。例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第七页，共44页。例3一个15.4cm×30.48cm的混凝土柱在加压实验中的应力-应变关系测试点的数据如表所示1.552.472.933.03已知应力-应变关系可以用一条指数曲线来描述,即假设式中,表示应力,单位是N/m2;表示应变.2.89第八页，共44页。已知应力-应变关系可以用一条指数曲线来描述,即假设式中,表示应力,单位是N/m2;表示应变.解选取指数函数作拟合时,在拟合前需作变量代换,化为k1,k2的线性函数.于是,令即第九页，共44页。在命令窗口输入:x=[500*1.0e-61000*1.0e-61500*1.0e-62000*1.0e-62375*1.0e-6]y=[3.103*1.0e+32.465*1.0e+31.953*1.0e+31.517*1.0e+31.219*1.0e+3]z=log(y)a=polyfit(x,z,1)k1=exp(8.3009)w=[1.552.472.933.032.89]plot(x,w,'*')y1=exp(8.3009)*x.*exp(-494.5209*x)plot(x,w,'*',x,y1,'r-')第十页，共44页。已知应力-应变关系可以用一条指数曲线来描述,即假设式中,表示应力,单位是N/m2;表示应变.拟合曲线为:令则求得于是第十一页，共44页。在实际应用中常见的拟合曲线有:直线多项式一般n=2,3,不宜过高.双曲线(一支)指数曲线第十二页，共44页。2.非线性曲线拟合:lsqcurvefit.功能:x=lsqcurvefit(fun,x0,xdata,ydata)[x,resnorm]=lsqcurvefit(fun,x0,xdata,ydata)根据给定的数据xdata,ydata(对应点的横,纵坐标),按函数文件fun给定的函数,以x0为初值作最小二乘拟合,返回函数fun中的系数向量x和残差的平方和resnorm.第十三页，共44页。例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与已知数据点在最小二乘意义上充分接近.首先编写存储拟合函数的函数文件.functionf=nihehanshu(x,xdata)f=x(1)*exp(xdata)+x(2)*xdata.^2+x(3)*xdata.^3保存为文件nihehanshu.m第十四页，共44页。例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)第十五页，共44页。编写下面的程序调用拟合函数.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.9404resnorm=0.0912第十六页，共44页。例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.94x3.此时的残差是:0.0912.第十七页，共44页。f(x)=3ex+4.03x2+0.94x3.拟合函数为:第十八页，共44页。练习: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+bsinx+clnx与已知数据点在最小二乘意义上充分接近.第十九页，共44页。插值问题已知n+1个节点其中互不相同,不妨设求任一插值点处的插值节点可视为由产生,g

表达式复杂,甚至无表达式第二十页，共44页。1.分段线性插值xjxj-1xj+1x0xn实用插值方法机翼下轮廓线2.三次样条插值细木条:样条第二十一页，共44页。输入:节点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作插值计算第二十二页，共44页。例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,'*')第二十三页，共44页。例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'y2']plot(x,z,'r',x,y,'k:',x,y1,'b',x,y2,'g')gtext('Piece.-linear.'),gtext('Spline'),gtext('y=1/(1+x^2)')第二十四页，共44页。01.00001.00001.00000.50000.80000.75000.82051.00000.50000.50000.50001.50000.30770.35000.29732.00000.20000.20000.20002.50000.13790.15000.14013.00000.10000.10000.10003.50000.07550.07940.07454.00000.05880.05880.05884.50000.04710.04860.04845.00000.03850.03850.0385例6对在[-5,5]上,用n=11个等距分点作分段线性插值和三次样条插值,用m=21个插值点作图,比较结果.xyy1y2第二十五页，共44页。解在命令窗口输入:例7在一天24h内,从零点开始每间隔2h测得的环境温度为12,9,9,10,18,24,28,27,25,20,18,15,13(单位:)推测在每1s时的温度.并描绘温度曲线.t=0:2:24T=[129910182428272520181513]plot(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-')第二十六页，共44页。例8在飞机的机翼加工时,由于机翼尺寸很大,通常在图纸上只能标出部分关键点的数据.某型号飞机的机翼上缘轮廓线的部分数据如下:x04.749.051938577695114133y05.238.111.9716.1517.116.3414.6312.166.69x152171190y7.033.990第二十七页，共44页。例8在飞机的机翼加工时,由于机翼尺寸很大,通常在图纸上只能标出部分关键点的数据.某型号飞机的机翼上缘轮廓线的部分数据如下:x=[04.749.051938577695114133152171190]y=[05.238.111.9716.1517.116.3414.6312.169.697.033.990]xi=[0:0.001:190]yi=interp1(x,y,xi,'spline')plot(xi,yi)第二十八页，共44页。例9天文学家在1914年8月份的7次观测中,测得地球与金星之间距离(单位:m),并取其常用对数值与日期的一组历史数据如下所示,试推断何时金星与地球的距离(单位:m)的对数值为9.9352.日期18202224262830距离对数9.96189.95449.94689.93919.93129.92329.9150解由于对数值9.9352位于24和26两天所对应的对数值之间,所以对上述数据用三次样条插值加细为步长为1的数据:第二十九页，共44页。解由于对数值9.9352位于24和26两天所对应的对数值之间,所以对上述数据用三次样条插值加细为步长为1的数据:x=[18:2:30]y=[9.96189.95449.94689.93919.93129.92329.9150]xi=[18:1:30]yi=interp1(x,y,xi,'spline')A=[xi;yi]A=18.000019.000020.000021.000022.000023.000024.000025.000026.000027.000028.000029.000030.00009.96189.95819.95449.95069.94689.94309.93919.93529.93129.92729.92329.91919.9150第三十页，共44页。练习:1.设在区间[-2,2]上用10等分点作为节点,分别用三种插值方法:(1)计算并输出在该区间的20等分点的函数值.(2)输出这个函数及两个插值函数的图形.(3)对输出的数据和图形进行分析.第三十一页，共44页。1.设在区间[-2,2]上用10等分点作为节点,分别用三种插值方法:(1)计算并输出在该区间的20等分点的函数值.zi=0.01830.03870.07730.14110.23690.36850.52730.69800.85210.95991.00000.95990.85210.69800.52730.36850.23690.14110.07730.03870.0183第三十二页，共44页。1.设在区间[-2,2]上用10等分点作为节点,分别用两种插值方法:(2)输出这个函数及两个插值函数的图形.第三十三页，共44页。练习:2.已知某型号飞机的机翼断面下缘轮廓线上的部分数据如表所示:假设需要得到x坐标每改变0.1时的y坐标,分别用两种插值方法对机翼断面下缘轮廓线上的部分数据加细,并作出插值函数的图形.xy0031.251.772.092.1112.0121.8131.2141.0151.6第三十四页，共44页。例5给药方案。

一种新药用于临床之前,必须设计给药方案.在快速静脉注射的给药方式下,所谓给药方案是指,每次注射剂量多大,间隔时间多长.药物进入机体后随血液输送到全身,在这个过程中不断地被吸收,分布,代谢,最终排除体外.药物在血液中的浓度,即单位体积血液中的药物含量,称血药浓度.在最简单的一室模型中,将整个机体看作一个房室,称中心室,室内的血药浓度是均匀的.快速静脉注射后,浓度立即上升;然后逐渐下降.当浓度太低时,达不到预期的治疗效果;血药浓度太高,又可能导致药物中毒或副作用太强.临床上,每种药物有一个最小有效浓度c1和一个最大治疗浓度c2.设计给药方案时,要使血药浓度保持在c1-c2之间.设本题所研究药物的最小有效浓度c1=10,最大治疗浓度c2=25第三十五页，共44页。例5给药方案。

显然,要设计给药方案,必须知道给药后血药浓度随时间变化的规律.为此,从实验和理论两方面着手.在实验方面,对某人用快速静脉注射方式一次注入该药物300mg后,在一定时刻t(小时)采集血样,测得血药浓度c如表:

血药浓度c(t)的测试数据t0.250.511.523468c19.2118.1515.3614.1012.899.327.455.243.01第三十六页，共44页。