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[],
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_KNOTS, float knots[] (Input)
This option requires the user to provide the 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 3.4 — 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_XDATA_INTERLACING | The i-th smallest element of xdata (xj) must satisfy tj ≤ xj < tj+order where t is the knot sequence. |
IMSL_XDATA_TOO_LARGE | The array xdata must satisfy xdataj ≤ tndata, for i = 1, …, ndata. |
IMSL_XDATA_TOO_SMALL | The array xdata must satisfy xdataj ≥ torder-1, for i = 1, …, ndata. |
