The IMSL_CSSHAPE function computes a shape-preserving cubic spline.

This routine requires an IDL Advanced Math and Stats license. For more information, contact your sales or technical support representative.

The IMSL_CSSHAPE function computes a C1 cubic spline interpolant to a set of data points (xi, fi) for the following:

i = 0, ..., (N_ELEMENTS(Xdata) – 1) = (n – 1)

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, de Boor 1978).

If the CONCAVE keyword is set, then this function computes a cubic spline interpolant to the data. For ease of explanation, xi < xi + 1 is assumed, 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 Micchelli 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 (for example, C4) function f, i.e., fi = f(xi), then the error behaves in a predictable fashion. Let ξ be the breakpoint vector for either of the above spline interpolants. Then, the maximum absolute error satisfies:


and ξm is the last breakpoint.

The returned value for this function is a structure. 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:

y = IMSL_SPVALUE(x, spline)


Example 1

In this example, a cubic spline interpolant to function values is computed. Evaluations of the computed spline are plotted along with the original data values.

x = FINDGEN(10)/9
; Define the abscissas. f = FLTARR(10)
f(0:4) = 0.25 f(5:9) = 0.75
; Define the function values. 
pp = IMSL_CSSHAPE(x, f)
; Compute the interpolant.
ppval = IMSL_SPVALUE(FINDGEN(100)/99, pp)
; Evaluate the interpolant at 100 values in [0,1].
PLOT, FINDGEN(100)/99, ppval
; Plot the results.
OPLOT, x, f, Psym = 6

Example 2

This example compares interpolants computed by IMSL_CSINTERP and IMSL_CSSHAPE with the keyword CONCAVE, as shown in the figure below.

x = [0, .1, .2, .3, .4, .5, .6, .8, 1]
y = [0, .9, .95, .9, .1, .05, .05, .2, 1]
; Define the data set and compute interpolant from IMSL_CSINTERP. pp1 = IMSL_CSINTERP(x, y)
pp2 = IMSL_CSSHAPE(x, y, /Concave)
; Compute the interpolant from IMSL_CSSHAPE with keyword Concave. 
x2 = FINDGEN(100)/99
PLOT, x2, IMSL_SPVALUE(x2, pp1), Linestyle = 2
OPLOT, x2, IMSL_SPVALUE(x2, pp2) OPLOT, x, y, Psym = 6
XYOUTS, .4, .9, 'IMSL_CSINTERP', Charsize = 1.2
OPLOT, [.73, .85], [.925, .925], Linestyle = 2
XYOUTS, .4, .8, 'IMSL_CSSHAPE !cwith CONCAVE', Charsize = 1.2
OPLOT, [.73, .85], [.8, .8]
XYOUTS, .4, .6, 'Original data', Charsize = 1.2
OPLOT, [.73], [.622], Psym = 6


Warning Errors

MATH_MAX_ITERATIONS_REACHED: Maximum number of iterations has been reached. The best approximation is returned.

Fatal Errors

MATH_DUPLICATE_XDATA_VALUES: The Xdata values must be distinct.


Result = IMSL_CSSHAPE(Xdata, Fdata [, /CONCAVE] [, /DOUBLE] [, ITMAX=value])

Return Value

A structure that represents the cubic spline interpolant.



One-dimensional array containing the abscissas of the interpolation problem.


One-dimensional array containing the ordinates for the interpolation problem.


CONCAVE (optional)

If present and nonzero, IMSL_CSSHAPE produces a cubic interpolant that preserves the concavity of the data.

DOUBLE (optional)

If present and nonzero, then double precision is used.


Allows the user to set the maximum number of iterations of Newton’s Method. To use ITMAX, the keyword CONCAVE must also be set. Default: 25.

Version History