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∥≤c∥f (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. |
