用MATLAB进行数据插值资料

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第十一讲 数据插值,刘北战,我们经常会遇到大量的数据需要处理，而处理数据的关键就在于这些算法，例如数据拟合、参数估计、插值等数据处理算法。此类问题在,MATLAB,中有很多现成的函数可以调用，熟悉,MATLAB,，这些方法都能游刃有余的用好。,一、,概述,数据拟合在很多赛题中有应用，与图形处理有关的问题很多与插值和拟合有关系，例如,98,年美国赛,A,题，生物组织切片的三维插值处理，,94,年,A,题逢山开路，山体海拔高度的插值计算，,2003,年吵的沸沸扬扬的“非典”问题也要用到数据拟合算法，观察数据的走向进行处理，,2005,年的雨量预报的评价的插值计算。,2001,年的公交车调度拟合问题，,2003,年的饮酒驾车拟合问题。,插值问题雨量预报的评价,预测点和实测点的图形,插值后的图形,拟合问题饮酒驾车,喝两瓶酒的拟合曲线,喝,1-5,瓶酒的拟合曲线,在实际中，常常要处理由实验或测量所得到的一些离散数据。插值与拟合方法就是要通过这些数据去确定某一类已知函数的参数或寻求某个近似函数，使所得到的近似函数与已知数据有较高的拟合精度。,如果要求这个近似函数（曲线或曲面）经过所已知的所有数据点，则称此类问题为,插值问题,。 （不需要函数表达式）,二、,基本概念,如果不要求近似函数通过所有数据点，而是要求它能较好地反映数据变化规律的近似函数的方法称为,数据拟合,。（必须有函数表达式）,近似函数不一定（曲线或曲面）通过所有的数据点。,1,、联系,都是根据实际中一组已知数据来构造一个能够反映数据变化规律的近似函数的方法。,2,、区别,插值问题,不一定得到近似函数的表达形式，仅通过插值方法找到未知点对应的值。,数据拟合,要求得到一个具体的近似函数的表达式。,三、插值与,拟合的区别和联系,插 值,当数据量不够，需要补充，且认定已有数据可信时,通常利用函数插值方法。,实际问题当中碰到的函数,f,(,x,),是各种各样的，有的表达式很复杂，有的甚至给不出数学的式子，只提供了一些离散数据，警如，某些点上的函数值和导数值。,拉格朗日插值,分段线性插值,三次样条插值,一 维 插 值,一、,插值的定义,二、插值的方法,三、用,Matlab,解插值问题,返回,返回,二维插值,一、,二维插值定义,二、网格节点插值法,三、用,Matlab,解插值问题,最邻近插值,分片线性插值,双线性插值,网格节点数据的插值,散点数据的插值,一维插值的定义,已知,n+1,个节点,其中,互不相同，不妨设,求任一插值点,处的插值,节点可视为由,产生,，,表达式复杂,，,或无封闭形式,，,或未知,.,。,构造一个,(,相对简单的,),函数,通过全部节点,即,再用,计算插值，即,返回,称为,拉格朗日插值基函数,。,已知函数,f,(,x,),在,n,+1,个点,x,0,x,1,x,n,处的函数值为,y,0,y,1,y,n,。求一,n,次多项式函数,P,n,(,x,),，使其满足：,P,n,(,x,i,)=,y,i,i,=0,1,n,.,解决此问题的拉格朗日插值多项式公式如下,其中,L,i,(,x,),为,n,次多项式：,拉格朗日,(Lagrange),插值,拉格朗日,(Lagrange),插值,特别地,:,两点一次,(,线性,),插值多项式,:,三点二次,(,抛物,),插值多项式,:,拉格朗日多项式插值的,这种振荡现象叫,Runge,现象,采用拉格朗日多项式插值：选取不同插值节点个数,n,+1,，其中,n,为插值多项式的次数，当,n,分别取,2,4,6,8,10,时，绘出插值结果图形,.,例,返回,To Matlab,lch(larg1),分段线性插值,n,越大，误差越小,.,x,j,x,j-1,x,j+1,x,0,x,n,x,o,y,To MATLAB,xch11,，,xch12,，,xch13,，,xch14,返回,例,用分段线性插值法求插值,并观察插值误差,.,1.,在,-6,6,中平均选取,5,个点作插值,(xch11),4.,在,-6,6,中平均选取,41,个点作插值,(xch14),2.,在,-6,6,中平均选取,11,个点作插值,(xch12),3.,在,-6,6,中平均选取,21,个点作插值,(xch13),比分段线性插值更光滑。,x,y,x,i-1,x,i,a,b,在数学上，光滑程度的定量描述是：函数,(,曲线,),的,k,阶导数存在且连续，则称,该曲线具有,k,阶光滑性,。,光滑性的阶次越高，则越光滑。是否存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例子。,三次样条插值,三次样条插值,g,(,x,),为被插值函数,。,例,用三次样条插值选取,11,个基点计算插值,(ych),返回,To MATLAB,ych,用,MATLAB,作插值计算,一维插值函数：,yi=interp1(x,，,y,，,xi,，,method),插值方法,被插值点,插值节点,xi,处的插值结果,nearest,：最邻近插值,linear,：,线性插值；,spline,： 三次样条插值；,cubic,： 立方插值。,缺省时： 分段线性插值。,注意：所有的插值方法都要求,x,是单调的，并且,xi,不能够超过,x,的范围。,（,1,）,nearest,方法速度最快，占用内存最小，但一般来说误差最大，插值结果最不光滑；,(2) liner,分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函数预测，,MATLAB,中,interp1,的默认方法,（,3,）,spline,三次样条插值是所有插值方法中运行耗时最长的，其插值函数以及插值函数的一阶、二阶导函数都连续，因此是最光滑的插值方法，占用内存上比,cubic,方法小，但当已知数据点不均匀分布时可能出现异常结果。,（,4,）,cubic,三次多项式插值法中插值函数及其一阶导数都是连续的，因此其插值结果也比较光滑，运算速度比,spline,方法略快，但占用内存最多。在实际的使用中，应根据实际需求和运算条件选择合适的算法。,例 用其他一维插值方法对以下,7,个离散数据点,(1,3.5),、,(2,2.1),、,(3,1.3),、,(4.0.8),、,(5,2.9),、,(6,4.2),、,(7,5.7),进行一维插值方法。,解：,在,MATLAB,命令窗口中输入以下命令：, x=1 2 3 4 5 6 7;, y=3.5 2.1 1.3 0.8 2.9 4.2 5.7;, xx=1:0.5:7;, y1=interp1(x,y,xx,nearest);, y2=interp1(x,y,xx,spline);, y3=interp1(x,y,xx,cubic);, plot(x,y,o,xx,y1,-,xx,y2,-.,xx,y3,:),例：在,1-12,的,11,小时内，每隔,1,小时测量一次温度，测得的温度依次为：,5,，,8,，,9,，,15,，,25,，,29,，,31,，,30,，,22,，,25,，,27,，,24,。试估计每隔,1/10,小时的温度值。,To MATLAB,(temp),hours=1:12;,temps=5 8 9 15 25 29 31 30 22 25 27 24;,h=1:0.1:12;,t=interp1(hours,temps,h,spline); (,直接输出数据将是很多的,),plot(hours,temps,+),hold on,plot(h,t,hours,temps,r:) %,作图,xlabel(Hour),ylabel(Degrees Celsius),x,y,机翼下轮廓线,例 已知飞机下轮廓线上数据如下，求,x,每改变,0.1,时的,y,值。,To MATLAB(plane),返回,二维插值的定义,x,y,O,第一种（网格节点）：,已知,m,n,个节点,其中,互不相同，不妨设,构造一个二元函数,通过全部已知节点,即,再用,计算插值，即,第二种（散乱节点）：,y,x,0,已知,n,个节点,其中,互不相同，,构造一个二元函数,通过全部已知节点,即,再用,计算插值，即,返回,注意：,最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。,最邻近插值,x,y,(,x,1,y,1,),(,x,1,y,2,),(,x,2,y,1,),(,x,2,y,2,),O,二维或高维情形的最邻近插值，与被插值点最邻近的,节点的函数值即为所求。,返回,将四个插值点（矩形的四个顶点）处的函数值依次简记为：,分片线性插值,x,y,(,x,i,y,j,),(,x,i,y,j,+1,),(,x,i,+1,y,j,),(,x,i,+1,y,j,+1,),O,f,(,x,i,y,j,)=,f,1,，,f,(,x,i,+1,y,j,)=,f,2,，,f,(,x,i,+1,y,j,+1,)=,f,3,，,f,(,x,i,y,j,+1,)=,f,4,插值函数为：,第二片,(,上三角形区域,),：,(,x,y,),满足,插值函数为：,注意,：,(,x,y,),当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；,分两片的函数表达式如下：,第一片,(,下三角形区域,),：,(,x,y,),满足,返回,双线性插值是一片一片的空间二次曲面构成。,双线性插值函数的形式如下：,其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。,双线性插值,x,y,(,x,1,y,1,),(,x,1,y,2,),(,x,2,y,1,),(,x,2,y,2,),O,返回,要求,x0,y0,单调；,x,，,y,可取,为矩阵，或,x,取行向量，,y,取为列向量，,x,y,的值分别不能超出,x0,y0,的范围。,z=interp2(x0,y0,z0,x,y,method),被插值点,插值方法,用,MATLAB,作网格节点数据的插值,插值节点,被插值点的函数值,nearest,最邻近插值,linear,双线性插值,cubic,双三次插值,缺省时,双线性插值,例：测得平板表面,3*5,网格点处的温度分别为：,82 81 80 82 84 79 63 61 65 81 84 84 82 85 86,试作出平板表面的温度分布曲面,z=f(x,y),的图形。,输入以下命令：,x=1:5;,y=1:3;,temps=82 81 80 82 84;79 63 61 65 81;84 84 82 85 86;,mesh(x,y,temps),1.,先在三维坐标画出原始数据，画出粗糙的温度分布曲图,.,2,以平滑数据,在,x,、,y,方向上每隔,0.2,个单位的地方进行插值,.,再输入以下命令,:,xi=1:0.2:5;,yi=1:0.2:3;,zi=interp2(x,y,temps,xi,yi,cubic);,mesh(xi,yi,zi),画出插值后的温度分布曲面图,.,To MATLAB,(wendu),例, x,y=meshgrid(-3:.6:3,-2:.4:2);, z=(x.2-2*x).*exp,(-x.2-y.2-x.*y);, surf(x,y,z),选较密的插值点，用默认的线性插值算法进行插值, x1,y1=meshgrid(-3:.2:3,-2:.2:2);, z0=interp2(x,y,z,x1,y1);, surf(x1,y1,z0),立方和样条插值：, z1=interp2(x,y,z,x1,y1,cubic);, z2=interp2(x,y,z,x1,y1,spline);, surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5), figure;surf(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5),算法误差比较, z=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);, surf(x1,y1,abs(z-z1), figure;surf(x1,y1,abs(z-z2), figure;surf(x1,y1,abs(z-z0),山区地形地貌图,已知某处山区地形选点测量坐标数据为：,x=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5,y=0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6,海拔高度数据为：,z=89 90 87 85 92 91 96 93 90 87 82,92 96 98 99 95 91 89 86 84 82 84,96 98 95 92 90 88 85 84 83 81 85,80 81 82 89 95 96 93 92 89 86 86,82 85 87 98 99 96 97 88 85 82 83,82 85 89 94 95 93 92 91 86 84 88,88 92 93 94 95 89 87 86 83 81 92,92 96 97 98 96 93 95 84 82 81 84,85 85 81 82 80 80 81 85 90 93 95,84 86 81 98 99 98 97 96 95 84 87,80 81 85 82 83 84 87 90 95 86 88,80 82 81 84 85 86 83 82 81 80 82,87 88 89 98 99 97 96 98 94 92 87,山区地形地貌图 程序,原始地貌图程序：,x=0:.5:5;,y=0:.5:6;,xx,yy=meshgrid(x,y);,z=89 90 87 85 92 91 96 93 90 87 82,92 96 98 99 95 91 89 86 84 82 84,96 98 95 92 90 88 85 84 83 81 85,80 81 82 89 95 96 93 92 89 86 86,82 85 87 98 99 96 97 88 85 82 83,82 85 89 94 95 93 92 91 86 84 88,88 92 93 94 95 89 87 86 83 81 92,92 96 97 98 96 93 95 84 82 81 84,85 85 81 82 80 80 81 85 90 93 95,84 86 81 98 99 98 97 96 95 84 87,80 81 85 82 83 84 87 90 95 86 88,80 82 81 84 85 86 83 82 81 80 82,87 88 89 98 99 97 96 98 94 92 87;,mesh(xx,yy,z),加密后的地貌图,x=0:.5:5;y=0:.5:6;,z=89 90 87 85 92 91 96 93 90 87 82,92 96 98 99 95 91 89 86 84 82 84,96 98 95 92 90 88 85 84 83 81 85,80 81 82 89 95 96 93 92 89 86 86,82 85 87 98 99 96 97 88 85 82 83,82 85 89 94 95 93 92 91 86 84 88,88 92 93 94 95 89 87 86 83 81 92,92 96 97 98 96 93 95 84 82 81 84,85 85 81 82 80 80 81 85 90 93 95,84 86 81 98 99 98 97 96 95 84 87,80 81 85 82 83 84 87 90 95 86 88,80 82 81 84 85 86 83 82 81 80 82,87 88 89 98 99 97 96 98 94 92 87;,xi=linspace(0,5,50); %,加密横坐标数据到,50,个,yi=linspace(0,6,80); %,加密纵坐标数据到,60,个,xii,yii=meshgrid(xi,yi); %,生成网格数据,zii=interp2(x,y,z,xii,yii,cubic); %,插值,mesh(xii,yii,zii) %,加密后的地貌图,山区地形地貌图 结果,通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。,To MATLAB (moutain),返回,插值函数,griddata,格式为,:,cz,=griddata,（,x,，,y,，,z,，,cx,，,cy,，,method,）,用,MATLAB,作散点数据的插值计算,要求,cx,取行向量，,cy,取为列向量,。,被插值点,插值方法,插值节点,被插值点的函数值,nearest,最邻近插值,linear,双线性插值,cubic,双三次插值,v4- Matlab,提供的插值方法,缺省时,双线性插值, x=-3+6*rand(200,1);y=-2+4*rand(200,1);, z=(x.2-2*x).*exp(-x.2-y.2-x.*y);, x1,y1=meshgrid(-3:.2:3,-2:.2:2);, z1=griddata(x,y,z,x1,y1,cubic);, surf(x1,y1,z1),axis(-3,3,-2,2,-0.7,1.5), z2=griddata(x,y,z,x1,y1,v4);, figure;surf(x1,y1,z2),axis(-3,3,-2,2,-0.7,1.5),例：,在,x,为,3,，,3,，,y,为,2,，,2,矩形区域随机选择一组坐标，用,v4 ,与,cubic,插值法进行处理，并对误差进行比较。,误差分析, z0=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);, surf(x1,y1,abs(z0-z1),axis(-3,3,-2,2,0,0.15), figure;surf(x1,y1,abs(z0-z2),axis(-3,3,-2,2,0,0.15),例：,在,x,为,3,，,3,，,y,为,2,，,2,矩形区域随机选择一组坐标中，对分布不均匀数据，进行插值分析。, x=-3+6*rand(200,1); y=-2+4*rand(200,1);, z=(x.2-2*x).*exp(-x.2-y.2-x.*y); %,生成已知数据, plot(x,y,x) %,样本点的二维分布, figure, plot3(x,y,z,x), axis(-3,3,-2,2,-0.7,1.5),，,grid on,去除在,(-1,-1/2),点为圆心，以,0.5,为半径的圆内的点。, x=-3+6*rand(200,1); y=-2+4*rand(200,1); %,重新生成样本点, z=(x.2-2*x).*exp(-x.2-y.2-x.*y);, ii=find(x+1).2+(y+0.5).20.52); %,找出不在圆内的点坐标, x=x(ii); y=y(ii); z=z(ii); plot(x,y,x), t=0:.1:2*pi,2*pi; x0=-1+0.5*cos(t); y0=-0.5+0.5*sin(t);, line(x0,y0),新样本点拟合曲面, x1,y1=meshgrid(-3:.2:3, -2:.2:2);, z1=griddata(x,y,z,x1,y1,v4);, surf(x1,y1,z1), axis(-3,3,-2,2,-0.7,1.5),误差分析, z0=(x1.2-2*x1).*exp(-x1.2-y1.2-x1.*y1);, surf(x1,y1,abs(z0-z1), axis(-3,3,-2,2,0,0.1), contour(x1,y1,abs(z0-z1),30); hold on, plot(x,y,x); line(x0,y0),误差的二维等高线图,命令,3,interp3,三维网格生成用,meshgrid( ),函数,调用格式,:,x,y,z=meshgrid(x1,y1,z1),其中,x1,y1,z1,为这三维所需要的分割形式，应以向量形式给出，返回,x,y,z,为网格的数据生成，均为三维数组。,griddata3( ),三维非网格形式的插值拟合,命令,4,interpn,n,维网格生成用,ndgrid( ),函数,调用格式,:,x1,x2,xn=ndgridv1,v2,vn,griddatan( ) n,维非网格形式的插值拟合,interp3 ( ),、,interpn( ),调用格式同,interp2(),函数一致；,griddata3( ),、,griddatan( ),调用格式同,griddata( ),函数一致。,例：,通过函数生成一些网格型样本点，试根据样本点进行拟合，并给出拟合误差。, x,y,z=meshgrid(-1:0.2:1); x0,y0,z0=meshgrid(-1:0.05:1);, V=exp(x.2.*z+y.2.*x+z.2.*y,).*cos(x.2.*y.*z+z.2.*y.*x);, V0=exp(x0.2.*z0+y0.2.*x0,+z0.2.*y0).*cos(x0.2.*y0.*z0+z0.2.*y0.*x0);, V1=interp3(x,y,z,V,x0,y0,z0,spline); err=V1-V0; max(err(:),ans =,0.0419,slice(x0,y0,z0,V1,-0.5,0.3, 0.9,0.6,-0.1,-1,-0.5,0.5,1),title(Slives for Four Dim Figures),x,y,z=meshgrid(-3:.1:3,-4:.2:4,-4:.4:4);,v=x.*exp(-x.2-y.2-z.2);,xslice=-2.2,.5,3;,yslice=3,0.9;,zslice=-4,1;,slice(x,y,z,v,xslice,yslice,zslice);,例 在某海域测得一些点,(x,y),处的水深,z,由下表给出，船的吃水深度为,5,英尺，在矩形区域（,75,，,200,）*（,-50,，,150,）里的哪些地方船要避免进入。,To MATLAB hd1,返回,4.,作出水深小于,5,的海域范围,即,z=5,的等高线,.,%,程序一：插值并作海底曲面图,x =129.0 140.0 103.5 88.0 185.5 195.0 105.5 157.5 107.5 77.0 81.0 162.0 162.0 117.5 ;,y = 7.5 141.5 23.0 147.0 22.5 137.5 85.5 -6.5 -81 3.0 56.5 -66.5 84.0 -33.5 ;,z = 4 8 6 8 6 8 8 9 9 8 8 9 4 9 ;,x1=75:1:200;,y1=-50:1:150;,x1,y1=meshgrid(x1,y1);,z1=griddata(x,y,z,x1,y1,v4);,meshc(x1,y1,z1),海底曲面图,%,程序二：插值并作出水深小于,5,的海域范围。,x1=75:1:200;,y1=-50:1:150;,x1,y1=meshgrid(x1,y1);,z1=griddata(x,y,z,x1,y1,v4);,%,插值,z1(z1=5)=nan;,%,将水深大于,5,的置为,nan,，这样绘图就不会显示出来,meshc(x1,y1,z1),水深小于,5,的海域范围,拉格朗日插值法,拉格朗日插值法,是基于基函数的插值方法，插值多项式可表示为,其中 称为,i,次基函数：,在,MATLAB,中编程实现,拉格朗日插值法,函数为：,Language,。,功能：求已知数据点的拉格朗日多项式；,调用格式：,f=,Language(x,y),或,f=,Language(x,y,x0),。,其中，,x,为已知数据点的,x,坐标向量；,y,为已知数据点的,y,坐标向量；,x0,为插值点的,x,坐标；,f,为求得的拉格朗日多项式或,x0,处的插值。,function f = Language(x,y,x0),%,求已知数据点的拉格朗日多项式,%,已知数据点的,x,坐标向量,:x,%,已知数据点的,y,坐标向量,:y,%,为插值点的,x,坐标,:x0,%,求得的拉格朗日多项式或,x0,处的插值,:f,syms t;,if(length(x) = length(y),n = length(x);,else,disp(x,和,y,的维数不相等！,);,return;,end %,检错,f = 0.0;,for(i = 1:n),l = y(i);,for(j = 1:i-1),l = l*(t-x(j)/(x(i)-x(j);,end;,for(j = i+1:n),l = l*(t-x(j)/(x(i)-x(j); %,计算拉格朗日基函数,end;,f = f + l; %,计算拉格朗日插值函数,simplify(f); %,化简,if(i=n),if(nargin = 3),f = subs(f,t,x0); %,计算插值点的函数值,else,f = collect(f); %,将插值多项式展开,f = vpa(f,6); %,将插值多项式的系数化成,6,位精度的小数,end,end,end,例,4-3,根据下表的数据点求出其拉格朗日插值多项式，并计算当,x=1.6,时,y,的值。,解：,x=1 1.2 1.8 2.5 4;, y=0.8415 0.9320 0.9738 0.5985 -0.7568;, f=language(x,y),f =,1.05427*t-.145485e-1*t2-.204917*t3+.328112e-1*t4-.261189e-1, f=language(x,y,1.6),f =,0.9992,x,1,1.2,1.8,2. 5,4,y,0.8415,0.9320,0.9738,0.5985,-0.7568,利用均差的牛顿插值法,函数,f,的零阶均差定义为 ，一阶定义均差为,一般地，函数,f,的,k,阶均差定义为：,利用均差的牛顿插值法,多项式为,系数的计算过程如表所示,表,均差计算表格,一阶均差,二阶均差,三阶均差,n,阶均差,在,MATLAB,中编程实现,利用均差牛顿插值法,函数为：,Newton,。,功能：求已知数据点的均差形式的牛顿插值多项式；,调用格式：,f=,Newton(x,y),或,f=,Newton(x,y,x0),。,其中，,x,为已知数据点的,x,坐标向量；,y,为已知数据点的,y,坐标向量；,x0,为插值点的,x,坐标；,f,为求得的牛顿插值法多项式或,x0,处的插值。,function f = Newton(x,y,x0),%,求已知数据点的均差形式牛顿插值多项式,%,已知数据点的,x,坐标向量,:x,%,已知数据点的,y,坐标向量,:y,%,为插值点的,x,坐标,:x0,%,求得的均差形式牛顿插值多项式或,x0,处的插值,:f,syms t;,if(length(x) = length(y),n = length(x);,c(1:n) = 0.0;,else,disp(x,和,y,的维数不相等！,);,return;,end,f = y(1);,y1 = 0;,l = 1;,for(i=1:n-1),for(j=i+1:n),y1(j) = (y(j)-y(i)/(x(j)-x(i);,end,c(i) = y1(i+1);,l = l*(t-x(i);,f = f + c(i)*l;,simplify(f);,y = y1;,if(i=n-1),if(nargin = 3),f = subs(f,t,x0);,else,f = collect(f); %,将插值多项式展开,f = vpa(f, 6);,end,end,end,例 根据下表的数据点求出其均差形式牛顿插值多项式，并计算当,x=2.0,时,y,的值。,解：, x=1 1.2 1.8 2.5 4;, y=1 1.44 3.24 6.25 16;, f=Newton(x,y),f =,.182711e-14-.482154e-14*t+1.00000*t2-.169177e-14*t3+.211471e-15*t4, f=Newton(x,y,2.0),f =,4,x,1,1.2,1.8,2. 5,4,y,1,1.44,3.24,6.25,16,利用差分的牛顿插值法,差分分为向前差分、后向差分和中心差分三种，它们的记法及定义如下为：,其中： 代表向前差分；,代表向后差分；,代表向后差分。,假设 。为了方便，可构造如表所示的差分表（ ）。,表,差分计算表格,向前牛顿插值,向前牛顿插值,多项式可表示如下：,其中 叫做步长， ， 且 的取值范围为,。,在,MATLAB,中编程实现,向前牛顿插值法,函数为：,Newtonforward,。,功能：求已知数据点的向前牛顿插值法多项式；,调用格式：,f=,Newtonforward(x,y),或,f=,Newtonforward,(x,y,x0),。,其中，,x,为已知数据点的,x,坐标向量；,y,为已知数据点的,y,坐标向量；,x0,为插值点的,x,坐标；,f,为求得的向前牛顿插值法多项式或,x0,处的插值。,function f = Newtonforward(x,y,x0),%,求已知数据点的向前差分牛顿插值多项式,%,已知数据点的,x,坐标向量,:x,%,已知数据点的,y,坐标向量,:y,%,为插值点的,x,坐标,:x0,%,求得的向前差分牛顿插值多项式或,x0,处的插值,:f,syms t;,if(length(x) = length(y),n = length(x);,c(1:n) = 0.0;,else,disp(x,和,y,的维数不相等！,);,return;,end,f = y(1);,y1 = 0;,xx =linspace(x(1),x(n),(x(2)-x(1);,if(xx = x),disp(,节点之间不是等距的！,);,return;,end,for(i=1:n-1),for(j=1:n-i),y1(j) = y(j+1)-y(j);,end,c(i) = y1(1);,l = t;,for(k=1:i-1),l = l*(t-k);,end;,f = f + c(i)*l/factorial(i);,simplify(f);,y = y1;,if(i=n-1),if(nargin = 3),f = subs(f,t,(x0-x(1)/(x(2)-x(1);,else,f = collect(f);,f = vpa(f, 6);,end,end,end,向前牛顿插值,向后牛顿插值,多项式可表示如下：,其中 叫做步长， ， 且 的取值范围为,。,在,MATLAB,中编程实现,向后牛顿插值法,函数为：,Newtonback,。,功能：求已知数据点的向前牛顿插值法多项式；,调用格式：,f=,Newtonback,(x,y),或,f =,Newtonback (x,y,x0),。,其中，,x,为已知数据点的,x,坐标向量；,y,为已知数据点的,y,坐标向量；,x0,为插值点的,x,坐标；,f,为求得的向前牛顿插值法多项式或,x0,处的插值。,function f = Newtonback(x,y,x0),%,求已知数据点的向后差分牛顿插值多项式,%,已知数据点的,x,坐标向量,:x,%,已知数据点的,y,坐标向量,:y,%,为插值点的,x,坐标,:x0,%,求得的向前差分牛顿插值多项式或,x0,处的插值,:f,syms t;,if(length(x) = length(y),n = length(x);,c(1:n) = 0.0;,else,disp(x,和,y,的维数不相等！,);,return;,end,f = y(n);,y1 = 0;,xx =linspace(x(1),x(n),(x(2)-x(1);,if(xx = x),disp(,节点之间不是等距的！,);,return;,end,for(i=1:n-1),for(j=i+1:n),y1(j) = y(j)-y(j-1);,end,c(i) = y1(n);,l = t;,for(k=1:i-1),l = l*(t+k);,end;,f = f + c(i)*l/factorial(i);,simplify(f);,y = y1;,if(i=n-1),if(nargin = 3),f = subs(f,t,(x0-x(n)/(x(2)-x(1);,else,f = collect(f);,f = vpa(f, 6);,end,end,end,例,5,根据下表的数据点求出其差分形式的牛顿插值多项式，并计算当,x=1.55,时,y,的值。,解：,x=1 1.2 1.4 1.6 1.8 ;, y=0.8415 0.93200.9854 0.9996 0.9738;, f=Newtonforward(x,y),f =,.841500+.108025*t-.169042e-1*t2-.675000e-3*t3+.541667e-4*t4, f=Newtonforward(x,y,1.55),f =,0.9998,f=Newtonback(x,y),f =,.973800-.457417e-1*t-.198042e-1*t2+.191667e-3*t3+.541667e-4*t4, f=Newtonback(x,y,1.55),f =,0.9998,x,1,1.2,1.4,1.6,1.8,y,0.8415,0.9320,0.9854,0.9996,0.9738,Hermite,插值,Hermite,插值,满足在节点上等于给定函数值，而且在节点上的导数值也等于给定的导数值。对于高阶导数的情况，,Hermite,插值,多项式比较复杂，在实际中，常常遇到的是函数值与一阶导数给定的情况。在此情况下，,n,个节点,x,1,，,x,2,，,，,x,n,的,Hermite,插值,多项式的表达式如下：,其中 ， ，,，,在,MATLAB,中编程实现,Hermite,插值法,函数为：,Hermite,。,功能：求已知数据点的,Hermite,插值法多项式；,调用格式：,f=,Hermite,(x,y,y_1),或,f =,Hermite,(x,y, y_1,x0),。,其中，,x,为已知数据点的,x,坐标向量；,y,为已知数据点的,y,坐标向量；,y_1,为已知数据点导数向量；,x0,为插值点的,x,坐标；,f,为求得的,Hermite,插值法多项式或,x0,处的插值。,function f = Hermite(x,y,y_1,x0),%,求已知数据点的向后差分牛顿插值多项式,%,已知数据点的,x,坐标向量,:x,%,已知数据点的,y,坐标向量,:y,%,已知数据点的导数向量,:y_1,%,求得的,Hermite,插值多项式或,x0,处的插值,:f,syms t;,f = 0.0;,if(length(x) = length(y),if(length(y) = length(y_1),n = length(x);,else,disp(y,和,y,的导数的维数不相等！,);,return;,end,else,disp(x,和,y,的维数不相等！,);,return;,end,for i=1:n,h = 1.0;,a = 0.0;,for j=1:n,if( j = i),h = h*(t-x(j)2/(x(i)-x(j)2);,a = a + 1/(x(i)-x(j);,end,end,f = f + h*(x(i)-t)*(2*a*y(i)-y_1(i)+y(i);,if(i=n),if(nargin = 4),f = subs(f,t,x0);,else,f = vpa(f,6);,end,end,end,例,6,根据下表的数据点求出其,Hermite,插值多项式，并计算当,x=1.144,时,y,的值。,解：, x=1:0.2:1.8;, y=1 1.0954 1.1832 1.2649 1.3416;, y_1=0.5 0.4564 0.4226 0.3953 0.3727;, f=Hermite(x,y,y_1),f =,678.168*(t-1.20000)2*(t-1.40000)2*(t-1.60000)2*(t-1.80000)2*(-20.3333+21.3333*t)+10850.7*(t-1.)2*(t-1.40000)2*(t-1.60000)2*(t-1.80000)2*(-10.4063+9.58473*t)+24414.1*(t-1.)2*(t-1.20000)2*(t-1.60000)2*(t-1.80000)2*(.591560+.422600*t)+10850.7*(t-1.)2*(t-1.20000)2*(t-1.40000)2*(t-1.80000)2*(17.4978-10.1455*t)+678.168*(t-1.)2*(t-1.20000)2*(t-1.40000)2*(t-1.60000)2*(50.9807-27.5773*t), f=Hermite(x,y,y_1,1.44),f =,1.2000,x,1,1.2,1.4,1.6,1.8,y,1,1.0954,1.1832,1.2649,1.3416,y,0.5000,0.4564,0.4226,0.3953,0.3727,spline,三次样条插值,Spline,插值为分段三次插值，即：在每一个小区间上是不超过三次的多项式且具有二阶连续导数，利用具有一阶导数边界条件的源程序为：,naspline.m,function m=naspline(x,y,dy0,dyn,xx),%,用途：三阶样条插值（一阶导数边界条件）,%,格式：,x,为节点向量；,y,为数据；,dyo,，,dyn,为左右两端点的一阶导数,%,如果,xx,缺省，则输出各节点的一阶导数值，否则,m,为,xx,的三阶样条插值,n=length(x)-1; %,计算小区间的个数,h=diff(x);,lemda=h(2:n)./(h(1:n-1)+h(2:n);,mu=1-lemda;,g=3*(lemda.*diff(y(1:n)./h(1:n-1)+mu.*diff(y(2:n+1)./h(2:n);,g(1)=g(1)-lemda(1)*dy0;,g(n-1)=g(n-1)-mu(n-1)*dyn;,%,求解三对角方程,dy=nachase(lemda,2*ones(1:n-1),mu,g);,%,若给插值节点，计算插值,m=dy0;dy;dyn;,if nargin=5,s=zeros(size(xx);,for i=1:n,if i=1,kk=find(xx=x(n);,else,kk=find(xxx(i),end,xbar=(xx(kk)-x(i)/h(i);,s(kk)=alpha0(xbar)*y(i)+alpha1(xbar)*y(i+1)+.,h(i)*beta0(xbar)*m(i)+h(i)*beta1(xbar)*m(i+1);,end,m=s;,end,%,追赶法,function x=nachase(a,b,c,d),n=length(a);,for k=2:n,b(k)=b(k)-a(k)/b(k-1)*c(k-1);,d(k)=d(k)-a(k)/b(k-1)*d(k-1);,end,x(n)=d(n)/b(n);,for k=n-1:-1:1,x(k)=(d(k)-c(k)*x(k+1)/b(k);,end,x=x(:);,%,基函数,function y=alpha0(x),y=2*x.3-3*x.2+1;,function y=alpha1(x),y=-2*x.3+3*x.2;,function y=beta0(x),y=x.3-2*x.2+x;,function y=beta1(x),y=x.3-x.2;,