精华区文章阅读

发信人: bage (最近比较烦), 信区: Matlab
标  题: 精通Matlab(八)
发信站: 哈工大紫丁香 (Sun Feb  4 14:03:58 2001), 转信

发信人: Security (淼水), 信区: MathTools
发信站: BBS 水木清华站 (Tue Jun  1 22:37:21 1999)

第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(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所含的面积为：

三次样条下的面积容易求得为：

    式中xk<=x<= xk+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.
%
%  PPI=SPINTGRL(PP) returns the piecewise polynomial vector PPI
%  describing the integral of the cubic spline described by
%  the piecewise polynomial in PP. PP is returned by the function
%  SPLINE and is a data vector containing all information to
%  evaluate and manipulate a spline.
%
%  YI=SPINTGRL(PP, XI) integrates the cubic splines given by
%  the piecewise polynomial PP, and returns the values of the
%  integral evaluated at the points in XI.
%
%  See also SPLINE, PPVAL, MKPP, UNMKPP, SPDERIV

%  Copyright  (c) 1996 by Prentice-Hall, Inc.
if nargin==3
pp=spline(x, y)
else
        pp=x;
end
[br, co, npy, nco]=unmkpp(pp);  %  take apart pp
if pp(1)==10
        error(' Spline data does not have the correct form. ')
end
sf=nco : -1 : 1;  %  scale factors for integration
ico=[co ./ sf(ones(npy, 1), : ) zeros(npy, 1)];  %  integral coefficients
nco=nco+1;  %  spline order increases
for k=2 : npy  %  find constant terms
        ico(k, nco)=polyval(ico(k-1, : ),br(k)-br(k-1));
end
ppi=mkpp(br, ico);  %  build pp form for integral
if nargin==1
        z=ppi;
elseif nargin==2
        z=ppval(ppi, y);
else
        z=ppval(ppi, xi);
end

考虑如下运用spintgrl的例子：

>>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所示，它定性地证明了恒等式：

图12.2  函数的积分图形

12.4  微分

如同人们对样条积分感兴趣一样，一个由样条所描述的函数的微分或斜率也是很有用的。给定第k个三次多项式为：

sk(x)=ak(x-xk)3+ bk(x-xk)2+ ck(x-xk)+dk,  xk<=x<= xk+1

容易写出其导数为：

如同样条积分，样条微分也是一种样条。不过，在这种情况下，这个多项式为二阶，所以它是二次样条。
基于上述描述，《精通MATLAB工具箱》中的函数spderiv实施样条微分。spderiv的主体为：

function z=spderiv(x, y, xi)
%  SPDERIV Cubic Spline Derivative Interpolation
%  YI=SPDERIV(X, Y, XI) uses cubic spline interpolation to fit the
%  data in X and Y, differentiates the spline and returns values
%  of the spline derivatives evaluated at the points in XI.
%
%  PPD=SPDERIV(PP) returns the piecewise polynomial vector PPD
%  describing the cubic spline derivative of the curve described by
%  the piecewise polynomial in PP. PP is returned by the function
%  SPLINE and is a data vector containing all information to
%  evaluate and manipulate a spline.
%
%  YI=SPDERIV(PP, XI) differentiates the cubic spline given by
%  the piecewise polynomial PP, and returns the value of the
%  spline derivatives evaluated at the points in XI.
%
%  See also SPLINE, PPVAL, MKPP, UNMKPP, SPINTGRL

%  Copyright (c) 1996 by Prentice-Hall, Inc.

if nargin==3
        pp=spline(x, y);
else
        pp=x;
end
[br, co, npy, nco]=unmkpp(pp);  %  take apart pp
if nco==1 | pp(1)~=10
        error(' Spline data does not have the correct PP form. ')
end
sf=nco-1 : -1 : 1;  %  scale factors for differentiation
dco=sf(ones(npy, 1),:).*co( : , 1 : nco-1);  %  derivative coefficients
ppd=mkpp(br, dco);  %  build pp form for derivative
if nargin==1
        z=ppd;
elseif nargin==2
        z=ppval(ppd, y);
else
        z=ppval(ppd, xi);
end

为演示spderiv的使用，考虑如下例子：

>>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定性地证明了恒等式：

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(pp)        给定y=f(x)的分段多项式的表示，返回y=f(x)积分的分段多项式的表示。
yi=spintgrl(pp,xi)      给定y=f(x)的分段多项式的表示，计算点xi的值，寻找 y=f(x)积分的分段多项式的描述。
yi=spderiv(x,y,xi)      在xi各点对y=f(x)微分的三次样条的插值
ppi=spderiv(pp) 给定y=f(x)的分段多项式的表示，返回y=f(x)微分的分段多项式的表示。
yi=spderiv(pp,xi)       给定y=f(x)的分段多项式的表示，计算点xi的值，寻找y=f(x)微分的分段多项式的表示。

--

--
☆ 来源:．哈工大紫丁香 bbs.hit.edu.cn．[FROM: bage.bbs@smth.org]

Matlab 版 (精华区)