matlab数值分析插值.ppt

第三章第三章 插值插值 插值插值：插值就是定义一个在特定点取给定值的函数的过程本章的重点是介绍两个紧密相关的插值函数：分段三次样条函数和保形分段三次插值函数（称为“pchip”）matlab数值分析插值•3.1 插值多项式插值多项式•3.2 分段线性插值分段线性插值•3.3 分段三次埃米特插值分段三次埃米特插值•3.4 保形分段三次插值保形分段三次插值•3.5 三次样条三次样条•3.6 pchiptx，，splinetx•3.7interpguimatlab数值分析插值3.1插值多项式插值多项式 平面上的任意两点（平面上的任意两点（x x1 1，，y y1 1）和（）和（x x2 2，，y y2 2））, ,只只要要x x1 1≠x≠x2 2, ,就为以确定一个关于就为以确定一个关于x x的一次多项式的一次多项式, ,其图形经过这两点对于这个多项式，有多种其图形经过这两点对于这个多项式，有多种不同的公式表达，但它们都对应同一个直线图不同的公式表达，但它们都对应同一个直线图形matlab数值分析插值两个点时，假定给定区间两个点时，假定给定区间[x[xk k,x,xk+1k+1] ]及端点函数值及端点函数值y yk k=f=f（（x xk k），），y yk+1k+1=f=f（（x xk+1k+1），要求线性插值多项式），要求线性插值多项式L L1 1（（x x），使它满足），使它满足 L L1 1（（x xk k））=y=yk k，，L L1 1（（x xk+1k+1））=y=yk+1k+1 y=L y=L1 1（（x x）的几何意义就是通过两点）的几何意义就是通过两点(x(xk k，，y yk k) )与与（（x xk+1 k+1 ，，y yk+1k+1) ) 的直线，的直线， L L1 1（（x x）的表达式可由几何意义）的表达式可由几何意义直接给出直接给出 matlab数值分析插值•由两点式可以看出，由两点式可以看出，L1（（x）是由两个线性函数）是由两个线性函数•的线性组合得到，其系数分别为的线性组合得到，其系数分别为yk和和yk+1，即，即•显然，显然， lk（（x）及）及lk+1（（x）也是线性插值多项式，）也是线性插值多项式，在节点在节点xk及及xk+1上满足条件上满足条件•我们称函数我们称函数lk（（x）与）与lk+1（（x）为线性插值基函数。

为线性插值基函数matlab数值分析插值•这种用插值基函数表示的方法推广到一般情形，以这种用插值基函数表示的方法推广到一般情形，以下讨论如何构造通过下讨论如何构造通过n+1个节点个节点x0

的多项式，使其图形经过这些点 很容易看出，数据点的数目很容易看出，数据点的数目n+1n+1也是多项式系数也是多项式系数的个数尽管，一些首项的系数可能是零，但多项式的个数尽管，一些首项的系数可能是零，但多项式的次数实际上也小于的次数实际上也小于n n同样，这个多项式，有多种同样，这个多项式，有多种不同的公式表达，但它们都对应同一个直线图形不同的公式表达，但它们都对应同一个直线图形matlab数值分析插值 这样的多项式称为这样的多项式称为插值插值多项式，它可以准确的多项式，它可以准确的重新计算出初始给定的数据：重新计算出初始给定的数据： 表示插值多项式的最紧凑的方式是拉格朗日形式表示插值多项式的最紧凑的方式是拉格朗日形式matlab数值分析插值例如，考虑下面一组数据例如，考虑下面一组数据x=0:3;x=0:3;y=[-5 -6 -1 16];y=[-5 -6 -1 16];输入命令输入命令disp([x;y])disp([x;y])其输出为其输出为 0 1 2 30 1 2 3-5 -6 -1 16-5 -6 -1 16这些数据的拉格朗日形式的多项式插值为这些数据的拉格朗日形式的多项式插值为matlab数值分析插值•一个多项式通常不用拉格朗日形式表示，它更常见的写成类似•的形式。

其中简单的x的次方项称为单项式单项式，而多项式的这种形式称为使用幂形式幂形式的多项式•插值多项式使用幂形式表示为matlab数值分析插值•其中的系数，原则上可以通过求解下面的线性代数其中的系数，原则上可以通过求解下面的线性代数方程组得到方程组得到•这个线性方程组的系数矩阵记为这个线性方程组的系数矩阵记为V，也被称为范德，也被称为范德尔蒙（尔蒙（Vandermonde）矩阵，该矩阵的各个元素）矩阵，该矩阵的各个元素为为 •上述范德尔蒙矩阵的各列，有时也按相反的顺序排上述范德尔蒙矩阵的各列，有时也按相反的顺序排列，但在列，但在MATLAB中，多项式系数向量，通常按中，多项式系数向量，通常按从高次幂到低次幂排列从高次幂到低次幂排列matlab数值分析插值•MATLAB中的函数中的函数vander可生成范德尔蒙矩阵，可生成范德尔蒙矩阵，例如对于前面的那组数据，例如对于前面的那组数据，V=vander(x)V=vander(x)生成生成V =V = 0 0 0 1 0 0 0 1 1 1 1 1 1 1 1 1 8 4 2 1 8 4 2 1 27 9 3 1 27 9 3 1然后，输入命令然后，输入命令 c=V\y' c=V\y'计算出插值系数计算出插值系数c =c = 1.0000 1.0000 0.0000 0.0000 -2.0000 -2.0000 -5.0000 -5.0000matlab数值分析插值能实现各种插值算法的能实现各种插值算法的MATLAB函数函数它们都采用下面的调用格式它们都采用下面的调用格式 v=interp v=interp（（x x，，y y，，u u））前两个参数，前两个参数，x x和和y y，是长度相同的向量，它们定义了插值，是长度相同的向量，它们定义了插值点。

第三个参数点第三个参数u u，为要计算函数值的范围上的点组成的，为要计算函数值的范围上的点组成的向量输出向量向量输出向量v v和和u u长度相等，长度相等，其分量其分量v v（（k k））=interp=interp（（x x，，y y，，u u（（k k））第一个这样的插值函数是第一个这样的插值函数是polyinterppolyinterp，它基于拉格朗日形式，它基于拉格朗日形式matlab数值分析插值•为了解释为了解释polyinterp函数的功能，先构造一个间隔很函数的功能，先构造一个间隔很密的求值点向量密的求值点向量•u = -.25:.01:3.25;•然后输入命令然后输入命令•v = polyinterp(x,y,u);•plot(x,y,'o',u,v,'-')•可生成图可生成图3-1matlab数值分析插值•函数函数polyinterppolyinterp也可以处理符号变量，例如创建符号也可以处理符号变量，例如创建符号变量变量 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的幂形式的幂形式 P = x^3-2*x-5计算机显示插值多项式的符号形式matlab数值分析插值另外一个例子，使用的是本章另一种方法所用的例子另外一个例子，使用的是本章另一种方法所用的例子x = 1:6;x = 1:6;y = [16 18 21 17 15 12];y = [16 18 21 17 15 12];disp([x; y])disp([x; y])u = .75:.05:6.25;u = .75:.05:6.25;v = polyinterp(x,y,u);v = polyinterp(x,y,u);plot(x,y,'o',u,v,'-');plot(x,y,'o',u,v,'-');其运行后的结果为其运行后的结果为 1 2 3 4 5 61 2 3 4 5 6 16 18 21 17 15 12 16 18 21 17 15 12同时输出同时输出3-2。

matlab数值分析插值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,'-');matlab数值分析插值•在生成图在生成图3-33-3所示的图线时，所示的图线时，MATLAB图像处理函数使图像处理函数使用了分段线性插值用了分段线性插值•这个分段线性插值算法是其他更复杂算法的基础，他这个分段线性插值算法是其他更复杂算法的基础，他用了三个量用了三个量•首先要确定间隔序号（首先要确定间隔序号（interval index））k，使得，使得• 第二个量是局部变量（第二个量是局部变量（local variable））s，其定义为，其定义为•最后一个量是一次均差（最后一个量是一次均差（first divided difference））•定义了这三个量，则插值基函数可表示为定义了这三个量，则插值基函数可表示为• 显然，这是通过点(xk; yk)和(xk+1; yk+1) 点的线性函数matlab数值分析插值•点点x有时也被称为有时也被称为断点断点。

有上述基函数构成的分段线有上述基函数构成的分段线性插值基函数性插值基函数L(x)是关于是关于x的连续函数，但它的一阶导的连续函数，但它的一阶导数数L’(x)，则不连续在每个，则不连续在每个x的子区间上，导数值为的子区间上，导数值为常数常数δδ，但在断点上，它的值发生跳变但在断点上，它的值发生跳变•用用piecelin.mpiecelin.m函数可实现分段线性插值，输入的参数函数可实现分段线性插值，输入的参数u u，可以是需要计算的点构成的向量下标，可以是需要计算的点构成的向量下标k k实际上是一实际上是一个由序号组成的向量个由序号组成的向量matlab数值分析插值matlab数值分析插值3.3分段三次埃米特插值分段三次埃米特插值•许多最有效的插值技术都基于分段三次多项式令许多最有效的插值技术都基于分段三次多项式令hk为第为第k段子区间的长度：段子区间的长度：•那么一次均差那么一次均差δδk k由下面的公式给出由下面的公式给出•令令dk为插值基函数在点为插值基函数在点xk处的斜率，即：处的斜率，即：•对于分段线性插值基函数，对于分段线性插值基函数，dk=δδk k或或δδ k+1，但对于，但对于更高次的插值多项式不一定成立。

更高次的插值多项式不一定成立matlab数值分析插值•考虑一个定义在区间考虑一个定义在区间xk≤x ≤xk+1的函数，采用局的函数，采用局部变量部变量s= x -xk并令并令h=hk，它可表示为：，它可表示为：•这是一个关于，也即的三次多项式它满足四个这是一个关于，也即的三次多项式它满足四个插值条件，其中两个关于函数值，两个关于函数插值条件，其中两个关于函数值，两个关于函数的导数值：的导数值：matlab数值分析插值•那些满足关于导数值插值条件的函数称为那些满足关于导数值插值条件的函数称为埃米埃米特特（（hermite）或）或密切密切（（osculatory）插值基）插值基函数，因为这些函数在插值点上保持高阶的连函数，因为这些函数在插值点上保持高阶的连续性（在拉丁文中续性（在拉丁文中“密切密切”一词的本意为一词的本意为“亲亲吻吻”）•如果正好给定了一系列数据点上的函数值和一如果正好给定了一系列数据点上的函数值和一阶导数值，那么就可以用埃米特插值拟合这些阶导数值，那么就可以用埃米特插值拟合这些数据但是如果没有给出这些导数值，那么需数据但是如果没有给出这些导数值，那么需要用一些方法来限定斜率要用一些方法来限定斜率dk，我们在下一节中，我们在下一节中讨论两种可能的办法，即在讨论两种可能的办法，即在MATLAB中的函中的函数数pchip和和spline。

matlab数值分析插值3.4 保形分段三次插值保形分段三次插值•pchip实际是实际是“分段三次埃米特插值多项式分段三次埃米特插值多项式”（（piecewise cubic Hermite interpolating polynominal））的的英文首字母缩写有意思的是，根据这个名字并不能英文首字母缩写有意思的是，根据这个名字并不能确定它到底是哪一种分段三次埃米特插值多项式，因确定它到底是哪一种分段三次埃米特插值多项式，因为样条插值函数实际也是分段三次埃米特插值多项式，为样条插值函数实际也是分段三次埃米特插值多项式，只是对斜率的限制条件不同而已只是对斜率的限制条件不同而已•在这里，我们说的在这里，我们说的pchip实际上是一个最近才引入实际上是一个最近才引入MATLAB、保形的（、保形的（shape-preserving）且看上去不）且看上去不错的特定插值函数它基于一个由错的特定插值函数它基于一个由Fritsch和和Carlson编写的旧的编写的旧的Fortran程序，在程序，在Kahaner、、Moler和和Nash的书【的书【33】中可以找到相关的介绍中可以找到相关的介绍matlab数值分析插值•对于前面的那个例子数据，图对于前面的那个例子数据，图3-4显示了显示了pchip插值出插值出来的结果。

来的结果matlab数值分析插值 关键思想是如何确定斜率关键思想是如何确定斜率dk，使得函数值不会过度地，使得函数值不会过度地偏离（至少在局部）给定的数据偏离（至少在局部）给定的数据 第一，如果第一，如果δδk k和和δδ k-1和的正负号相反，或者他们和的正负号相反，或者他们中有一个为零，那么在处函数为离散的极大或极小，中有一个为零，那么在处函数为离散的极大或极小，于是可以令于是可以令dk=0 关于它的解释可见右图关于它的解释可见右图 图中实线为分段线性插值，它在中间断点两侧的斜率图中实线为分段线性插值，它在中间断点两侧的斜率符号相反因此，图中虚线斜率为零图中的曲线为符号相反因此，图中虚线斜率为零图中的曲线为由两个三次多项式组成的保形插值函数，这两个三次由两个三次多项式组成的保形插值函数，这两个三次多项式在中间断点处相接，在那一点，两条曲线的导多项式在中间断点处相接，在那一点，两条曲线的导数都为零然而，在断点处的二阶导数值存在跳变数都为零然而，在断点处的二阶导数值存在跳变matlab数值分析插值 第二，如果第二，如果δδk k和和δδ k-1和的正负号相同，并且两个子和的正负号相同，并且两个子区间长度相等，则区间长度相等，则dk令为两侧两个斜率的调和平均数：令为两侧两个斜率的调和平均数： 这中埃米特插值函数，在断点处斜率的倒数为两侧这中埃米特插值函数，在断点处斜率的倒数为两侧分段线性插值函数斜率导数的平均。

这种情况如下图分段线性插值函数斜率导数的平均这种情况如下图所示 在断点处，分段线性插值函数的斜率的倒数从在断点处，分段线性插值函数的斜率的倒数从1到到变到变到5，因此图中虚线斜率的倒数为，因此图中虚线斜率的倒数为1和和5的平均值，的平均值，即即3这个保形插值函数由两个三次多项式组成，它这个保形插值函数由两个三次多项式组成，它们在中间断点处相接，并且在那一点处的导数都为们在中间断点处相接，并且在那一点处的导数都为1/3matlab数值分析插值•如果如果δδk k和和δδk-1和的正负号相同，并且两个子区和的正负号相同，并且两个子区间长度不等，则间长度不等，则dk为加权的调和平均，权重由为加权的调和平均，权重由两个子区间的长度决定：两个子区间的长度决定：•其中，其中，•上面介绍了函数上面介绍了函数pchip中如何确定中间断点处的中如何确定中间断点处的斜率，而对于整个数据区间的两个端点处的斜斜率，而对于整个数据区间的两个端点处的斜率率d1和和dn，需要用一个稍许不同的、单方向分，需要用一个稍许不同的、单方向分析的方法加以计算，有关细节请参考文件析的方法加以计算，有关细节请参考文件pchiptx.m。

matlab数值分析插值。