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

第三章 插值 插值：插值就是定义一个在特定 点取给定值的函数的过程。本章的重 点是介绍两个紧密相关的插值函数： 分段三次样条函数和保形分段三次插 值函数（称为“pchip”） 3.1 插值多项式 3.2 分段线性插值 3.3 分段三次埃米特插值 3.4 保形分段三次插值 3.5 三次样条 3.6 pchiptx，splinetx 3.7interpgui 3.1插值多项式 平面上的任意两点（x1，y1）和（x2，y2）,只 要x1x2,就为以确定一个关于x的一次多项式, 其图形经过这两点。对于这个多项式，有多种 不同的公式表达，但它们都对应同一个直线图 形。 两个点时，假定给定区间xk,xk+1及端点函数值 yk=f（xk），yk+1=f（xk+1），要求线性插值多项式 L1（x），使它满足 L1（xk）=yk，L1（xk+1）=yk+1 y=L1（x）的几何意义就是通过两点(xk，yk)与 （xk+1 ，yk+1) 的直线， L1（x）的表达式可由几何意义 直接给出 由两点式可以看出，L1（x）是由两个线性函数 的线性组合得到，其系数分别为yk和yk+1，即 显然， lk（x）及lk+1（x）也是线性插值多项式， 在节点xk及xk+1上满足条件 我们称函数lk（x）与lk+1（x）为线性插值基函数。 这种用插值基函数表示的方法推广到一般情形，以 下讨论如何构造通过n+1个节点x0x1xn的n次 插值多项式Ln（x），假定它满足条件 若n次多项式lj（x）（j=0,1,n）在n+1个节点 x0x1xn 上满足条件 就称这n+1个n次多项式l0(x), l1(x), ln(x)为节点 x0，x1，xn上的n次插值基函数。 用类似的推导方法，可得到n次插值基函数为 显然它满足条件（1）式。于是满足条件（1）的插值 多项式Ln（x）可表示为 由lk（x）的定义知 形如（3）式的插值多项式Ln（x）称为拉格朗日插值 多项式。 则对于平面上有着不同xk值的n+1个点,（xk，yk ）， 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） 定义了这三个量，则插值基函数可表示为 显然，这是通过 点(xk; yk)和 (xk+1; yk+1) 点的 线性函数 点x有时也被称为断点。有上述基函数构成的分段线 性插值基函数L(x)是关于x的连续函数，但它的一阶导 数L(x)，则不连续。在每个x的子区间上，导数值为 常数，但在断点上，它的值发生跳变。 用piecelin.m函数可实现分段线性插值，输入的参数u ，可以是需要计算的点构成的向量。下标k实际上是一 个由序号组成的向量。 3.3分段三次埃米特插值 许多最有效的插值技术都基于分段三次多项式。令 hk为第k段子区间的长度： 那么一次均差k由下面的公式给出 令dk为插值基函数在点xk处的斜率，即： 对于分段线性插值基函数，dk=k或 k+1，但对 于更高次的插值多项式不一定成立。 考虑一个定义在区间xkx xk+1的函数，采用局 部变量s= x -xk并令h=hk，它可表示为： 这是一个关于，也即的三次多项式。它满足四个 插值条件，其中两个关于函数值，两个关于函数 的导数值： 那些满足关于导数值插值条件的函数称为埃米 特（hermite）或密切（osculatory）插值基 函数，因为这些函数在插值点上保持高阶的连 续性（在拉丁文中“密切”一词的本意为“亲吻” ）。 如果正好给定了一系列数据点上的函数值和一 阶导数值，那么就可以用埃米特插值拟合这些 数据。但是如果没有给出这些导数值，那么需 要用一些方法来限定斜率dk，我们在下一节中 讨论两种可能的办法，即在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插值出 来的结果。 关键思想是如何确定斜率dk，使得函数值不会过度地 偏离（至少在局部）给定的数据。 第一，如果k和 k-1和的正负号相反，或者他们 中有一个为零，那么在处函数为离散的极大或极小， 于是可以令dk=0 关于它的解释可见右图 图中实线为分段线性插值，它在中间断点两侧的斜率 符号相反。因此，图中虚线斜率为零。图中的曲线为 由两个三次多项式组成的保形插值函数，这两个三次 多项式在中间断点处相接，在那一点，两条曲线的导 数都为零。然而，在断点处的二阶导数值存在跳变。 第二，如果k和 k-1和的正负号相同，并且两个子 区间长度相等，则dk令为两侧两个斜率的调和平均数 ： 这中埃米特插值函数，在断点处斜率的倒数为两侧 分段线性插值函数斜率导数的平均。这种情况如下图 所示。 在断点处，分段线性插值函数的斜率的倒数从1到 变到5，因此图中虚线斜率的倒数为1和5的平均值，