MATLAB 三次样条.pdf

第 12 章 三次样条 众所周知 使用高阶多项式的插值常常产生病态的结果 目前 有多种消除病态的方 法 在这些方法中 三次样条是最常用的一种 在 MATLAB 中 实现基本的三次样条插值 的函数有 spline ppval mkpp 和 unmkpp 在这些函数中 仅 spline 在 MATLAB 参考 指南 中有说明 下面几节 将展示在 M 文件函数中实现三次样条的基本特征 12 1 基本特征 在三次样条中 要寻找三次多项式 以逼近每对数据点间的曲线 在样条术语中 这 些数据点称之为断点 因为 两点只能决定一条直线 而在两点间的曲线可用无限多的三 次多项式近似 因此 为使结果具有唯一性 在三次样条中 增加了三次多项式的约束条 件 通过限定每个三次多项式的一阶和二阶导数 使其在断点处相等 就可以较好地确定 所有内部三次多项式 此外 近似多项式通过这些断点的斜率和曲率是连续的 然而 第 一个和最后一个三次多项式在第一个和最后一个断点以外 没有伴随多项式 因此必须通 过其它方法确定其余的约束 最常用的方法 也是函数 spline 所采用的方法 就是采用非 扭结 not a knot 条件 这个条件强迫第一个和第二个三次多项式的三阶导数相等 对最后 一个和倒数第二个三次多项式也做同样地处理 基于上述描述 人们可能猜想到 寻找三次样条多项式需要求解大量的线性方程 实 际上 给定 N 个断点 就要寻找 N 1 个三次多项式 每个多项式有 4 个未知系数 这样 所求解的方程组包含有 4 N 1 个未知数 把每个三次多项式列成特殊形式 并且运用各种 约束 通过求解 N 个具有 N 个未知系数的方程组 就能确定三次多项式 这样 如果有 50 个断点 就有 50 个具有 50 个未知系数的方程组 幸好 用稀疏矩阵 这些方程式能够简 明地列出并求解 这就是函数 spline 所使用的计算未知系数的方法 12 2 分段多项式 在最简单的用法中 spline 获取数据 x 和 y 以及期望值 xi 寻找拟合 x 和 y 的三次样 条内插多项式 然后 计算这些多项式 对每个 xi 的值 寻找相应的 yi 例如 x 0 12 y tan pi x 25 xi linspace 0 12 yi spline x y xi plot x y o xi yi title Spline fit 见图 12 1 样条拟合 这种方法适合于只需要一组内插值的情况 不过 如果需要从相同数据集里获取另一 组内插值 再次计算三次样条系数是没有意义的 在这种情况下 可以调用仅带前两个参 量的 spline 图 12 1 样条拟合 pp spline x y pp Columns 1 through 7 10 0000 1 0000 12 0000 0 1 0000 2 0000 3 0000 Columns 8 through 14 4 0000 5 0000 6 0000 7 0000 8 0000 9 0000 10 0000 Columns 15 through 21 11 0000 12 0000 4 0000 0 0007 0 0007 0 0010 0 0012 Columns 22 through 28 0 0024 0 0019 0 0116 0 0083 0 1068 0 1982 1 4948 Columns 29 through 35 1 4948 0 0001 0 0020 0 0042 0 0072 0 0109 0 0181 Columns 36 through 42 0 0237 0 0586 0 0336 0 3542 0 2406 4 2439 0 1257 Columns 43 through 49 0 1276 0 1339 0 1454 0 1635 0 1925 0 2344 0 3167 Columns 50 through 56 0 4089 0 7967 0 9102 4 9136 0 0 1263 0 2568 Columns 57 through 63 0 3959 0 5498 0 7265 0 9391 1 2088 1 5757 2 1251 Columns 64 through 65 3 0777 5 2422 当采用这种方式调用时 spline 返回一个称之为三次样条的 pp 形式形式或分段多项式形式分段多项式形式 的数组 这个数组包含了对于任意一组所期望的内插值和计算三次样条所必须的全部信息 给定 pp 形式 函数 ppval 计算该三次样条 例如 yi ppval pp xi 计算先前计算过的同样的 yi 类似地 xi2 linspace 10 12 yi2 ppval pp xi2 运用 pp 形式 在限定的更细区间 10 12 内 再次计算该三次样条 xi3 10 15 yi3 ppval pp xi3 yi3 3 0777 5 2422 15 8945 44 0038 98 5389 188 4689 它表明 可在计算三次多项式所覆盖的区间外 计算三次样条 当数据出现在最后一个断 点之后或第一个断点之前时 则分别运用最后一个或第一个三次多项式来寻找内插值 上述给定的三次样条 pp 形式 存储了断点和多项式系数 以及关于三次样条表示的其 它信息 因为 所有信息都被存储在单个向量里 所以这种形式在 MATLAB 中是一种方便 的数据结构 当要计算三次样条表示时 必须把 pp 形式分解成它的各个表示段 在 MATLAB 中 通过函数 unmkpp 完成这一过程 运用上述 pp 形式 该函数给出如下结果 break coefs npolys ncoefs unmkpp pp breaks Columns 1 through 12 0 1 2 3 4 5 6 7 8 9 10 11 Column 13 12 coefs 0 0007 0 0001 0 1257 0 0 0007 0 0020 0 1276 0 1263 0 0010 0 0042 0 1339 0 2568 0 0012 0 0072 0 1454 0 3959 0 0024 0 0109 0 1635 0 5498 0 0019 0 0181 0 1925 0 7265 0 0116 0 0237 0 2344 0 9391 0 0083 0 0586 0 3167 1 2088 0 1068 0 0336 0 4089 1 5757 0 1982 0 3542 0 7967 2 1251 1 4948 0 2406 0 9102 3 0777 1 4948 4 2439 4 9136 5 2422 npolys 12 ncoefs 4 这里 break 是断点 coefs 是矩阵 它的第 i 行是第 i 个三次多项式 npolys 是多项式 的数目 ncoefs 是每个多项式系数的数目 注意 这种形式非常一般 样条多项式不必是 三次 这对于样条的积分和微分是很有益的 给定上述分散形式 函数 mkpp 恢复了 pp 形式 pp mkpp break coefs pp Columns 1 through 7 10 0000 1 0000 12 0000 0 1 0000 2 0000 3 0000 Columns 8 through 14 4 0000 5 0000 6 0000 7 0000 8 0000 9 0000 10 0000 Columns 15 through 21 11 0000 12 0000 4 0000 0 0007 0 0007 0 0010 0 0012 Columns 22 through 28 0 0024 0 0019 0 0116 0 0083 0 1068 0 1982 1 4948 Columns 29 through 35 1 4948 0 0001 0 0020 0 0042 0 0072 0 0109 0 0181 Columns 36 through 42 0 0237 0 0586 0 0336 0 3542 0 2406 4 2439 0 1257 Columns 43 through 49 0 1276 0 1339 0 1454 0 1635 0 1925 0 2344 0 3167 Columns 50 through 56 0 4089 0 7967 0 9102 4 9136 0 0 1263 0 2568 Columns 57 through 63 0 3959 0 5498 0 7265 0 9391 1 2088 1 5757 2 1251 Columns 64 through 65 3 0777 5 2422 因为矩阵 coefs 的大小确定了 npolys 和 neofs 所以 mkpp 不需要 npolys 和 ncoefs 去 重 构 pp 形 式 pp 形 式 的 数 据 结 构 仅 在 mkpp 中 给 定 为 pp 10 1 npolys break ncoefs coefs 前两个元素出现在所有的 pp 形式中 它们作为确认 pp 形式 向量的一种方法 12 3 积分 在大多数情况下 需要知道由三次样条所描述的 自变量为 x 的函数所包含的面积 也就是 如果这个函数记为 y f x 我们感兴趣的是计算 S xs x dx x x 1 其中 是 s x1 0 式中的 x1 是第一个样条的断点 因为 s x 由被连接的三次多项式组成 其中第 k 个三 次多项式为 sk x ak x xk 3 bk x xk 2 ck x xk dk xk x xk 1 并且该函数在区间xk x xk 1所含的面积为 Sxsx dxaxxbxxcxxdxx kk x x kkkkkkkk 1 432 432 三次样条下的面积容易求得为 S xs x dxsx dx i x x i k k x x i i k 1 1 1 式中xk x xk 1 或者 S xS xSx iik i k 1 1 1 式中xk x x 0 1 1 2 pi y sin x create rough data pp spline x y pp form fitting rough data ppi spintgrl pp pp form of integral xi linspace 0 2 pi finer points for interpolation yi ppval pp xi evaluate curve yyi ppval ppi xi evaluate integral plot x y o xi yi xi yyi plot results 注意 如图 12 2 所示 它定性地证明了恒等式 sin cos x dxx x x 1 1 图 12 2 函数的积分图形 12 4 微分 如同人们对样条积分感兴趣一样 一个由样条所描述的函数的微分或斜率也是很有用 的 给定第 k 个三次多项式为 sk x ak x xk 3 bk x xk 2 ck x xk dk xk x x 0 1 1 2 pi same data as earlier y sin x pp spline x y pp form fitting rough data ppd spderiv pp pp form of derivative xi linspace 0 2 pi finer points for interpolation yi ppval pp xi evaluate curve yyd ppval ppd xi evaluate derivative plot x y o xi yi xi yyd plot results 图 12 3 函数的微分图形 注意 图 12 3 定性地证明了恒等式 d dx xxsin cos 12 5 小结 下面的两个表总结了本章所讨论的样条函数 表 12 1 三次样条函数 yi spline x y xi y f x 在 xi 中各点的三次样条插值 pp spline x y 返回 y f x 的分段多项式的表示 yi ppval pp xi 计算 xi 中各点的分段多项式 break coefs npolys ncoefs unmkpp pp 分解分段多项式的表示 pp mkpp break coefs 形成分段多项式的表示 表 12 2 精通 MATLAB 的样条函数 yi spintgrl x y xi 在 xi 各点对 y f x 积分的三次样条的插值 ppi spintgrl