spline_interp

Compute a spline interpolant.

Synopsis

#include <imsl.h>

Imsl_f_spline *imsl_f_spline_interp (int ndata, float xdata[], float fdata[], …, 0)

The type Imsl_d_spline function is imsl_d_spline_interp.

Required Arguments

int ndata (Input)
Number of data points.

float xdata[] (Input)
Array with ndata components containing the abscissas of the interpolation problem.

float fdata[] (Input)
Array with ndata components containing the ordinates of the interpolation problem.

Return Value

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

Synopsis with Optional Arguments

#include <imsl.h>

Imsl_f_spline *imsl_f_spline_interp (int ndata, float xdata[], float fdata[],

IMSL_ORDER, int order,

IMSL_KNOTS, float knots[],

Optional Arguments

IMSL_ORDER, int order (Input)
The order of the spline subspace for which the knots are desired. This option is used to communicate the order of the spline subspace.
Default: order = 4, i.e., cubic splines

IMSL_KNOTS, float knots[] (Input)
An array of size ndata + order containing the user-provided knots.
Default: knots are selected by the function imsl_f_spline_knots using its defaults.

Description

Given the data points x = xdata, f = fdata, and the number n = ndata of elements in xdata and fdata, the default action of imsl_f_spline_interp computes a cubic (k = 4) spline interpolant s to the data using the default knot sequence generated by imsl_f_spline_knots.

The optional argument IMSL_ORDER allows the user to choose the order of the spline interpolant. The optional argument IMSL_KNOTS allows user specification of knots.

The function imsl_f_spline_interp is based on the routine SPLINT by de Boor (1978, p. 204).

First, imsl_f_spline_interp sorts the xdata vector and stores the result in x. The elements of the fdata vector are permuted appropriately and stored in f, yielding the equivalent data (xi, fi) for i = 0 to n − 1.

The following preliminary checks are performed on the data. We verify that

The first test checks to see that the abscissas are distinct. The second and third inequalities verify that a valid knot sequence has been specified.

In order for the interpolation matrix to be nonsingular, we also check tk-1 ≤ xi ≤ tn for i = 0 to n − 1. This first inequality in the last check is necessary since the method used to generate the entries of the interpolation matrix requires that the k possibly nonzero B-splines at xi,

Bj-k+1, …, Bj where j satisfies tj ≤ xi < tj+1

be well-defined (that is, j − k + 1 ≥ 0).

General conditions are not known for the exact behavior of the error in spline interpolation; however, if t and x are selected properly and the data points arise from the values of a smooth (say Ck) function f, i.e. fj = f(xj), then the error will behave in a predictable fashion. The maximum absolute error satisfies

where

For more information on this problem, see de Boor (1978, Chapter 13) and his reference. This function can be used in place of the IMSL function imsl_f_cub_spline_interp.

The return value for this function is a pointer of type Imsl_f_spline. The calling program must receive this in a pointer Imsl_f_spline *sp. This structure contains all the information to determine the spline (stored as a linear combination of B-splines) 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_f_spline_value (x, sp, 0)

Three spline interpolants of order 2, 3, and 5 are plotted. These splines use the default knots.

Figure 1, Three Spline Interpolants

Examples

Example 1

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. Since the default settings are used, the interpolant is determined by the “not-a-knot” condition (see de Boor 1978).

#include <imsl.h>

#include <stdio.h>

#include <math.h>

#define NDATA 11

/* Define function */

#define F(x) (float)(sin(15.0*x))

int main()

{

int i;

float xdata[NDATA], fdata[NDATA], x, y;

Imsl_f_spline *sp;

/* Set up a grid */

for (i = 0; i < NDATA; i++) {

xdata[i] = (float)i /((float)(NDATA-1));

fdata[i] = F(xdata[i]);

}

/* Compute cubic spline interpolant */

sp = imsl_f_spline_interp (NDATA, xdata, fdata, 0);

/* Print results */

printf(" x F(x) Interpolant Error\n");

for (i = 0; i < 2*NDATA-1; i++){

x = (float) i /(float)(2*NDATA-2);

y = imsl_f_spline_value(x, sp, 0);

printf(" %6.3f %10.3f %10.3f %10.4f\n", x, F(x), y,

fabs(F(x)-y));

}

Output

x F(x) Interpolant Error

0.000 0.000 0.000 0.0000

0.050 0.682 0.809 0.1270

0.100 0.997 0.997 0.0000

0.150 0.778 0.723 0.0552

0.200 0.141 0.141 0.0000

0.250 -0.572 -0.549 0.0228

0.300 -0.978 -0.978 0.0000

0.350 -0.859 -0.843 0.0162

0.400 -0.279 -0.279 0.0000

0.450 0.450 0.441 0.0093

0.500 0.938 0.938 0.0000

0.550 0.923 0.903 0.0199

0.600 0.412 0.412 0.0000

0.650 -0.320 -0.315 0.0049

0.700 -0.880 -0.880 0.0000

0.750 -0.968 -0.938 0.0295

0.800 -0.537 -0.537 0.0000

0.850 0.183 0.148 0.0347

0.900 0.804 0.804 0.0000

0.950 0.994 1.086 0.0926

1.000 0.650 0.650 0.0000

Example 2

Recall that in the first example, a cubic spline interpolant to a function f is computed. The values of this spline are then compared with the exact function values. This example chooses to use a quadratic (k = 3) and a quintic k = 6 spline interpolant to the data instead of the default values.

#include <imsl.h>

#include <stdio.h>

#include <math.h>

#define NDATA 11

/* Define function */

#define F(x) (float)(sin(15.0*x))

int main()

{

int i, order;

float fdata[NDATA], xdata[NDATA], x, y;

Imsl_f_spline *sp;

/* Set up a grid */

for (i = 0; i < NDATA; i++) {

xdata[i] = (float)i /((float)(NDATA-1));

fdata[i] = F(xdata[i]);

}

for (order =3; order<7; order += 3) {

/* Compute cubic spline interpolant */

sp = imsl_f_spline_interp (NDATA, xdata, fdata,

IMSL_ORDER, order,

0);

/* Print results */

printf("\nThe order of the spline is %d\n", order);

printf(" x F(x) Interpolant Error\n");

for (i = NDATA/2; i < 3*NDATA/2; i++){

x = (float) i /(float)(2*NDATA-2);

y = imsl_f_spline_value(x,sp,0);

printf(" %6.3f %10.3f %10.3f %10.4f\n", x, F(x), y,

fabs(F(x)-y));

}

Output

The order of the spline is 3

x F(x) Interpolant Error

0.250 -0.572 -0.542 0.0299

0.300 -0.978 -0.978 0.0000

0.350 -0.859 -0.819 0.0397

0.400 -0.279 -0.279 0.0000

0.450 0.450 0.429 0.0210

0.500 0.938 0.938 0.0000

0.550 0.923 0.879 0.0433

0.600 0.412 0.412 0.0000

0.650 -0.320 -0.305 0.0149

0.700 -0.880 -0.880 0.0000

0.750 -0.968 -0.922 0.0459

The order of the spline is 6

x F(x) Interpolant Error

0.250 -0.572 -0.573 0.0016

0.300 -0.978 -0.978 0.0000

0.350 -0.859 -0.856 0.0031

0.400 -0.279 -0.279 0.0000

0.450 0.450 0.448 0.0020

0.500 0.938 0.938 0.0000

0.550 0.923 0.922 0.0003

0.600 0.412 0.412 0.0000

0.650 -0.320 -0.322 0.0025

0.700 -0.880 -0.880 0.0000

0.750 -0.968 -0.959 0.0090

Warning Errors

IMSL_ILL_COND_INTERP_PROB

The interpolation matrix is ill-conditioned. The solution might not be accurate.

Fatal Errors

IMSL_DUPLICATE_XDATA_VALUES	The xdata values must be distinct.
IMSL_KNOT_MULTIPLICITY	Multiplicity of the knots cannot exceed the order of the spline.
IMSL_KNOT_NOT_INCREASING	The knots must be nondecreasing.
IMSL_KNOT_DATA_INTERLACING	The i-th smallest element of xdata (xi) must satisfy ti ≤ xi < ti+order where t is the knot sequence.
IMSL_XDATA_TOO_LARGE	The array xdata must satisfy xdatai ≤ tndata, for i = 0, …, ndata-1.
IMSL_XDATA_TOO_SMALL	The array xdata must satisfy xdatai ≥ torder-1, for i = 0, …, ndata-1.