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