spline_2d_interp
Computes a two-dimensional, tensor-product spline interpolant from two-dimensional, tensor-product data.
Synopsis
#include <imsl.h>
Imsl_f_spline *imsl_f_spline_2d_interp (int num_xdata, float xdata[], int num_ydata, float ydata[], float fdata[], …, 0)
The type Imsl_d_spline function is imsl_d_spline_2d_interp.
Required Arguments
int num_xdata (Input)
Number of data points in the X direction.
float xdata[] (Input)
Array with num_xdata components containing the data points in the X direction.
int num_ydata (Input)
Number of data points in the Y direction.
float ydata[] (Input)
Array with num_ydata components containing the data points in the Y direction.
float fdata[] (Input)
Array of size num_xdata × num_ydata containing the values to be interpolated. fdata[i][j] is the value at (xdata[i], ydata[j]).
Return Value
A pointer to the structure that represents the tensor-product 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_2d_interp (int num_xdata, float xdata[], int num_ydata, float ydata[], float fdata[],
IMSL_ORDER, int xorder, int yorder,
IMSL_KNOTS, float xknots[], float yknots[],
IMSL_FDATA_COL_DIM, int fdata_col_dim,
0)
Optional Arguments
IMSL_ORDER, int xorder, int yorder (Input)
This option is used to communicate the order of the spline subspace.
Default: xorder, yorder = 4, (i.e., tensor-product cubic splines)
IMSL_KNOTS, float xknots[], float yknots[] (Input)
User-provided arrays containing the knots. Array xknots must be of size num_xdata + xorder and array yknots of size num_ydata + yorder. The default knots are selected by the function imsl_f_spline_knots using its defaults.
IMSL_FDATA_COL_DIM, int fdata_col_dim (Input)
The column dimension of the matrix fdata.
Default: fdata_col_dim = num_ydata
Description
The function imsl_f_spline_2d_interp computes a tensor-product spline interpolant. The tensor-product spline interpolant to data {(xj, yj, fjj)}, where 0 ≤ i ≤ nx − 1 and 0 ≤ j ≤ ny − 1 has the form
where kx and ky are the orders of the splines. These numbers are defaulted to be 4, but can be set to any positive integer using the keyword, IMSL_ORDER. Likewise, tx and ty are the corresponding knot sequences (xknots and yknots). These values are defaulted to the knots returned by imsl_f_spline_knots. The algorithm requires that
Tensor-product spline interpolants in two dimensions can be computed quite efficiently by solving (repeatedly) two univariate interpolation problems.
The computation is motivated by the following observations. It is necessary to solve the system of equations
Setting
note that for each fixed i from 1 to nx − 1, we have ny − 1 linear equations in the same number of unknowns as can be seen below:
The same matrix appears in all of the equations above:
Thus, only factor this matrix once and then apply this factorization to the nx right-hand sides. Once this is done and hmi is computed, then solve for the coefficients cnm using the relation
for m from 0 to ny − 1, which again involves one factorization and ny solutions to the different right-hand sides. The function imsl_f_spline_2d_interp is based on the routine SPLI2D by de Boor (1978, p. 347).
The return value for this function is a pointer to the structure 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 in B-spline format) that is computed by this procedure. For example, the following code sequence evaluates this spline at (x,y) and returns the value in z.
z = imsl_f_spline_2d_value (x, y, sp, 0);
Examples
Example 1
In this example, a tensor-product spline interpolant to a function f is computed. The values of the interpolant and the error on a 4 × 4 grid are displayed.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define NDATA 11
#define OUTDATA 2
/* Define function */
#define F(x, y) (float)(x*x*x+y*y)
int main()
{
int i, j, num_xdata, num_ydata;
float fdata[NDATA][NDATA], xdata[NDATA], ydata[NDATA];
float x, y, z;
Imsl_f_spline *sp;
/* Set up grid */
for (i = 0; i < NDATA; i++) {
xdata[i] = ydata[i] = (float)i / ((float)(NDATA-1));
}
for (i = 0; i < NDATA; i++) {
for (j = 0; j < NDATA; j++) {
fdata[i][j] = F(xdata[i], ydata[j]);
}
}
num_xdata = num_ydata = NDATA;
/* Compute tensor-product interpolant */
sp = imsl_f_spline_2d_interp(num_xdata, xdata, num_ydata,
ydata, (float*) fdata, 0);
/* Print results */
printf(" x y F(x, y) Interpolant Error \n");
for (i = 0; i < OUTDATA; i++) {
x = (float) i / (float) (OUTDATA);
for (j = 0; j < OUTDATA; j++) {
y = (float) j / (float) (OUTDATA);
z = imsl_f_spline_2d_value(x, y, sp, 0);
printf(" %6.3f %6.3f %10.3f %10.3f %10.4f\n",
x, y, F(x,y), z, fabs(F(x,y)-z));
}
}
}
Output
x y F(x, y) Interpolant Error
0.000 0.000 0.000 0.000 0.0000
0.000 0.500 0.250 0.250 0.0000
0.500 0.000 0.125 0.125 0.0000
0.500 0.500 0.375 0.375 0.0000
Example 2
Recall that in the first example, a tensor-product spline interpolant to a function f is computed. The values of the interpolant and the error on a 4 × 4 grid are displayed. Notice that the first interpolant with order = 3 does not reproduce the cubic data, while the second interpolant with order = 6 does reproduce the data.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define NDATA 7
#define OUTDATA 4
/* Define function */
#define F(x,y) (float)(x*x*x+y*y)
int main()
{
int i, j, num_xdata, num_ydata, order;
float fdata[NDATA][NDATA], xdata[NDATA], ydata[NDATA];
float x, y, z;
Imsl_f_spline *sp;
/* Set up grid */
for (i = 0; i < NDATA; i++) {
xdata[i] = ydata[i] = (float) i / ((float) (NDATA - 1));
}
for (i = 0; i < NDATA; i++) {
for (j = 0; j < NDATA; j++) {
fdata[i][j] = F(xdata[i], ydata[j]);
}
}
num_xdata = num_ydata = NDATA;
for(order = 3; order < 7; order += 3) {
/* Compute tensor-product interpolant */
sp = imsl_f_spline_2d_interp(num_xdata, xdata, num_ydata,
ydata, (float *)fdata,
IMSL_ORDER, order, order,
0);
/* Print results */
printf("\nThe order of the spline is %d \n", order);
printf(" x y F(x, y) Interpolant Error\n");
for (i = 0; i < OUTDATA; i++) {
x = (float) i / (float) (OUTDATA);
for (j = 0; j < OUTDATA; j++) {
y = (float) j / (float) (OUTDATA);
z = imsl_f_spline_2d_value(x, y, sp, 0);
printf(" %6.3f %6.3f %10.3f %10.3f %10.4f \n",
x, y, F(x,y), z, fabs(F(x,y)-z));
}
}
}
}
Output
The order of the spline is 3
x y F(x, y) Interpolant Error
0.000 0.000 0.000 0.000 0.0000
0.000 0.250 0.062 0.063 0.0000
0.000 0.500 0.250 0.250 0.0000
0.000 0.750 0.562 0.562 0.0000
0.250 0.000 0.016 0.016 0.0002
0.250 0.250 0.078 0.078 0.0002
0.250 0.500 0.266 0.266 0.0002
0.250 0.750 0.578 0.578 0.0002
0.500 0.000 0.125 0.125 0.0000
0.500 0.250 0.188 0.188 0.0000
0.500 0.500 0.375 0.375 0.0000
0.500 0.750 0.688 0.687 0.0000
0.750 0.000 0.422 0.422 0.0002
0.750 0.250 0.484 0.484 0.0002
0.750 0.500 0.672 0.672 0.0002
0.750 0.750 0.984 0.984 0.0002
The order of the spline is 6
x y F(x, y) Interpolant Error
0.000 0.000 0.000 0.000 0.0000
0.000 0.250 0.062 0.063 0.0000
0.000 0.500 0.250 0.250 0.0000
0.000 0.750 0.562 0.562 0.0000
0.250 0.000 0.016 0.016 0.0000
0.250 0.250 0.078 0.078 0.0000
0.250 0.500 0.266 0.266 0.0000
0.250 0.750 0.578 0.578 0.0000
0.500 0.000 0.125 0.125 0.0000
0.500 0.250 0.188 0.188 0.0000
0.500 0.500 0.375 0.375 0.0000
0.500 0.750 0.688 0.688 0.0000
0.750 0.000 0.422 0.422 0.0000
0.750 0.250 0.484 0.484 0.0000
0.750 0.500 0.672 0.672 0.0000
0.750 0.750 0.984 0.984 0.0000
Warning Errors
IMSL_ILL_COND_INTERP_PROB | The interpolation matrix is ill-conditioned. The solution might not be accurate. |
Fatal Errors
IMSL_XDATA_NOT_INCREASING | The xdata values must be strictly increasing. |
IMSL_YDATA_NOT_INCREASING | The ydata values must be strictly increasing. |
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 the data arrays xdata and ydata must satisfy ti ≤ datai < ti+order, where t is the knot sequence. |
IMSL_DATA_TOO_LARGE | The data arrays xdata and ydata must satisfy datai ≤ tnum_data, for i = 1, …, num_data. |
IMSL_DATA_TOO_SMALL | The data arrays xdata and ydata must satisfy datai ≥ torder-1, for i = 1, …, num_data. |
