PDA_BSPDOC

Documentation for working with piecewise polynomial functions in B-representation.

Origin

SLATEC / CAMSUN
        SUBROUTINE PDA_BSPDOC
  
  
  ***BEGIN PROLOGUE  PDA_BSPDOC
  ***PURPOSE  Documentation for BSPLINE, a package of subprograms for
              working with piecewise polynomial functions
              in B-representation.
  ***LIBRARY   SLATEC
  ***CATEGORY  E, E1A, K, Z
  ***TYPE      ALL (PDA_BSPDOC-A)
  ***KEYWORDS  B-SPLINE, DOCUMENTATION, SPLINES
  ***AUTHOR  Amos, D. E., (SNLA)
  ***DESCRIPTION
  
       Abstract
           PDA_BSPDOC is a non-executable, B-spline documentary routine.
           The narrative describes a B-spline and the routines
           necessary to manipulate B-splines at a fairly high level.
           The basic package described herein is that of reference
           5 with names altered to prevent duplication and conflicts
           with routines from reference 3.  The call lists used here
           are also different.  Work vectors were added to ensure
           portability and proper execution in an overlay environ-
           ment.  These work arrays can be used for other purposes
           except as noted in BSPVN.  While most of the original
           routines in reference 5 were restricted to orders 20
           or less, this restriction was removed from all routines
           except the quadrature routine BSQAD.  (See the section
           below on differentiation and integration for details.)
  
           The subroutines referenced below are single precision
           routines.  Corresponding double precision versions are also
           part of the package, and these are referenced by prefixing
           a D in front of the single precision name.  For example,
           BVALU and PDA_DBVALU are the single and double precision
           versions for evaluating a B-spline or any of its deriva-
           tives in the B-representation.
  
                  ****Description of B-Splines****
  
       A collection of polynomials of fixed degree K-1 defined on a
       subdivision (X(I),X(I+1)), I=1,...,M-1 of (A,B) with X(1)=A,
       X(M)=B is called a B-spline of order K.  If the spline has K-2
       continuous derivatives on (A,B), then the B-spline is simply
       called a spline of order K.  Each of the M-1 polynomial pieces
       has K coefficients, making a total of K(M-1) parameters.  This
       B-spline and its derivatives have M-2 jumps at the subdivision
       points X(I), I=2,...,M-1.  Continuity requirements at these
       subdivision points add constraints and reduce the number of free
       parameters.  If a B-spline is continuous at each of the M-2 sub-
       division points, there are K(M-1)-(M-2) free parameters; if in
       addition the B-spline has continuous first derivatives, there
       are K(M-1)-2(M-2) free parameters, etc., until we get to a
       spline where we have K(M-1)-(K-1)(M-2) = M+K-2 free parameters.
       Thus, the principle is that increasing the continuity of
       derivatives decreases the number of free parameters and
       conversely.
  
       The points at which the polynomials are tied together by the
       continuity conditions are called knots.  If two knots are
       allowed to come together at some X(I), then we say that we
       have a knot of multiplicity 2 there, and the knot values are
       the X(I) value.  If we reverse the procedure of the first
       paragraph, we find that adding a knot to increase multiplicity
       increases the number of free parameters and, according to the
       principle above, we thereby introduce a discontinuity in what
       was the highest continuous derivative at that knot.  Thus, the
       number of free parameters is N = NU+K-2 where NU is the sum
       of multiplicities at the X(I) values with X(1) and X(M) of
       multiplicity 1 (NU = M if all knots are simple, i.e., for a
       spline, all knots have multiplicity 1.)  Each knot can have a
       multiplicity of at most K.  A B-spline is commonly written in the
       B-representation
  
                 Y(X) = sum( A(I)*B(I,X), I=1 , N)
  
       to show the explicit dependence of the spline on the free
       parameters or coefficients A(I)=BCOEF(I) and basis functions
       B(I,X).  These basis functions are themselves special B-splines
       which are zero except on (at most) K adjoining intervals where
       each B(I,X) is positive and, in most cases, hat or bell-
       shaped.  In order for the nonzero part of B(1,X) to be a spline
       covering (X(1),X(2)), it is necessary to put K-1 knots to the
       left of A and similarly for B(N,X) to the right of B.  Thus, the
       total number of knots for this representation is NU+2K-2 = N+K.
       These knots are carried in an array T(*) dimensioned by at least
       N+K.  From the construction, A=T(K) and B=T(N+1) and the spline is
       defined on T(K).LE.X.LE.T(N+1).  The nonzero part of each basis
       function lies in the  Interval (T(I),T(I+K)).  In many problems
       where extrapolation beyond A or B is not anticipated, it is common
       practice to set T(1)=T(2)=...=T(K)=A and T(N+1)=T(N+2)=...=
       T(N+K)=B.  In summary, since T(K) and T(N+1) as well as
       interior knots can have multiplicity K, the number of free
       parameters N = sum of multiplicities - K.  The fact that each
       B(I,X) function is nonzero over at most K intervals means that
       for a given X value, there are at most K nonzero terms of the
       sum.  This leads to banded matrices in linear algebra problems,
       and references 3 and 6 take advantage of this in con-
       structing higher level routines to achieve speed and avoid
       ill-conditioning.
  
                       ****Basic Routines****
  
       The basic routines which most casual users will need are those
       concerned with direct evaluation of splines or B-splines.
       Since the B-representation, denoted by (T,BCOEF,N,K), is
       preferred because of numerical stability, the knots T(*), the
       B-spline coefficients BCOEF(*), the number of coefficients N,
       and the order K of the polynomial pieces (of degree K-1) are
       usually given.  While the knot array runs from T(1) to T(N+K),
       the B-spline is normally defined on the interval T(K).LE.X.LE.
       T(N+1).  To evaluate the B-spline or any of its derivatives
       on this interval, one can use
  
                    Y = BVALU(T,BCOEF,N,K,ID,X,INBV,WORK)
  
       where ID is an integer for the ID-th derivative, 0.LE.ID.LE.K-1.
       ID=0 gives the zero-th derivative or B-spline value at X.
       If X.LT.T(K) or X.GT.T(N+1), whether by mistake or the result
       of round off accumulation in incrementing X, BVALU gives a
       diagnostic.  INBV is an initialization parameter which is set
       to 1 on the first call.  Distinct splines require distinct
       INBV parameters.  WORK is a scratch vector of length at least
       3*K.
  
       When more conventional communication is needed for publication,
       physical interpretation, etc., the B-spline coefficients can
       be converted to piecewise polynomial (PP) coefficients.  Thus,
       the breakpoints (distinct knots) XI(*), the number of
       polynomial pieces LXI, and the (right) derivatives C(*,J) at
       each breakpoint XI(J) are needed to define the Taylor
       expansion to the right of XI(J) on each interval XI(J).LE.
       X.LT.XI(J+1), J=1,LXI where XI(1)=A and XI(LXI+1)=B.
       These are obtained from the (T,BCOEF,N,K) representation by
  
                  CALL BSPPP(T,BCOEF,N,K,LDC,C,XI,LXI,WORK)
  
       where LDC.GE.K is the leading dimension of the matrix C and
       WORK is a scratch vector of length at least K*(N+3).
       Then the PP-representation (C,XI,LXI,K) of Y(X), denoted
       by Y(J,X) on each interval XI(J).LE.X.LT.XI(J+1), is
  
       Y(J,X) = sum( C(I,J)*((X-XI(J))**(I-1))/factorial(I-1), I=1,K)
  
       for J=1,...,LXI.  One must view this conversion from the B-
       to the PP-representation with some skepticism because the
       conversion may lose significant digits when the B-spline
       varies in an almost discontinuous fashion.  To evaluate
       the B-spline or any of its derivatives using the PP-
       representation, one uses
  
                  Y = PPVAL(LDC,C,XI,LXI,K,ID,X,INPPV)
  
       where ID and INPPV have the same meaning and usage as ID and
       INBV in BVALU.
  
       To determine to what extent the conversion process loses
       digits, compute the relative error ABS((Y1-Y2)/Y2) over
       the X interval with Y1 from PPVAL and Y2 from BVALU.  A
       major reason for considering PPVAL is that evaluation is
       much faster than that from BVALU.
  
       Recall that when multiple knots are encountered, jump type
       discontinuities in the B-spline or its derivatives occur
       at these knots, and we need to know that BVALU and PPVAL
       return right limiting values at these knots except at
       X=B where left limiting values are returned.  These values
       are used for the Taylor expansions about left end points of
       breakpoint intervals.  That is, the derivatives C(*,J) are
       right derivatives.  Note also that a computed X value which,
       mathematically, would be a knot value may differ from the knot
       by a round off error.  When this happens in evaluating a dis-
       continuous B-spline or some discontinuous derivative, the
       value at the knot and the value at X can be radically
       different.  In this case, setting X to a T or XI value makes
       the computation precise.  For left limiting values at knots
       other than X=B, see the prologues to BVALU and other
       routines.
  
                       ****Interpolation****
  
       BINTK is used to generate B-spline parameters (T,BCOEF,N,K)
       which will interpolate the data by calls to BVALU.  A similar
       interpolation can also be done for cubic splines using BINT4
       or the code in reference 7.  If the PP-representation is given,
       one can evaluate this representation at an appropriate number of
       abscissas to create data then use BINTK or BINT4 to generate
       the B-representation.
  
                 ****Differentiation and Integration****
  
       Derivatives of B-splines are obtained from BVALU or PPVAL.
       Integrals are obtained from BSQAD using the B-representation
       (T,BCOEF,N,K) and PPQAD using the PP-representation (C,XI,LXI,
       K).  More complicated integrals involving the product of a
       of a function F and some derivative of a B-spline can be
       evaluated with BFQAD or PFQAD using the B- or PP- represen-
       tations respectively.  All quadrature routines, except for PPQAD,
       are limited in accuracy to 18 digits or working precision,
       whichever is smaller.  PPQAD is limited to working precision
       only.  In addition, the order K for BSQAD is limited to 20 or
       less.  If orders greater than 20 are required, use BFQAD with
       F(X) = 1.
  
                        ****Extrapolation****
  
       Extrapolation outside the interval (A,B) can be accomplished
       easily by the PP-representation using PPVAL.  However,
       caution should be exercised, especially when several knots
       are located at A or B or when the extrapolation is carried
       significantly beyond A or B.  On the other hand, direct
       evaluation with BVALU outside A=T(K).LE.X.LE.T(N+1)=B
       produces an error message, and some manipulation of the knots
       and coefficients are needed to extrapolate with BVALU.  This
       process is described in reference 6.
  
                  ****Curve Fitting and Smoothing****
  
       Unless one has many accurate data points, direct inter-
       polation is not recommended for summarizing data.  The
       results are often not in accordance with intuition since the
       fitted curve tends to oscillate through the set of points.
       Monotone splines (reference 7) can help curb this undulating
       tendency but constrained least squares is more likely to give an
       acceptable fit with fewer parameters.  Subroutine FC, des-
       cribed in reference 6, is recommended for this purpose.  The
       output from this fitting process is the B-representation.
  
                **** Routines in the B-Spline Package ****
  
                        Single Precision Routines
  
           The subroutines referenced below are SINGLE PRECISION
           routines. Corresponding DOUBLE PRECISION versions are also
           part of the package and these are referenced by prefixing
           a D in front of the single precision name. For example,
           BVALU and PDA_DBVALU are the SINGLE and DOUBLE PRECISION
           versions for evaluating a B-spline or any of its deriva-
           tives in the B-representation.
  
       BINT4 - interpolates with splines of order 4
       BINTK - interpolates with splines of order k
       BSQAD - integrates the B-representation on subintervals
       PPQAD - integrates the PP-representation
       BFQAD - integrates the product of a function F and any spline
               derivative in the B-representation
       PFQAD - integrates the product of a function F and any spline
               derivative in the PP-representation
       BVALU - evaluates the B-representation or a derivative
       PPVAL - evaluates the PP-representation or a derivative
       INTRV - gets the largest index of the knot to the left of x
       BSPPP - converts from B- to PP-representation
       BSPVD - computes nonzero basis functions and derivatives at x
       BSPDR - sets up difference array for BSPEV
       BSPEV - evaluates the B-representation and derivatives
       BSPVN - called by BSPEV, BSPVD, BSPPP and BINTK for function and
               derivative evaluations
                          Auxiliary Routines
  
         BSGQ8,PPGQ8,BNSLV,BNFAC,PDA_XERMSG,DBSGQ8,DPPGQ8,PDA_DBNSLV,PDA_DBNFAC
  
                      Machine Dependent Routines
  
                        PDA_I1MACH, R1MACH, PDA_D1MACH
  
  ***REFERENCES  1. D. E. Amos, Computation with splines and
                   B-splines, Report SAND78-1968, Sandia
                   Laboratories, March 1979.
                 2. D. E. Amos, Quadrature subroutines for splines and
                   B-splines, Report SAND79-1825, Sandia Laboratories,
                   December 1979.
                 3. Carl de Boor, A Practical Guide to Splines, Applied
                   Mathematics Series 27, Springer-Verlag, New York,
                   1978.
                 4. Carl de Boor, On calculating with B-Splines, Journal
                   of Approximation Theory 6, (1972), pp. 50-62.
                 5. Carl de Boor, Package for calculating with B-splines,
                   SIAM Journal on Numerical Analysis 14, 3 (June 1977),
                   pp. 441-472.
                 6. R. J. Hanson, Constrained least squares curve fitting
                   to discrete data using B-splines, a users guide,
                   Report SAND78-1291, Sandia Laboratories, December
                   1978.
                 7. F. N. Fritsch and R. E. Carlson, Monotone piecewise
                   cubic interpolation, SIAM Journal on Numerical Ana-
                   lysis 17, 2 (April 1980), pp. 238-246.
  ***ROUTINES CALLED  (NONE)
  ***REVISION HISTORY  (YYMMDD)
     810223  DATE WRITTEN
     861211  REVISION DATE from Version 3.2
     891214  Prologue converted to Version 4.0 format.  (BAB)
     900723  PURPOSE section revised.  (WRB)
     920501  Reformatted the REFERENCES section.  (WRB)
  ***END PROLOGUE  PDA_BSPDOC