第12章三次样条 众所周知,使用高阶多项式的插值常常产生病态的结果。目前,有多种消除病态的方 法。在这些方法中,三次样条是最常用的一种。在 MATLAB中,实现基本的三次样条插值 的函数有 spline,ppva,mkpp和 unmkpp。在这些函数中,仅 spline在《 MATLAB参考 指南》中有说明。下面几节,将展示在M文件函数中实现三次样条的基本特征 121基本特征 在三次样条中,要寻找三次多项式,以逼近每对数据点间的曲线。在样条术语中,这 些数据点称之为断点。因为,两点只能决定一条直线,而在两点间的曲线可用无限多的三次 多项式近似。因此,为使结果具有唯一性。在三次样条中,增加了三次多项式的约束条件。 通过限定每个三次多项式的一阶和二阶导数,使其在断点处相等,就可以较好地确定所有内 部三次多项式。此外,近似多项式通过这些断点的斜率和曲率是连续的。然而,第一个和最 后一个三次多项式在第一个和最后一个断点以外,没有伴随多项式。因此必须通过其它方法 确定其余的约束。最常用的方法,也是函数 spline所采用的方法,就是采用非扭结(not-akno 条件。这个条件强迫第一个和第二个三次多项式的三阶导数相等。对最后一个和倒数第二个 三次多项式也做同样地处理。 基于上述描述,人们可能猜想到,寻找三次样条多项式需要求解大量的线性方程。实 际上,给定N个断点,就要寻找N-1个三次多项式,每个多项式有4个未知系数。这样, 所求解的方程组包含有4*(N-1)个未知数。把每个三次多项式列成特殊形式,并且运用各种 约束,通过求解N个具有N个未知系数的方程组,就能确定三次多项式。这样,如果有50 个断点,就有50个具有50个未知系数的方程组。幸好,用稀疏矩阵,这些方程式能够简明 地列出并求解,这就是函数 spline所使用的计算未知系数的方法。 122分段多项式 在最简单的用法中, spline获取数据x和y以及期望值xi,寻找拟合x和y的三次样条 内插多项式,然后,计算这些多项式,对每个xi的值,寻找相应的yi。例如: >>plot(x, y, 1), title(" Spline fit) (见图12.1样条拟合)
第 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 样条拟合）
这种方法适合于只需要一组内插值的情况。不过,如果需要从相同数据集里获取另一 组内插值,再次计算三次样条系数是没有意义的。在这种情况下,可以调用仅带前两个参量 的spli Spline fit 16 8 4 图12.1样条拟合 pp Columns 1 through 7 1000001.000012.0000 0100002.00003.0000 Columns 8 through 14 4.00005.00006.0000700008.00009.000010.0000 Columns 15 through 21 11.000012.00004.00000.00070.00070.00100.0012 Columns 22 through 28 0.00240.00190.0116-0.00830.1068-0.19821.4948 Columns 29 through 35 149480.00010.00200.00420.00720.01090.0181 Columns 36 through 42 0.02370.05860.03360.3542-0.24064.24390.1257 Columns 43 through 49 0.12760.13390.14540.16350.19250.23440.3167 Columns 50 through 56 0.40890.79670.91024.9136 00.12630.2568 Columns 57 through 63 0.39590.5498.72650.93911.2088 57572.1251 Columns 64 through 65
这种方法适合于只需要一组内插值的情况。不过，如果需要从相同数据集里获取另一 组内插值，再次计算三次样条系数是没有意义的。在这种情况下，可以调用仅带前两个参量 的 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.07775.2422 当采用这种方式调用时, spline返回一个称之为三次样条的pp形式或分段多项式形式 的数组。这个数组包含了对于任意一组所期望的内插值和计算三次样条所必须的全部信息。 给定pp形式,函数 ppva计算该三次样条。例如, >>yi=ppva(pp, xi) 计算先前计算过的同样的yi 类似地 >>X12= linspace(10,12) >>y12=ppva(pp, x12) 运用pp形式,在限定的更细区间[10,12]内,再次计算该三次样条 >>xi3=10:15 >>yi3=ppva(pp, x13 3.07775.242215.894544.003898.53891884689 它表明,可在计算三次多项式所覆盖的区间外,计算三次样条。当数据出现在最后一个断点 之后或第一个断点之前时,则分别运用最后一个或第一个三次多项式来寻找内插值 上述给定的三次样条pp形式,存储了断点和多项式系数,以及关于三次样条表示的其 它信息。因为,所有信息都被存储在单个向量里,所以这种形式在 MATLAB中是一种方便 的数据结构。当要计算三次样条表示时,必须把p形式分解成它的各个表示段在 MATLAB 中,通过函数 unmkpp完成这一过程。运用上述pp形式,该函数给出如下结果 break, coefs, polys, ncoefs]=unmkpp(pp) breaks= Columns I through 12 Column 13 coefs 0.0007-0.00010.1257 0 0.00070.00200.12760.1263 0.00100.00420.13390.2568 0.00120.00720.14540.3959 0.00240.01090.16350.5498 0.00190.01810.19250.7265
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.01160.02370.23440.9391 -0.00830.05860.31671.2088 0.10680.03360.40891.5757 0.19820.35420.79672.1251 149480.24060.91023.0777 14948424394.91365.2422 ncoefs= 这里 break是断点,coes是矩阵,它的第i行是第i个三次多项式, polys是多项式 的数目, ncoefs是每个多项式系数的数目。注意,这种形式非常一般,样条多项式不必是三 次。这对于样条的积分和微分是很有益的 给定上述分散形式,函数mkpp恢复了pp形式。 >>pp=mkpp(break, coefs) Columns 1 through 7 1.000012.0000 01.00002.00003.0000 Columns 8 through 14 4.00005.00006.00007.00008.00009.000010.0000 Columns 15 through 21 11.00001200004.00000.00070.00070.00100.0012 Columns 22 through 28 0.00240.00190.0116-0.00830.10680.19821.4948 Columns 29 through 35 1.4948-0.00010.00200.00420.00720.010900181 Columns 36 through 42 0.02370.05860.03360.3542-0.24064.24390.1257 Columns 43 through 49 0.12760.13390.14540.16350.19250.23440.3167 Columns 50 through 56 0.40890.79670.91024.9136 00.12630.2568 Columns 57 through 63 0.39590.54980.72650.93911.20881.575721251 Columns 64 through 65 3.07775.2422 因为矩阵 coefs的大小确定了 polys和 neos,所以mkpp不需要 polys和 ncoefs 重构pp形式。pp形式的数据结构仅在mkpp中给定为pp={101 polys break() ncoefs
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(()l前两个元素出现在所有的pp形式中,它们作为确认pp形式向量的一种方法 123积分 在大多数情况下,需要知道由三次样条所描述的、自变量为x的函数所包含的面积。 也就是,如果这个函数记为y=f(x),我们感兴趣的是计算: 其中,是sx 式中的x1是第一个样条的断点。因为sx)由被连接的三次多项式组成,其中第k个三 次多项式为: Sk(x=ak(x-xk+ bx(x-xxF+ck(x-xkH+dk, XK<=X<=Xk+ 并且该函数在区间xk<=x<=xk+1所含的面积为: ()=s(x)x=a24x-x)+b/Xx-x)2+c/2(x-x1)2+d(x 三次样条下的面积容易求得为 S(x)=∑(xx+s(x)dx式中x<X=x+ 或者 S(x)=∑S(x1+1)+S(x) 式中xk<=x<=xk+1 上式相加项是所有处理的三次多项式面积的累加和。照此,因为Sk(x)是一个多项式 所以该面积很容易计算,且在描述S(x)的多项式中可形成常数项。有了上述理解,积分本身 可写成样条形式。在这种情况下,因为单个多项式具有4阶,所以积分是四次样条 MATLAB中使用的pp形式支持任意阶的样条,所以在《精通 MATLAB工具箱》中的 函数 spintgrl体现了上述样条积分。该函数的主体如下: function z=spintgrl(x, y, xi) SPINTGRL Cubic Spline Integral Interpolation YI=SPINTRGL(X, Y, XI) uses cubic spline interpolation to fit the date in X and Y Integrates the spline and returns values of the integral evaluated at the points in XI
coefs(:)‘]。前两个元素出现在所有的 pp 形式中，它们作为确认 pp 形式向量的一种方法。 12.3 积分 在大多数情况下，需要知道由三次样条所描述的、自变量为 x 的函数所包含的面积。 也就是，如果这个函数记为 y=f(x)，我们感兴趣的是计算： S x s 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 所含的面积为： Sk x sk x dx a x x b x x c x x d x x x x k k k k k k k k ( ) = ( ) = / ( − ) + / ( − ) + / ( − ) + ( − ) 1 4 3 2 4 3 2 三次样条下的面积容易求得为： S x si x dx s x dx x x i k k x x i i k ( ) = ( ) + ( ) = = − 1 1 1 式中 xk<=x<= xk+1 或者 S x Si xi Sk x i k ( ) = ( + ) + ( ) = − 1 1 1 式中 xk<=x<= xk+1 上式相加项是所有处理的三次多项式面积的累加和。照此，因为 Sk(x)是一个多项式， 所以该面积很容易计算，且在描述 S(x)的多项式中可形成常数项。有了上述理解，积分本身 可写成样条形式。在这种情况下，因为单个多项式具有 4 阶，所以积分是四次样条。 MATLAB 中使用的 pp 形式支持任意阶的样条，所以在《精通 MATLAB 工具箱》中的 函数 spintgrl 体现了上述样条积分。该函数的主体如下： function z=spintgrl(x, y, xi) % SPINTGRL Cubic Spline Integral Interpolation % YI=SPINTRGL(X, Y, XI) uses cubic spline interpolation to fit the date in X and Y, integrates % the spline and returns values of the integral evaluated at the points in XI