第12章三次样条

第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(‘Splinefit‘)（见图12.1样条拟合）这种方法适合于只需要一组内插值的情况。不过，如果需要从相同数据集里获取另一组内插值，再次计算三次样条系数是没有意义的。在这种情况下，可以调用仅带前两个参量的spline：图12.1样条拟合pp=spline(x,y)pp=Columns1through710.00001.000012.000001.00002.00003.0000Columns8through144.00005.00006.00007.00008.00009.000010.0000Columns15through2111.000012.00004.00000.00070.00070.00100.0012Columns22through280.00240.00190.0116-0.00830.1068-0.19821.4948Columns29through351.4948-0.00010.00200.00420.00720.01090.0181Columns36through420.02370.05860.03360.3542-0.24064.24390.1257Columns43through490.12760.13390.14540.16350.19250.23440.3167Columns50through560.40890.79670.91024.913600.12630.2568Columns57through630.39590.54980.72650.93911.20881.57572.1251Columns64through653.07775.2422当采用这种方式调用时，spline返回一个称之为三次样条的pp形式或分段多项式形式的数组。这个数组包含了对于任意一组所期望的内插值和计算三次样条所必须的全部信息。给定pp形式，函数ppval计算该三次样条。例如，yi=ppval(pp,xi);计算先前计算过的同样的yi。类似地，xi2=linspace(10,12);yi2=ppval(pp,xi2);运用pp形式，在限定的更细区间[10，12]内，再次计算该三次样条。xi3=10:15yi3=ppval(pp,xi3)yi3=3.07775.242215.894544.003898.5389188.4689它表明，可在计算三次多项式所覆盖的区间外，计算三次样条。当数据出现在最后一个断点之后或第一个断点之前时，则分别运用最后一个或第一个三次多项式来寻找内插值。上述给定的三次样条pp形式，存储了断点和多项式系数，以及关于三次样条表示的其它信息。因为，所有信息都被存储在单个向量里，所以这种形式在MATLAB中是一种方便的数据结构。当要计算三次样条表示时，必须把pp形式分解成它的各个表示段。在MATLAB中，通过函数unmkpp完成这一过程。运用上述pp形式，该函数给出如下结果：[break,coefs,npolys,ncoefs]=unmkpp(pp)breaks=Columns1through1201234567891011Column1312coefs=0.0007-0.00010.125700.00070.00200.12760.12630.00100.00420.13390.25680.00120.00720.14540.39590.00240.01090.16350.54980.00190.01810.19250.72650.01160.02370.23440.9391-0.00830.05860.31671.20880.10680.03360.40891.5757-0.19820.35420.79672.12511.4948-0.24060.91023.07771.49484.24394.91365.2422npolys=12ncoefs=4这里break是断点，coefs是矩阵，它的第i行是第i个三次多项式，npolys是多项式的数目，ncoefs是每个多项式系数的数目。注意，这种形式非常一般，样条多项式不必是三次。这对于样条的积分和微分是很有益的。给定上述分散形式，函数mkpp恢复了pp形式。pp=mkpp(break,coefs)pp=Columns1through710.00001.000012.000001.00002.00003.0000Columns8through144.00005.00006.00007.00008.00009.000010.0000Columns15through2111.000012.00004.00000.00070.00070.00100.0012Columns22through280.00240.00190.0116-0.00830.1068-0.19821.4948Columns29through351.4948-0.00010.00200.00420.00720.01090.0181Columns36through420.02370.05860.03360.3542-0.24064.24390.1257Columns43through490.12760.13390.14540.16350.19250.23440.3167Columns50through560.40890.79670.91024.913600.12630.2568Columns57through630.39590.54980.72650.93911.20881.57572.1251Columns64through653.07775.2422因为矩阵coefs的大小确定了npolys和neofs，所以mkpp不需要npolys和ncoefs去重构pp形式。pp形式的数据结构仅在mkpp中给定为pp=[101npolysbreak(:)‘ncoefscoefs(:)‘]。前两个元素出现在所有的pp形式中，它们作为确认pp形式向量的一种方法。12.3积分在大多数情况下，需要知道由三次样条所描述的、自变量为x的函数所包含的面积。也就是，如果这个函数记为y=f(x)，我们感兴趣的是计算：Sxsxdxxx()()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所含的面积为：Sxsxdxaxxbxxcxxdxxkkxxkkkkkkkk()()/()/()/()()1432432三次样条下的面积容易求得为：Sxsxdxsxdxixxikkxxiik()()()111式中xk=x=xk+1或者SxSxSxiikik()()()111式中xk=x=xk+1上式相加项是所有处理的三次多项式面积的累加和。照此，因为Sk(x)是一个多项式，所以该面积很容易计算，且在描述S(x)的多项式中可形成常数项。有了上述理解，积分本身可写成样条形式。在这种情况下，因为单个多项式具有4阶，所以积分是四次样条。MATLAB中使用的pp形式支持任意阶的样条，所以在《精通MATLAB工具箱》中的函数spintgrl体现了上述样条积分。该函数的主体如下：functionz=spintgrl(x,y,xi)%SPINTGRLCubicSplineIntegralInterpolation%YI=SPINTRGL(X,Y,XI)usescubicsplineinterpolationtofitthedateinXandY,integrates%thesplineandreturnsvaluesoftheintegralevaluatedatthepointsinXI.%%PPI=SPINTGRL(PP)returnsthepiecewisepolynomialvectorPPI%describingtheintegralofthecubicsplinedescribedby%thepiecewisepolynomialinPP.PPisreturnedbythefunction%SPLINEandisadatavectorcontainingallinformationto%evaluateandmanipulateaspline.%%YI=SPINTGRL(PP,XI)integratesthecubicsplinesgivenby%thepiecewisepolynomialPP,andreturnsthevaluesofthe%integralevaluatedatthepointsinXI.%%SeealsoSPLINE,PPVAL,MKPP,UNMKPP,SPDERIV%Copyright(c)1996byPrentice-Hall,Inc.ifnargin==3pp=spline(x,y)elsepp=x;end[br,co,npy,nco]=unmkpp(pp);%takeapartppifpp(1)==10error(‘Splinedatadoesnothavethecorrectform.‘)endsf=nco:-1:1;%scalefactorsforintegrationico=[co./sf(ones(npy,1),:)zeros(npy,1)];%integralcoefficientsnco=nco+1;%splineorderincreasesfork=2:npy%findconstanttermsico(k,nco)=polyval(ico(k-1,:),b