matlab数值分析第三章插值.ppt

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第三章 插值,插值,：,插值就是定义一个在特定点取给定值的函数的过程。本章的重点是介绍两个紧密相关的插值函数：分段三次样条函数和保形分段三次插值函数（称为,“,pchip,”,）,3.1,插值多项式,3.2,分段线性插值,3.3,分段三次埃米特插值,3.4,保形分段三次插值,3.5,三次样条,3.6 pchiptx,，,splinetx,3.7interpgui,3.1,插值多项式,平面上的任意两点（,x,1,，,y,1,）和（,x,2,，,y,2,）,只要,x,1,x,2,就为以确定一个关于,x,的一次多项式,其图形经过这两点。对于这个多项式，有多种不同的公式表达，但它们都对应同一个直线图形。,两个点时，假定给定区间,x,k,x,k+1,及端点函数值,y,k,=f,（,x,k,），,y,k+1,=f,（,x,k+1,），要求线性插值多项式,L,1,（,x,），使它满足,L,1,（,x,k,）,=y,k,，,L,1,（,x,k+1,）,=y,k+1,y=L,1,（,x,）的几何意义就是通过两点,(x,k,，,y,k,),与,（,x,k+1,，,y,k+1,),的直线，,L,1,（,x,）的表达式可由几何意义直接给出,由两点式可以看出，,L,1,（,x,）是由两个线性函数,的线性组合得到，其系数分别为,y,k,和,y,k+1,，即,显然，,l,k,（,x,）及,l,k+1,（,x,）也是线性插值多项式，在节点,x,k,及,x,k+1,上满足条件,我们称函数,l,k,（,x,）与,l,k+1,（,x,）为线性插值基函数。,这种用插值基函数表示的方法推广到一般情形，以下讨论如何构造通过,n+1,个节点,x,0,x,1,x,n,的,n,次插值多项式,Ln,（,x,），假定它满足条件,若,n,次多项式,lj,（,x,）（,j=0,1,n,）在,n+1,个节点,x,0,x,1,x,n,上满足条件,就称这,n+1,个,n,次多项式,l,0,(x),l,1,(x),l,n,(x),为节点,x,0,，,x,1,，,，,x,n,上的,n,次插值基函数,。,用类似的推导方法，可得到,n,次插值基函数为,显然它满足条件（,1,）式。于是满足条件（,1,）的插值多项式,Ln,（,x,）可表示为,由,l,k,（,x,）的定义知,形如（,3,）式的插值多项式,Ln,（,x,）称为,拉格朗日插值多项式。,则对于平面上有着不同,x,k,值的,n+1,个点,（,x,k,，,y,k,），,k=0,，,1,n,存在唯一一个关于,x,的次数小于,n+1,的多项式，使其图形经过这些点。,很容易看出，数据点的数目,n+1,也是多项式系数的个数。尽管，一些首项的系数可能是零，但多项式的次数实际上也小于,n,。同样，这个多项式，有多种不同的公式表达，但它们都对应同一个直线图形。,这样的多项式称为,插值,多项式，它可以准确的重新计算出初始给定的数据：,表示插值多项式的最紧凑的方式是拉格朗日形式,例如，考虑下面一组数据,x=0:3;,y=-5-6-1 16;,输入命令,disp(x;y),其输出为,0 1 2 3,-5 -6 -1 16,这些数据的拉格朗日形式的多项式插值为,一个多项式通常不用拉格朗日形式表示，它更常见的写成类似,的形式。其中简单的,x,的次方项称为,单项式,，而多项式的这种形式称为使用,幂形式,的多项式。,插值多项式使用幂形式表示为,其中的系数，原则上可以通过求解下面的线性代数方程组得到。,这个线性方程组的系数矩阵记为,V,，也被称为范德尔蒙（,Vandermonde,）矩阵，该矩阵的各个元素为,上述范德尔蒙矩阵的各列，有时也按相反的顺序排列，但在,MATLAB,中，多项式系数向量，通常按从高次幂到低次幂排列。,MATLAB,中的函数,vander,可生成范德尔蒙矩阵，例如对于前面的那组数据，,V=vander(x),生成,V=,0 0 0 1,1 1 1 1,8 4 2 1,27 9 3 1,然后，输入命令,c=Vy,计算出插值系数,c=,1.0000,0.0000,-2.0000,-5.0000,能实现各种插值算法的,MATLAB,函数,它们都采用下面的调用格式,v=interp,（,x,，,y,，,u,）,前两个参数，,x,和,y,，是长度相同的向量，它们定义了插值,点。第三个参数,u,，为要计算函数值的范围上的点组成的,向量。输出向量,v,和,u,长度相等，,其分量,v,（,k,）,=interp,（,x,，,y,，,u,（,k,）。,第一个这样的插值函数是,polyinterp,，它基于拉格朗日形式,。,为了解释,polyinterp,函数的功能，先构造一个间隔很密的求值点向量。,u=-.25:.01:3.25;,然后输入命令,v=polyinterp(x,y,u);,plot(x,y,o,u,v,-),可生成图,3-1,。,函数,polyinterp,也可以处理符号变量，例如创建符号变量,symx=sym(x),命令：,P=polyinterp(x,y,symx),pretty(P),其输出结果为,-5(-1/3 x+1)(-1/2 x+1)(-x+1)-6(-1/2 x+3/2)(-x+2)x,-1/2(-x+3)(x-1)x+16/3(x-2)(1/2 x-1/2)x,这个表达式是插值多项式的拉格朗日形式。,命令：,P=simplify(P),将其进行简化，从而得到,P,的幂形式,P=,x3-2*x-5,计算机显示插值多项式的符号形式,另外一个例子，使用的是本章另一种方法所用的例子,x=1:6;,y=16 18 21 17 15 12;,disp(x;y),u=.75:.05:6.25;,v=polyinterp(x,y,u);,plot(x,y,o,u,v,-);,其运行后的结果为,1 2 3 4 5 6,16 18 21 17 15 12,同时输出,3-2,。,3.2,分段线性插值,通过两步操作可以绘制出一个简单的图形：,第一步用圆圈在坐标系中标出个数据点,plot(x,y,o);,，,第二步用直线段依次连接这些数据点,plot(x,y-);,。,下面的语句执行这样的操作，生成图,3-3.,x=1:6;,y=16 18 21 17 15 12;,plot(x,y,o,x,y,-);,在生成图,3-3,所示的图线时，,MATLAB,图像处理函数使用了分段线性插值。,这个分段线性插值算法是其他更复杂算法的基础，他用了三个量。,首先要确定间隔序号（,interval index,）,k,，使得,第二个量是局部变量（,local variable,）,s,，其定义为,最后一个量是一次均差（,first divided difference,）,定义了这三个量，则插值基函数可表示为,显然，这是通过点,(,x,k,;y,k,),和,(,x,k,+1,;y,k,+1,),点的线性函数,点,x,有时也被称为,断点,。有上述基函数构成的分段线性插值基函数,L(x),是关于,x,的连续函数，但它的一阶导数,L(x),，则不连续。在每个,x,的子区间上，导数值为常数,，但在断点上，它的值发生跳变。,用,piecelin.m,函数可实现分段线性插值，输入的参数,u,，可以是需要计算的点构成的向量。下标,k,实际上是一个由序号组成的向量。,3.3,分段三次埃米特插值,许多最有效的插值技术都基于分段三次多项式。令,h,k,为第,k,段子区间的长度：,那么一次均差,k,由下面的公式给出,令,d,k,为插值基函数在点,x,k,处的斜率，即：,对于分段线性插值基函数，,d,k,=,k,或,k+1,，但对于更高次的插值多项式不一定成立。,考虑一个定义在区间,x,k,x x,k+1,的函数，采用局部变量,s=x-x,k,并令,h=h,k,，它可表示为：,这是一个关于，也即的三次多项式。它满足四个插值条件，其中两个关于函数值，两个关于函数的导数值：,那些满足关于导数值插值条件的函数称为,埃米特,（,hermite,）或,密切,（,osculatory,）插值基函数，因为这些函数在插值点上保持高阶的连续性（在拉丁文中“密切”一词的本意为“亲吻”）。,如果正好给定了一系列数据点上的函数值和一阶导数值，那么就可以用埃米特插值拟合这些数据。但是如果没有给出这些导数值，那么需要用一些方法来限定斜率,d,k,，我们在下一节中讨论两种可能的办法，即在,MATLAB,中的函数,pchip,和,spline,。,3.4,保形分段三次插值,pchip,实际是“分段三次埃米特插值多项式,”（,piecewise cubic Hermite interpolating polynominal,）,的英文首字母缩写。有意思的是，根据这个名字并不能确定它到底是哪一种分段三次埃米特插值多项式，因为样条插值函数实际也是分段三次埃米特插值多项式，只是对斜率的限制条件不同而已。,在这里，我们说的,pchip,实际上是一个最近才引入,MATLAB,、保形的（,shape-preserving,）且看上去不错的特定插值函数。它基于一个由,Fritsch,和,Carlson,编写的旧的,Fortran,程序，在,Kahaner,、,Moler,和,Nash,的书,【33】,中可以找到相关的介绍。,对于前面的那个例子数据，图,3-4,显示了,pchip,插值出来的结果。,关键思想是如何确定斜率,d,k,，使得函数值不会过度地偏离（至少在局部）给定的数据。,第一，如果,k,和,k-1,和的正负号相反，或者他们中有一个为零，那么在处函数为离散的极大或极小，于是可以令,d,k,=0,关于它的解释可见右图,图中实线为分段线性插值，它在中间断点两侧的斜率符号相反。因此，图中虚线斜率为零。图中的曲线为由两个三次多项式组成的保形插值函数，这两个三次多项式在中间断点处相接，在那一点，两条曲线的导数都为零。然而，在断点处的二阶导数值存在跳变。,第二，如果,k,和,k-1,和的正负号相同，并且两个子区间长度相等，则,d,k,令为两侧两个斜率的调和平均数：,这中埃米特插值函数，在断点处斜率的倒数为两侧分段线性插值函数斜率导数的平均。这种情况如下图所示。,在断点处，分段线性插值函数的斜率的倒数从,1,到变到,5,，因此图中虚线斜率的倒数为,1,和,5,的平均值，即,3,。这个保形插值函数由两个三次多项式组成，它们在中间断点处相接，并且在那一点处的导数都为,1/3,。,如果,k,和,k-1,和的正负号相同，并且两个子区间长度不等，则,d,k,为加权的调和平均，权重由两个子区间的长度决定：,其中，,上面介绍了函数,pchip,中如何确定中间断点处的斜率，而对于整个数据区间的两个端点处的斜率,d,1,和,d,n,，需要用一个稍许不同的、单方向分析的方法加以计算，有关细节请参考文件,pchiptx.m,。,