spline_knots

Computes the knots for a spline interpolant

Synopsis

#include <imsl.h>

float *imsl_f_spline_knots (int ndata, float xdata[], 0)

The type double function is imsl_d_spline_knots.

Required Arguments

int ndata (Input)
Number of data points.

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

Return Value

A pointer to the knots. If the knots cannot be computed, then NULL is returned. To release this space, use imsl_free.

Synopsis with Optional Arguments

#include <imsl.h>

float *imsl_f_spline_knots (int ndata, float xdata[],

IMSL_ORDER, int order,

IMSL_OPT,

IMSL_OPT_ITMAX, int itmax,

IMSL_RETURN_USER, float knots[],

0)

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_OPT
This option produces knots that satisfy an optimality criterion.

IMSL_OPT_ITMAX, int itmax (Input)
This option allows the user to set the maximum number of iterations of Newton’s method.
Default: itmax = 10

IMSL_RETURN_USER, float knots[] (Output)
A user-provided array of size ndata + order containing the knots. For example, the user could declare float knots[100]; and pass in knots. The return value is then also set to knots.

Description

Given the data points x = xdata, the order of the spline k = order, and the number n = ndata of elements in xdata, the default action of imsl_f_spline_knots returns a pointer to a knot sequence that is appropriate for interpolation of data on x by splines of order k (the default order is k = 4). The knot sequence is contained in its first n + k positions. If k is even, and we assume that the entries in the input vector x are increasing, then the resulting knot sequence t is returned as

 

There is some discussion concerning this selection of knots in de Boor (1978, p. 211). If k is odd, then t is returned as

 

It is not necessary to sort the values in xdata.

If the option IMSL_OPT is selected, then the knot sequence returned minimizes the constant c in the error estimate

f - s∥≤cf (k)

In the above formula, f is any function in Ck, and s is the spline interpolant to f at the abscissas x with knot sequence t.

The algorithm is based on a routine described in de Boor (1978, p. 204), which in turn is based on a theorem of Micchelli et al. (1976).

Examples

 

Example 1

In this example, knots for a cubic spline are generated and printed. Notice that the knots are stacked at the endpoints and that the second and next to last data points are not knots.

 

#include <imsl.h>

#include <stdio.h>

#include <math.h>

 

#define NDATA 6

 

int main()

{

int i;

float *knots, xdata[NDATA];

 

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

xdata[i] = i;

knots = imsl_f_spline_knots(NDATA, xdata, 0);

imsl_f_write_matrix("The knots for the cubic spline are:\n",

1, NDATA+4, knots,

IMSL_COL_NUMBER_ZERO,

0);

}

Output

 

The knots for the cubic spline are:

 

0 1 2 3 4 5

0 0 0 0 2 3

 

6 7 8 9

5 5 5 5

Example 2

This is a continuation of the examples for imsl_f_spline_interp. Recall that in these examples, a cubic spline interpolant to a function f is computed first. The values of this spline are then compared with the exact function values. The second example uses a quadratic (k = 3) and a quintic (k = 6) spline interpolant to the data. Now, instead of using the default knots, select the “optimal” knots as described above. Notice that the error is actually worse in this case.

 

#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], *knots, 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) {

knots = imsl_f_spline_knots(NDATA, xdata, IMSL_ORDER, order,

IMSL_OPT,

0);

/* Compute spline interpolant */

sp = imsl_f_spline_interp (NDATA, xdata,fdata,

IMSL_ORDER, order,

IMSL_KNOTS, knots,

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.543 0.0290

0.300 -0.978 -0.978 0.0000

0.350 -0.859 -0.819 0.0401

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.0150

0.700 -0.880 -0.880 0.0000

0.750 -0.968 -0.920 0.0478

 

The order of the spline is 6

x F(x) Interpolant Error

0.250 -0.572 -0.578 0.0061

0.300 -0.978 -0.978 0.0000

0.350 -0.859 -0.854 0.0054

0.400 -0.279 -0.279 0.0000

0.450 0.450 0.448 0.0019

0.500 0.938 0.938 0.0000

0.550 0.923 0.920 0.0022

0.600 0.412 0.412 0.0000

0.650 -0.320 -0.317 0.0020

0.700 -0.880 -0.880 0.0000

0.750 -0.968 -0.966 0.0023

Warning Errors

IMSL_NO_CONV_NEWTON

Newton’s method iteration did not converge.

Fatal Errors

IMSL_DUPLICATE_XDATA_VALUES

The xdata values must be distinct.

IMSL_ILL_COND_LIN_SYS

Interpolation matrix is singular. The xdata values may be too close together.