matlab插值曲线拟合省公共课一等奖全国赛课获奖课件.pptx

三、插三、插值、曲、曲线拟合合第1页插值就是已知一组离散数据点集，在集合内部某两个点之间预测函数值方法。插值法是实用数值方法，是函数迫近主要方法。在生产和科学试验中，自变量x与因变量y函数y=f(x)关系式有时不能直接写出表示式，而只能得到函数在若干个点函数值或导数值。当要求知道观察点之外函数值时，需要预计函数值在该点值。第2页Matlab惯用数据插值函数及功效函数格式功效Interp1yi=interp1(x,y,xi)一维插值。对一组点(x,y)进行插值，计算插值点xi函数值。yi=interp1(y,xi)此格式默认x=1:n，其中n是向量y长度，y为矩阵时，n=size(y,1)yi=interp1(x,y,xi,method)Method为指定插值算法。Nearest最近邻点插值，linear线性插值(默认方式)，spline三次样条函数插值，cubic三次插值。第3页四种插四种插值方法比方法比较第4页函数格式功效Interp2zi=interp2(x,y,z,xi,yi)二维插值。Z为由已知点值组成矩阵，参量x与y是与z同维已知点矩阵，且必须是单调。xi与yi为需要插值点。若xi与yi中有在x与y范围之外点，则对应地返回NaN。zi=interp2(z,xi,yi)此格式默认x=1:n、y=1:m，其中m,n=size(z)。再按zi=interp2(x,y,z,xi,yi)情形进行计算。zi=interp2(x,y,z,xi,yi,method)用指定算法method计算二维插值。Linear双线性插值(默认方式)，Nearest最邻近插值，spline三次样条插值，cubic双三次插值。第5页函数格式功效Interpfty=interpft(x,n)一维Fourier插值。适适用于周期函数生成数据插值。X为抽样序列，n为需要计算等距点数，y为n点插值计算结果。y=interpft(x,n,dim)沿着指定方向dim进行计算。splineyy=spline(x,y,xx)三次样条数据插值。用三次样条插值计算出由向量x与y确定一元函数在点xx处值yy。pp=spline(x,y)返回由向量x与y确定分段样条多项式系数矩阵pp。第6页注意：注意：（1）只）只对已知数据点集内部点已知数据点集内部点进行插行插值运算称运算称为内插，可内插，可比比较准确估准确估测插插值点上函数点上函数值。（2）当插）当插值点落在已知数据集外部点落在已知数据集外部时插插值称称为外插，要外插，要预计外插函数外插函数值极极难。MATLAB对已知数据集外部点上函数已知数据集外部点上函数值预测都返回都返回NaN，但，但可可经过为interp1函数添加函数添加extrap参数指明也用于外插。参数指明也用于外插。MATLAB外插外插结果偏差果偏差较大。大。第7页例例1 对 在在-1,1上上,用用n=20等距分点等距分点进行线性插值进行线性插值,绘制绘制 f(x)及插值函数图形及插值函数图形.第8页解解 在命令窗口输入在命令窗口输入: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,*)第9页例例2.在普通在普通V带设计中，中，带轮包角包角与包角系数与包角系数ka之之间关系如表所关系如表所表示。求表示。求=133.5时包角系数包角系数ka。包角与包角系数包角()90100110120125130135140包角系数0.690.740.780.820.840.860.880.89包角()145150155160165170175180包角系数0.910.920.930.950.960.980.991第10页a1=90,100,110,120,125,130,135,140,145,150,155,160,165,170,175,180;a2=0.69,0.74,0.78,0.82,0.84,0.86,0.88,0.89,0.91,0.92,0.93,0.95,0.96,0.98,0.99,1;ka=interp1(a1,a2,133.5)ka=0.8740第11页例例3.已知已知试验数据如表数据如表。x00.250.500.751.00y0.91620.81090.69310.55960.4055计算插算插值点点x0=0.6处函数函数值y0。第12页 x=0 0.25 0.50 0.75 1.00;y=0.9162 0.8109 0.6931 0.5596 0.4055;x0=0.6;y01=interp1(x,y,x0);y02=interp1(x,y,x0,nearest);y03=interp1(x,y,x0,pchip);y04=interp1(x,y,x0,spline);y01,y02,y03,y04y01=0.639700000000000y02=0.693100000000000y03=0.641818996851421y04=0.64183100000第13页例例 4 对 在在-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)第14页 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.0.0.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例例 4 对 在在-5,5上上,用用n=11个等距分点作分段线个等距分点作分段线性插值和三次样条插值性插值和三次样条插值,用用m=21个插值点作图个插值点作图,比较结果比较结果.xyy1y2第15页解解 在命令窗口输入在命令窗口输入:例例 5 在一天在一天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-)第16页例例 6 在在飞机机翼加工机机翼加工时,因因为机翼尺寸很大机翼尺寸很大,通常在通常在图纸上只能上只能标出部分关出部分关键点数据点数据.某型号某型号飞机机翼上机机翼上缘轮廓廓线部分数据以下部分数据以下: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第17页例例 6 在在飞机机翼加工机机翼加工时,因因为机翼尺寸很大机翼尺寸很大,通常在通常在图纸上只能上只能标出部分关出部分关键点数据点数据.某型号某型号飞机机翼上机机翼上缘轮廓廓线部分数据以下部分数据以下: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)第18页例7 天文学家在198月份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数数据据:第19页解解 因为对数值因为对数值 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第20页 x=linspace(0,2*pi,11);y=sin(x).*exp(-x/5);xi=linspace(0,2*pi,21);yi=interpft(y,21);plot(x,y,o,xi,yi);legend(Original,Curve by interpft)第21页Lagrange插值n方法介绍 对给定n个插值点 及对应函数值 ，利用n次Lagrange插值多项式，则对插值区间内任意x函数值y可经过下式求：nMATLAB实现第22页nfunction yi=lagrange(x,y,xi)nyi=zeros(size(xi);nnp=length(y);nfor i=1:npn z=ones(size(xi);n for j=1:npn if i=jn z=z.*(xi-x(j)/(x(i)-x(j);n endn endn yi=yi+z*y(i);nend第23页例：给出f(x)=ln(x)数值表，用Lagrange计算ln(0.54)近似值。x=0.4:0.1:0.8;y=log(x);x0=0.54;y0=lagrange(x,y,x0);ys=log(x0);y0,ysy0=-0.616142610505419ys=-0.616186139423817第24页曲曲线拟合合第25页据人口据人口统计年年鉴，知我国从，知我国从1949年至年至1994年人口数据年人口数据资料以下：料以下：(人口数人口数单位位为：百万：百万)（1）在直角坐）在直角坐标系上作出人口数系上作出人口数图象。象。（2）建立人口数与年份函数关系，并估算）建立人口数与年份函数关系，并估算1999年人年人口数。口数。年份年份19491954 1959 1964 1969人口数人口数 541.67602.66 672.09 704.99 806.71 年份年份 1974 1979 1984 1989 1994人口数人口数 908.59 975.42 1034.751106.761176.74 第26页怎样确定怎样确定a,b?线性模型性模型第27页1 曲曲线拟合合问题提法提法:已知一组（二维）数据，即平面上已知一组（二维）数据，即平面上n个点个点),(iiyx，ixni,2,1L=互不相同，寻求一个函数（曲线）互不相同，寻求一个函数（曲线）)(xfy=，使使)(xf在某种准则下与全部数据点最为靠近，即曲线拟合在某种准则下与全部数据点最为靠近，即曲线拟合得最好，如图得最好，如图:xy0+确定确定f(x)使得使得 到达最小到达最小 最小二乘准最小二乘准则 第28页.用什么用什么样曲曲线拟合已知数据合已知数据?惯用曲用曲线函数系函数系类型：型：画图观察；画图观察；理论分析理论分析指数曲线：指数曲线：双曲线（一支双曲线（一支):):多项式：多项式：直线：直线：第29页 拟合函数组中系数确定第30页第31页na=polyfit（xdata，ydata，n）n其中其中n表示多表示多项式最高式最高阶数数 nxdata，ydata 为要要拟合数据，它是用向量方合数据，它是用向量方式式输入。入。n输出参数出参数a为拟合多合多项式式 y=a1xn+anx+an+1系数系数a=a1,an,an+1。n多多项式在式在x处值y可用下面程序可用下面程序计算。算。n y=polyval(a,x)第32页例例a 已知观察数据点如表所表已知观察数据点如表所表示示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程序以下程序以下:第33页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)第34页例例b 用切削机床进行金属品加工时用切削机床进行金属品加工时,为了适当地调整为了适当地调整机床机床,需要测定刀具磨损速度需要测定刀具磨损速度.在一定时间测量刀具厚在一定时间测量刀具厚度度,得数据如表所表示得数据如表所表示:切削时间切削时间 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第35页解解:在命令窗口输入在命令窗口输入: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,*)第36页解解:在命令窗口输入在命令窗口输入: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第37页例例b 用切削机床进行金属品加工时用切削机床进行金属品加工时,为了适当地调整为了适当地调整机床机床,需要测定刀具磨损速度需要测定刀具磨损速度.在一定时间测量刀具厚在一定时间测量刀具厚度度,得数据如表所表示得数据如表所表示:切削时间切削时间 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第38页例例c 一个一个15.4cm30.48cm混凝土柱在加压试验中应力混凝土柱在加压试验中应力-应变关系测试点数据如表所表示应变关系测试点数据如表所表示1.552.472.933.03已知应力已知应力-应变关系能够用一条指数曲线来描述应变关系能够用一条指数曲线来描述,即假设即假设式中式中,表示应力表示应力,单位是单位是 N/m2;表示应变表示应变.2.89第39页已知应力已知应力-应变关系能够用一条指数曲线来描述应变关系能够用一条指数曲线来描述,即假设即假设式中式中,表示应力表示应力,单位是单位是 N/m2;表示应变表示应变.解解 选取选取指数函数指数函数作拟合时作拟合时,在拟合前需作变量代换在拟合前需作变量代换,化为化为 k1,k2 线性函数线性函数.于是于是,令令即即第40页在命令窗口输入在命令窗口输入:x=500*1.0e-6 1000*1.0e-6 1500*1.0e-6*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-)第41页已知应力已知应力-应变关系能够用一条指数曲线来描述应变关系能够用一条指数曲线来描述,即假设即假设式中式中,表示应力表示应力,单位是单位是 N/m2;表示应变表示应变.拟合曲线为拟合曲线为:令令则则求得求得于是于是第42页在实际应用中常见拟合曲线有在实际应用中常见拟合曲线有:直线直线多项式多项式普通普通 n=2,3,不宜过高不宜过高.双曲线双曲线(一支一支)指数曲线指数曲线第43页例例:海底光海底光缆线长度度预测模型模型某某一一通通信信企企业在在一一次次施施工工 中中,需需要要在在水水面面宽为2 20 0m m河河沟沟底底沿沿线走走向向铺设一一条条沟沟底底光光缆.在在铺设光光缆之之前前需需要要对沟沟底底地地形形做做初初B2468101214161820986420ADC探探测到一到一组等分点位置深度数据以下表所表示等分点位置深度数据以下表所表示.步探步探测,从而从而预计所需光所需光缆长度度,为工程工程预算提算提供依据供依据.基本情况如基本情况如图所表示所表示.第44页10.9310.809.818.867.957.959.1510.2211.2912.6113.32201918171615141312111013.2812.2611.1810.139.058.027.967.968.969.01深度深度(m)9876543210分点分点21个等分点个等分点处深度深度(1)(1)预测经过这条河沟所需光条河沟所需光缆长度近似度近似值.(2)(2)作出作出铺设沟底光沟底光缆曲曲线图.第45页解：解：用用12次多次多项式函数式函数拟合光合光缆走走势曲曲线图以下以下仿真仿真结果表果表明明,拟合曲合曲线能能较准确地准确地反反应光光缆走走势图.The length of the label is L=26.3809(m)假假设所所铺设光光缆足足够柔柔软,在在铺设过程中光程中光缆触地走触地走势光滑光滑,紧贴地面地面,而且忽略水流而且忽略水流对光光缆冲冲击.第46页format longt=linspace(0,20,21);x=linspace(0,20,100);P=9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93;a=polyfit(t,P,12);yy=polyval(a,x);plot(x,yy,r*-,t,P,b+-);L=0;for i=2:100 L=L+sqrt(x(i)-x(i-1)2+(yy(i)-yy(i-1)2);endL第47页2.非线性曲线拟合非线性曲线拟合:lsqcurvefit.功效功效:x=lsqcurvefit(fun,x0,xdata,ydata)x,resnorm=lsqcurvefit(fun,x0,xdata,ydata)依据给定数据依据给定数据 xdata,ydata(对应点横对应点横,纵纵坐标坐标),按函数文件按函数文件 fun 给定函数给定函数,以以x0为初为初值作值作最小二乘拟合最小二乘拟合,返回函数返回函数 fun中系数向中系数向量量x和残差平方和和残差平方和resnorm.第48页例例d 已知观察数据点如表所表已知观察数据点如表所表示示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第49页例例d 已知观察数据点如表所表已知观察数据点如表所表示示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)第50页编写下面程序调用拟合函数编写下面程序调用拟合函数.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第51页例例d 已知观察数据点如表所表已知观察数据点如表所表示示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.第52页f(x)=3ex+4.03x2+0.94 x3.拟合函数为拟合函数为:第53页（酒后驾车）中给出某人在短时间内喝下两瓶啤酒后，间隔一定时间测量他血液中酒精含量y（毫克/百毫升），得到数据以下表所表示。时间(小时)0.25 0.50.75 11.522.533.544.55酒精含量306875828277686858515041时间(小时)678910111213141516酒精含量3835282518151210774酒精浓度与时间关系为：第54页function f=Example2_1(c,tdata)f=c(1)*(exp(-c(2)*tdata)-exp(-c(3)*tdata);cleartdata=0.25 0.5 0.75 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15 16;ydata=30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 8 15 12 10 7 7 4;c0=1 1 1;for i=1:50;c=lsqcurvefit(Example2_1,c0,tdata,ydata);c0=c;end得到最优解为：c=117.05,0.1930,1.9546第55页练习练习: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 与已知数据点与已知数据点在最小二乘意义上充分靠近在最小二乘意义上充分靠近.第56页