三、插值、曲线拟合插值就是已知一组离散的数据点集,在集合内部某两个点之间预测函数值的方法。插值法是实用的数值方法,是函数逼近的重要方法。在生产和科学实验中,自变量x与因变量y的函数y=f(x)的关系式有时不能直接写出表达式,而只能得到函数在若干个点的函数值或导数值。当要求知道观测点之外的函数值时,需要估计函数值在该点的值。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三次插值。四种插值方法比较函数格式功能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双三次插值。函数格式功能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。注意:(1)只对已知数据点集内部的点进行的插值运算称为内插,可比较准确的估测插值点上的函数值。(2)当插值点落在已知数据集的外部时的插值称为外插,要估计外插函数值很难。MATLAB对已知数据集外部点上函数值的预测都返回NaN,但可通过为interp1函数添加‘extrap’参数指明也用于外插。MATLAB的外插结果偏差较大。例1对在[-1,1]上,用n=20的等距分点进行线性插值,绘制f(x)及插值函数的图形.2119fxx解在命令窗口输入: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。包角与包角系数包角(°)90100110120125130135140包角系数0.690.740.780.820.840.860.880.89包角(°)145150155160165170175180包角系数0.910.920.930.950.960.980.991a1=[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例3.已知实验数据如表。x00.250.500.751.00y0.91620.81090.69310.55960.4055计算插值点x0=0.6处的函数值y0。x=[00.250.500.751.00];y=[0.91620.81090.69310.55960.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.641831200000000例4对在[-5,5]上,用n=11个等距分点作分段线性插值和三次样条插值,用m=21个插值点作图,比较结果.211yx解在命令窗口输入: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)')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例4对在[-5,5]上,用n=11个等距分点作分段线性插值和三次样条插值,用m=21个插值点作图,比较结果.211yxxyy1y2解在命令窗口输入:例5在一天24h内,从零点开始每间隔2h测得的环境温度为12,9,9,10,18,24,28,27,25,20,18,15,13C(单位:)推测在每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-')例6在飞机的机翼加工时,由于机翼尺寸很大,通常在图纸上只能标出部分关键点的数据.某型号飞机的机翼上缘轮廓线的部分数据如下:x04.749.051938577695114133y05.238.111.9716.1517.116.3414.6312.166.69x152171190y7.033.990例6在飞机的机翼加工时,由于机翼尺寸很大,通常在图纸上只能标出部分关键点的数据.某型号飞机的机翼上缘轮廓线的部分数据如下: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)例7天文学家在1914年8月份的7次观测中,测得地球与金星之间距离(单位:m),并取其常用对数值与日期的一组历史数据如下所示,试推断何时金星与地球的距离(单位:m)的对数值为9.9352.日期18202224262830距离对数9.96189.95449.94689.93919.93129.92329.9150解由于对数值9.9352位于24和26两天所对应的对数值之间,所以对上述数据用三次样条插值加细为步长为1的数据:解由于对数值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.9150x=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','Curvebyinterpft')Lagrange插值方法介绍对给定的n个插值点及对应的函数值,利用n次Lagrange插值多项式,则对插值区间内任意x的函数值y可通过下式求的:MATLAB实现12,,,nxxx12,,,nyyy11()()nnjkkjkjjkxxyxyxxfunctionyi=lagrange(x,y,xi)yi=zeros(size(xi));np=length(y);fori=1:npz=ones(size(xi));forj=1:npifi~=jz=z.*(xi-x(j))/(x(i)-x(j));endendyi=yi+z*y(i);end11()()nnjkkjkjjkxxyxyxx例:给出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曲线拟合据人口统计年鉴,知我国从1949年至1994年人口数据资料如下:(人口数单位为:百万)(1)在直角坐标系上作出人口数的图象。(2)建立人口数与年份的函数关系,并估算1999年的人口数。年份19491954195919641969人口数541.67602.66672.09704.99806.71年份19741979198419891994人口数908.59975.421034.751106.761176.74如何确定a,b?baxy线性模型1曲线拟合问题的提法:已知一组(二维)数据,即平面上的n个点),(iiyx,ixni,,,2,1L互不相同,寻求一个函数(曲线))(xfy,使)(xf在某种准则下与所有数据点最为接近,即曲线拟合得最好,如图:xy0++++++++i),(iiyx)(xfyniiiniiyxf1212])([确定f(x)使得达到最小最小二乘准则2.用什么样的曲线拟合已知数据?常用的曲线函数系类型:1.画图观察;2.理论分析xaeay21指数曲线:21axay双曲线(一支):11mmmaxaxayL多项式:21axay直线:3拟合函数组中系数的确定达到最小。niiiniiyaxaaaJ12211221])[),(为此,只需利用极值的必要条件)2