1、三、插值、曲线拟合,插值就是已知一组离散的数据点集，在集合内部某两个点之间预测函数值的方法。,插值法是实用的数值方法，是函数逼近的重要方法。在生产和科学实验中，自变量x与因变量y的函数y = f(x)的关系式有时不能直接写出表达式，而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时，需要估计函数值在该点的值。,Matlab常用数据插值函数及功能,四种插值方法比较,注意：（1）只对已知数据点集内部的点进行的插值运算称为内插，可比较准确的估测插值点上的函数值。（2）当插值点落在已知数据集的外部时的插值称为外插，要估计外插函数值很难。MATLAB对已知数据集外部点上函数值的预测
2、都返回NaN，但可通过为interp1函数添加extrap参数指明也用于外插。 MATLAB的外插结果偏差较大。,例1 对 在-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,*),例2. 在普通V带设计中，带轮的包角与包角系数ka之间的关系如表所示。求=133.5时的包角系数ka。,包角与包角系数,a1=90,100,110,120,125,130,135,140,145,150
3、,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,例3. 已知实验数据如表。,计算插值点x0=0.6处的函数值y0。, 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
4、=interp1(x,y,x0,pchip); y04=interp1(x,y,x0,spline); y01,y02,y03,y04,y01 =0.639700000000000 y02 =0.693100000000000 y03 =0.641818996851421 y04 =0.641831200000000,例 4 对 在-5, 5上, 用n=11个等距分点作分段线性插值和三次样条插值, 用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./(
5、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.00000.5000 0.8000 0.7500 0.82051.0000 0.5000 0.5000 0.50001.5000 0.3077 0.3500 0.29732.0000 0.2000 0.2000 0.20002.5000 0.1379 0.150
6、0 0.14013.0000 0.1000 0.1000 0.10003.5000 0.0755 0.0794 0.07454.0000 0.0588 0.0588 0.05884.5000 0.0471 0.0486 0.04845.0000 0.0385 0.0385 0.0385,例 4 对 在-5, 5上, 用n=11个等距分点作分段线性插值和三次样条插值, 用m=21个插值点作图,比较结果.,x,y,y1,y2,解 在命令窗口输入:,例 5 在一天24h内, 从零点开始每间隔2h测得的环境温度为,12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15
7、, 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,spline) plot(t,T,*,ti,T1i,r-,ti,T2i,g-),例 6 在飞机的机翼加工时, 由于机翼尺寸很大, 通常在图纸上只能标出部分关键点的数据. 某型号飞机的机翼上缘轮廓线的部分数据如下:,x 0 4.74 9.05 19 38 57 7
8、6 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,例 6 在飞机的机翼加工时, 由于机翼尺寸很大, 通常在图纸上只能标出部分关键点的数据. 某型号飞机的机翼上缘轮廓线的部分数据如下:,x=0 4.74 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
9、(x,y,xi,spline) plot(xi,yi),例7 天文学家在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 两
10、天所对应的对数值之间, 所以对上述数据用三次样条插值加细为步长为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.0000 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
11、 9.9352 9.9312 9.9272 9.9232 9.9191 9.9150, 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),Lagrange插值,方法介绍对给定的n个插值点 及对应的函数值 ，利用n次Lagrange插值多项式，则对插值区间内任意x的函数值y可通过下式求的：MATLAB实现,function yi=lagrange(x,y,xi) yi
12、=zeros(size(xi); np=length(y); for i=1:npz=ones(size(xi);for j=1:npif i=jz=z.*(xi-x(j)/(x(i)-x(j);endendyi=yi+z*y(i); end,例：给出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.616142610505419 ys =-0.616186139423817,曲线拟合,据人口统计年鉴，知我国
13、从1949年至1994年人口数据资料如下： (人口数单位为：百万),（1）在直角坐标系上作出人口数的图象。 （2）建立人口数与年份的函数关系，并估算1999年的人口数。,如何确定a,b?,线性模型,1 曲线拟合问题的提法:,x,y,0,+,+,+,+,+,+,+,+,确定f(x)使得,达到最小,最小二乘准则,. 用什么样的曲线拟合已知数据?,常用的曲线函数系类型：,画图观察； 理论分析,指数曲线：,双曲线（一支):,多项式：,直线：, 拟合函数组中系数的确定,a = polyfit（xdata，ydata，n） 其中n表示多项式的最高阶数 xdata，ydata 为要拟合的数据，它是用向量的方
14、式输入。 输出参数a为拟合多项式 y = a1xn + + anx + an+1的系数a = a1, , an, an+1。 多项式在x处的值y可用下面程序计算。y = polyval (a, x),例a 已知观测数据点如表所示,分别用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 16) p3=polyfit(x,y,3) p6=polyfit(x,y,6),编写Matlab程序如下:
15、,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=polyfit(x,y,6),例b 用切削机床进行金属品加工时, 为了适当地调整机床, 需要测定刀具的磨损速度.
16、 在一定的时间测量刀具的厚度, 得数据如表所示:,解: 在命令窗口输入:,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.1 25.7 25.3 24.8 24.0 plot(t,y,*),a =-0.3012 29.3804,hold on,plot(t,
17、y1), hold off,a=polyfit(t,y,1),y1=-0.3012*t+29.3804,例b 用切削机床进行金属品加工时, 为了适当地调整机床, 需要测定刀具的磨损速度. 在一定的时间测量刀具的厚度, 得数据如表所示:,切削时间 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.3,12,26.1,13,25.7,14,25.3,15,24.8,16,24.0,刀具厚度 y/cm,拟合曲线为:,y=-0.3012t+
18、29.3804,例c 一个15.4cm30.48cm的混凝土柱在加压实验中的应力-应变关系测试点的数据如表所示,1.55,2.47,2. 93,3. 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
19、*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-),已知应力-应变关系可以用一条指数曲线来描述, 即假设,式中, 表示应力, 单位是 N/m2; 表示应变.,拟合曲线为:,令,则,求得,于是,在实际应用中常见的拟合曲线有:,直线,
20、多项式,一般 n=2, 3, 不宜过高.,双曲线(一支),指数曲线,例: 海底光缆线长度预测模型,某一通信公司在一次施工中,需要在水面宽为20m的河沟底沿线走向铺设一条沟底光缆.在铺设光缆之前需要对沟底的地形做初,探测到一组等分点位置的深度数据如下表所示.,步探测,从而估计所需光缆的长度,为工程预算提供依据.基本情况如图所示.,(1) 预测通过这条河沟所需光缆长度的近似值.,(2) 作出铺设沟底光缆的曲线图.,解： 用12次多项式函数拟合光缆走势的曲线图如下,仿真结果表明,拟合曲线能较准确地反映光缆的走势图.,The length of the label is L= 26.3809 (m),
21、假设所铺设的光缆足够柔软,在铺设过程中光缆触地走势光滑,紧贴地面,并且忽略水流对光缆的冲击.,format long t=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:100L=L+
22、sqrt(x(i)-x(i-1)2+(yy(i)-yy(i-1)2); end L,2. 非线性曲线拟合: lsqcurvefit.,功能:,x=lsqcurvefit(fun, x0, xdata, ydata),x, resnorm=lsqcurvefit(fun, x0, xdata, ydata),根据给定的数据 xdata, ydata (对应点的横, 纵坐标), 按函数文件 fun 给定的函数, 以x0为初值作最小二乘拟合, 返回函数 fun中的系数向量x和残差的平方和resnorm.,例d 已知观测数据点如表所示,求三个参数 a, b, c的值, 使得曲线 f(x)=aex+bx
23、2+cx3 与已知数据点在最小二乘意义上充分接近.,首先编写存储拟合函数的函数文件.,function f=nihehanshu(x,xdata) f=x(1)*exp(xdata)+x(2)*xdata.2+x(3)*xdata.3,保存为文件 nihehanshu.m,例d 已知观测数据点如表所示,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 与已知数据点在最小二乘意义
24、上充分接近.,编写下面的程序调用拟合函数.,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,ydat
25、a),程序运行后显示,x =3.0022 4.0304 0.9404,resnorm =0.0912,例d 已知观测数据点如表所示,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 与已知数据点在最小二乘意义上充分接近.,说明: 最小二乘意义上的最佳拟合函数为,f(x)= 3ex+ 4.03x2 + 0.94 x3.,此时的残差是: 0.0912.,f(x)= 3ex+ 4.0
26、3x2 + 0.94 x3.,拟合函数为:,（酒后驾车）中给出某人在短时间内喝下两瓶啤酒后，间隔一定的时间测量他的血液中酒精含量y（毫克/百毫升），得到数据如下表所示。,酒精浓度与时间的关系为：,function f=Example2_1(c,tdata) f=c(1)*(exp(-c(2)*tdata)-exp(-c(3)*tdata);,clear tdata=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,练习:,1. 已知观测数据点如表所示,求用三次多项式进行拟合的曲线方程.,2. 已知观测数据点如表所示,求a, b, c的值, 使得曲线 f(x)=aex+bsin x+c lnx 与已知数据点在最小二乘意义上充分接近.,
