Computes a shape-preserving cubic spline.

The type Imsl_d_ppoly function is imsl_d_cub_spline_interp_shape.

A pointer to the structure that represents the cubic spline interpolant. If an interpolant cannot be computed, then NULL is returned. To release this space, use imsl_free.

The function imsl_f_cub_spline_interp_shape computes a C1 cubic spline interpolant to a set of data points(xi, fi) for i = 0, …, ndata − 1 = n. The breakpoints of the spline are the abscissas. This computation is based on a method by Akima (1970) to combat wiggles in the interpolant. Endpoint conditions are automatically determined by the program; see Akima (1970) or de Boor (1978).

If the optional argument IMSL_CONCAVE is chosen, then this function computes a cubic spline interpolant to the data. For ease of explanation, we will assume that xi < xi+1, although it is not necessary for the user to sort these data values. If the data are strictly convex, then the computed spline is convex, C2, and minimizes the expression

over all convex C1 functions that interpolate the data. In the general case, when the data have both convex and concave regions, the convexity of the spline is consistent with the data, and the above integral is minimized under the appropriate constraints. For more information on this interpolation scheme, refer to Michelli et al. (1985) and Irvine et al. (1986).

One important feature of the splines produced by this function is that it is not possible, a priori, to predict the number of breakpoints of the resulting interpolant. In most cases, there will be breakpoints at places other than data locations. This function should be used when it is important to preserve the convex and concave regions implied by the data.

Both methods are nonlinear, and although the interpolant is a piecewise cubic, cubic polynomials are not reproduced. (However, linear polynomials are reproduced.) This explains the theoretical error estimate below.

If the data points arise from the values of a smooth (say C4) function f, i.e. fi = f(xi), then the error will behave in a predictable fashion. Let ξ be the breakpoint vector for either of the above spline interpolants. Then, the maximum absolute error satisfies

where

and ξm is the last breakpoint.

The return value for this function is a pointer of the type Imsl_f_ppoly. The calling program must receive this in a pointer Imsl_f_ppoly *ppoly. This structure contains all the information to determine the spline (stored as a piecewise polynomial) that is computed by this function. For example, the following code sequence evaluates this spline at x and returns the value in y.

The difference between the convexity-preserving spline and Akima’s spline is illustrated in the following figure. Note that the convexity-preserving interpolant exhibits linear segments where the convexity constraints are binding.

In this example, a cubic spline interpolant to a function f is computed. The values of this spline are then compared with the exact function values.

In this example, a cubic spline interpolant to a function f is computed. The values of this spline are then compared with the exact function values.