第 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 xx( )( )=∫ 1其中，是 s(x1)=0 式中的 x1 是第一个样条的断点因为 s(x)由被连接的三次多项式组成，其中第 k 个三 次多项式为： sk(x)=ak(x-xk)3+ bk(x-xk)2+ ck(x-xk)+dk, xk>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); % evaluat。