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

第三章插值 插值 插值就是定义一个在特定点取给定值的函数的过程 本章的重点是介绍两个紧密相关的插值函数 分段三次样条函数和保形分段三次插值函数 称为 pchip 3 1插值多项式3 2分段线性插值3 3分段三次埃米特插值3 4保形分段三次插值3 5三次样条3 6pchiptx splinetx3 7interpgui 3 1插值多项式 平面上的任意两点 x1 y1 和 x2 y2 只要x1 x2 就为以确定一个关于x的一次多项式 其图形经过这两点 对于这个多项式 有多种不同的公式表达 但它们都对应同一个直线图形 两个点时 假定给定区间 xk xk 1 及端点函数值yk f xk yk 1 f xk 1 要求线性插值多项式L1 x 使它满足L1 xk yk L1 xk 1 yk 1y 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个节点x0 x1 xn的n次插值多项式Ln x 假定它满足条件若n次多项式lj x j 0 1 n 在n 1个节点x0 x1 xn上满足条件就称这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 116 输入命令disp x y 其输出为0123 5 6 116这些数据的拉格朗日形式的多项式插值为 一个多项式通常不用拉格朗日形式表示 它更常见的写成类似的形式 其中简单的x的次方项称为单项式 而多项式的这种形式称为使用幂形式的多项式 插值多项式使用幂形式表示为 其中的系数 原则上可以通过求解下面的线性代数方程组得到 这个线性方程组的系数矩阵记为V 也被称为范德尔蒙 Vandermonde 矩阵 该矩阵的各个元素为上述范德尔蒙矩阵的各列 有时也按相反的顺序排列 但在MATLAB中 多项式系数向量 通常按从高次幂到低次幂排列 MATLAB中的函数vander可生成范德尔蒙矩阵 例如对于前面的那组数据 V vander x 生成V 00011111842127931然后 输入命令c V y 计算出插值系数c 1 00000 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 3x 1 1 2x 1 x 1 6 1 2x 3 2 x 2 x 1 2 x 3 x 1 x 16 3 x 2 1 2x 1 2 x这个表达式是插值多项式的拉格朗日形式 命令 P simplify P 将其进行简化 从而得到P的幂形式P x 3 2 x 5 计算机显示插值多项式的符号形式 另外一个例子 使用的是本章另一种方法所用的例子 x 1 6 y 161821171512 disp x y u 75 05 6 25 v polyinterp x y u plot x y o u v 其运行后的结果为123456161821171512同时输出3 2 3 2分段线性插值 通过两步操作可以绘制出一个简单的图形 第一步用圆圈在坐标系中标出个数据点plot x y o 第二步用直线段依次连接这些数据点plot x y 下面的语句执行这样的操作 生成图3 3 x 1 6 y 161821171512 plot x y o x y 在生成图3 3所示的图线时 MATLAB图像处理函数使用了分段线性插值 这个分段线性插值算法是其他更复杂算法的基础 他用了三个量 首先要确定间隔序号 intervalindex k 使得第二个量是局部变量 localvariable s 其定义为最后一个量是一次均差 firstdivideddifference 定义了这三个量 则插值基函数可表示为 显然 这是通过点 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 但对于更高次的插值多项式不一定成立 考虑一个定义在区间xk x xk 1的函数 采用局部变量s x xk并令h hk 它可表示为 这是一个关于 也即的三次多项式 它满足四个插值条件 其中两个关于函数值 两个关于函数的导数值 那些满足关于导数值插值条件的函数称为埃米特 hermite 或密切 osculatory 插值基函数 因为这些函数在插值点上保持高阶的连续性 在拉丁文中 密切 一词的本意为 亲吻 如果正好给定了一系列数据点上的函数值和一阶导数值 那么就可以用埃米特插值拟合这些数据 但是如果没有给出这些导数值 那么需要用一些方法来限定斜率dk 我们在下一节中讨论两种可能的办法 即在MATLAB中的函数pchip和spline 3 4保形分段三次插值 pchip实际是 分段三次埃米特插值多项式 piecewisecubicHermiteinterpolatingpolynominal 的英文首字母缩写 有意思的是 根据这个名字并不能确定它到底是哪一种分段三次埃米特插值多项式 因为样条插值函数实际也是分段三次埃米特插值多项式 只是对斜率的限制条件不同而已 在这里 我们说的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的平均值 即3 这个保形插值函数由两个三次多项式组成 它们在中间断点处相接 并且在那一点处的导数都为1 3 如果 k和 k 1和的正负号相同 并且两个子区间长度不等 则dk为加权的调和平均 权重由两个子区间的长度决定 其中 上面介绍了函数pchip中如何确定中间断点处的斜率 而对于整个数据区间的两个端点处的斜率d1和dn 需要用一个稍许不同的 单方向分析的方法加以计算 有关细节请参考文件pchiptx m