cub_spline_interp_shape
Computes a shape-preserving cubic spline.
Synopsis
#include <imsl.h>
Imsl_f_ppoly *imsl_f_cub_spline_interp_shape (int ndata, float xdata[], float fdata[], …, 0)
The type Imsl_d_ppoly function is imsl_d_cub_spline_interp_shape.
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 for the interpolation problem.
Return Value
A pointer to the structure that represents the cubic 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_ppoly *imsl_f_cub_spline_interp_shape (int ndata, float xdata[], float fdata[],
IMSL_CONCAVE,
IMSL_CONCAVE_ITMAX, int itmax,
0)
Optional Arguments
IMSL_CONCAVE
This option produces a cubic interpolant that will preserve the concavity of the data.
IMSL_CONCAVE_ITMAX, int itmax (Input)
This option allows the user to set the maximum number of iterations of Newton’s Method.
Default: itmax = 25.
Description
The function imsl_f_cub_spline_interp_shape computes a C1 cubic spline interpolant to a set of data points(xi, fi) for i = 0, …, ndata − 1 = n. 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) or de Boor (1978).
If the optional argument IMSL_CONCAVE is chosen, then this function computes a cubic spline interpolant to the data. For ease of explanation, we will assume that xi < xi+1, 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 Michelli 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 (say C4) function f, i.e. fi = f(xi), then the error will behave in a predictable fashion. Let ξ be the breakpoint vector for either of the above spline interpolants. Then, the maximum absolute error satisfies
where
and ξm is the last breakpoint.
The return value for this function is a pointer of the type Imsl_f_ppoly. The calling program must receive this in a pointer Imsl_f_ppoly *ppoly. 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_f_cub_spline_value (x, ppoly, 0)
The difference between the convexity-preserving spline and Akima’s spline is illustrated in the following figure. Note that the convexity-preserving interpolant exhibits linear segments where the convexity constraints are binding.
Figure 3.2 — Two Shape-Preserving Splines
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.
#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 fdata[NDATA], xdata[NDATA], x, y;
Imsl_f_ppoly *pp;
/* Compute xdata and fdata */
for (i = 0; i < NDATA; i++) {
xdata[i] = (float)(i)/(NDATA-1);
fdata[i] = F(xdata[i]);
}
/* Compute cubic spline interpolant */
pp = imsl_f_cub_spline_interp_shape(NDATA, xdata, fdata, 0);
/* Print results */
printf(" x F(x) Interpolant Error\n\n");
for (i = 0; i < 2*NDATA-1; i++) {
x = (float) i /(float)(2*NDATA-2);
y = imsl_f_cub_spline_value(x, pp, 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.818 0.1360
0.100 0.997 0.997 0.0000
0.150 0.778 0.615 0.1635
0.200 0.141 0.141 0.0000
0.250 -0.572 -0.478 0.0934
0.300 -0.978 -0.978 0.0000
0.350 -0.859 -0.812 0.0464
0.400 -0.279 -0.279 0.0000
0.450 0.450 0.386 0.0645
0.500 0.938 0.938 0.0000
0.550 0.923 0.854 0.0683
0.600 0.412 0.412 0.0000
0.650 -0.320 -0.276 0.0433
0.700 -0.880 -0.880 0.0000
0.750 -0.968 -0.889 0.0789
0.800 -0.537 -0.537 0.0000
0.850 0.183 0.149 0.0338
0.900 0.804 0.804 0.0000
0.950 0.994 0.932 0.0613
1.000 0.650 0.650 0.0000
Example 2
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.
#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 fdata[NDATA], xdata[NDATA], x, y;
Imsl_f_ppoly *pp;
/* Compute xdata and fdata */
for (i = 0; i < NDATA; i++) {
xdata[i] = (float)(i)/(NDATA-1);
fdata[i] = F(xdata[i]);
}
/* Compute cubic spline interpolant */
pp = imsl_f_cub_spline_interp_shape(NDATA, xdata, fdata,
IMSL_CONCAVE,
0);
/* Print results */
printf(" x F(x) Interpolant Error\n\n");
for (i = 0; i < 2*NDATA-1; i++){
x = (float) i /(float)(2*NDATA-2);
y = imsl_f_cub_spline_value(x, pp, 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.667 0.0150
0.100 0.997 0.997 0.0000
0.150 0.778 0.761 0.0172
0.200 0.141 0.141 0.0000
0.250 -0.572 -0.559 0.0126
0.300 -0.978 -0.978 0.0000
0.350 -0.859 -0.840 0.0189
0.400 -0.279 -0.279 0.0000
0.450 0.450 0.440 0.0098
0.500 0.938 0.938 0.0000
0.550 0.923 0.902 0.0208
0.600 0.412 0.412 0.0000
0.650 -0.320 -0.311 0.0086
0.700 -0.880 -0.880 0.0000
0.750 -0.968 -0.952 0.0156
0.800 -0.537 -0.537 0.0000
0.850 0.183 0.200 0.0174
0.900 0.804 0.804 0.0000
0.950 0.994 0.892 0.1020
1.000 0.650 0.650 0.0000
Warning Errors
IMSL_MAX_ITERATIONS_REACHED | The maximum number of iterations has been reached. The best approximation is returned. |
Fatal Errors
IMSL_DUPLICATE_XDATA_VALUES | The xdata values must be distinct. |
