Introduction

Many functions in this chapter produce cubic piecewise polynomial or general spline functions that either interpolate or approximate given data or are support functions for the evaluation and integration of these functions. Three major subdivisions of functions are provided. The cubic spline functions begin with the prefix CS and use the piecewise polynomial representation. The spline functions begin with the prefix BS and use the B-spline representation. The third major subdivision includes functions that operate on the output of both the cubic spline and B-spline functions. Most spline functions are based on routines documented by de Boor (1978).

A univariate piecewise polynomial function, p, is specified by giving its breakpoint sequence

, the order k (degree k – 1) of its polynomial pieces, and the k × (n – 1) matrix C of its local polynomial coefficients. In terms of this information, the piecewise polynomial (ppoly) function is given by the following equation:

The breakpoint sequence ξ is assumed to be strictly increasing, and the ppoly function is extended to the entire real axis by extrapolation from the first and last intervals. This representation is redundant when the ppoly function is known to be smooth. For example, if p is known to be continuous, then c1, i+1 can be computed from the cji as follows:

B-splines provide a particularly convenient and suitable basis for a given class of smooth ppoly functions. Such a class is specified by giving its breakpoint sequence, its order k, and the required smoothness across each of the interior breakpoints. The corresponding B-spline basis is specified by giving its knot sequence

. The specification rule is the following: If the class is to have all derivatives up to and including the jth derivative continuous across the interior breakpoint ξi, then the number ξi should occur k – j – 1 times in the knot sequence. Assuming that ξ1 and ξn are the endpoints of the interval of interest, choose the first k knots equal to ξ1 and the last k knots equal to ξn. This can be done since the B-splines are defined to be right continuous near ξ1 and left continuous near ξn.

When the above construction is completed, a knot sequence t of length M is generated and m: = M – k B-splines of order k (for example, B0, ..., Bm – 1) span the ppoly functions on the interval with the indicated smoothness. That is, each ppoly function in this class has a unique representation:

as a linear combination of B-splines. A B-spline is a particularly compact ppoly function. The function Bi is a nonnegative function that is nonzero only on the interval [ti, ti + k ]. More precisely, the support of the ith B-spline is [ti, ti + k ]. No ppoly function in the same class (other than the zero function) has smaller support (i.e., vanishes on more intervals) than a B-spline. This makes B-splines particularly attractive basis functions since the influence of any particular B-spline coefficient extends only over a few intervals. When it is necessary to emphasize the dependence of the B-spline on its parameters, the notation:

Cubic splines are smooth (i.e., C0, C1 or C2), fourth-order ppoly functions. For historical and other reasons, cubic splines are the most frequently used ppoly functions. Therefore, special functions are provided for their construction and evaluation. These routines use the ppoly representation as described above for general ppoly functions (with k = 4).

Two cubic spline interpolation functions, CSINTERP and CSSHAPE, are provided. Function CSINTERP allows the user to specify various endpoint conditions (such as the value of the first or second derivative at the right and left points). This means that the natural cubic spline can be obtained using this function by setting the second derivative to zero at both endpoints. Function CSSHAPE is designed so that the shape of the curve matches the shape of the data. In particular, one option of this function preserves the convexity of the data while the default attempts to minimize oscillations. The CSINTERP function includes keywords which allow the user to specify tension, continuity and bias parameters at each data point.

It is possible that the cubic spline interpolation functions will produce unsatisfactory results. For example, the interpolant may not have the shape required by the user, or the data may be noisy and require a least-squares fit. The BSINTERP interpolation function is more flexible, as it allows the user to choose the knots and order of the spline interpolant. The user is encouraged to use this routine and exploit the flexibility provided.

The simplest method of obtaining multivariate interpolation and approximation functions is to take univariate methods and form a multivariate method via tensor products. In the case of two-dimensional spline interpolation, the derivation proceeds as follows: Let tx be a knot sequence for splines of order kx and ty be a knot sequence for splines of order ky. Let Nx + kx be the length of tx and Ny + ky be the length of ty. Then, the tensor-product spline has the following form:

for which the corresponding univariate interpolation problem can be solved, the tensor-product interpolation problem finds the coefficients cnm, so that the following is true:

This problem can be solved efficiently by repeatedly solving univariate interpolation problems as described in de Boor (1978, p. 347). Three-dimensional interpolation can be handled in an analogous manner. This chapter provides functions that compute the two-dimensional, tensor-product spline coefficients given two-dimensional interpolation data (BSINTERP) and functions that compute the two-dimensional, tensor-product spline coefficients for a tensor-product, least-squares problem (BSLSQ). In addition, evaluation, differentiation, and integration routines (SPVALUE and SPINTEG) are provided for the two-dimensional, tensor-product spline functions.

PV-WAVE IMSL Mathematics provides functions to interpolate and approximate scattered data in Rn for n ≥ 1. Function SCAT2DINTERP interpolates scattered data in the plane and is based on work by Akima (1978), which uses C1 piecewise quintics on a triangular mesh. Function RADBF can be used to either interpolate or approximate scattered data in Rn for n ≥ 1. The RADBF function computes approximations based on radial-basis functions. The fit computed by RADBF can be evaluated using function RADBE.

PV-WAVE IMSL Mathematics includes functions for smoothing noisy data. Function FCNLSQ computes regressions with user-supplied functions. Function BSLSQ computes a one- or two-dimensional, least-squares fit using splines with fixed knots or variable knots. This function produces cubic-spline, least-squares fit by default. Keywords allow the user to choose the order and the knot sequence. The function RADBF computes an approximation to scattered data in RN using radial-basis functions.

One “smoothing spline” function is provided. The default action of CSSMOOTH estimates a smoothing parameter by cross-validation, then returns the cubic spline that smooths the data. If the user chooses to supply a smoothing parameter, this function returns the appropriate cubic spline.

A spline can be viewed as a mapping with domain Rd and target Rr, where d and r are positive integers. For this version of PV-WAVE IMSL Mathematics, only r = 1 is supported. Thus, if s is a spline, then the following is true for some d and r:

s(x) = (s0(x), ..., sr – 1(x))

All the information for a spline is stored in a structure. Since, in general, the components of this structure are of varying lengths, an anonymous structure is defined for each spline. An example of the information returned by the INFO command with keyword Structures set and an argument containing a spline structure follows:

For ppoly functions, a ppoly is viewed as a mapping with domain Rd and target Rr, where d and r are positive integers. Thus, if p is a ppoly, then the following is true for some d and r:

For this version of PV-WAVE IMSL Mathematics, only r = 1 is supported. This implies that such a ppoly p must have d breakpoint sequences and orders (one for each domain dimension). Thus, associated with p, breakpoints and orders are as follows:

p(x) = (p0(x), ..., pr(x))