spline_2d_least_squares
Computes a two-dimensional, tensor-product spline approximant using least-squares.
Synopsis
#include <imsl.h>
Imsl_f_spline *imsl_f_spline_2d_least_squares (int num_xdata, float xdata[], int num_ydata, float ydata[], float fdata[], int x_spline_space_dim, int y_spline_space_dim, …, 0)
The type Imsl_d_spline function is imsl_d_spline_2d_least_squares.
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 approximated. fdata[i][j] is the value at (xdata[i], ydata[j]).
int x_spline_space_dim (Input)
The linear dimension of the spline subspace for the x variable. It should be smaller than num_xdata and greater than or equal to xorder whose default value is 4.
int y_spline_space_dim (Input)
The linear dimension of the spline subspace for the y variable. It should be smaller than num_ydata and greater than or equal to yorder whose default value is 4.
Return Value
This is 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_least_squares (int num_xdata, float xdata[], int num_ydata, float ydata[], float fdata[], int x_spline_space_dim, int y_spline_space_dim,
IMSL_SSE, float *sse,
IMSL_ORDER, int xorder, int yorder,
IMSL_KNOTS, float xknots[], float yknots[],
IMSL_FDATA_COL_DIM, int fdata_col_dim,
IMSL_WEIGHTS, float xweights[], float yweights[],
0)
Optional Arguments
IMSL_SSE, float *sse (Output)
This option places the weighted error sum of squares in the place pointed to by sse.
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)
This option requires the user to provide the knots.
Default: The default knots are equally spaced in the x and y dimensions.
IMSL_FDATA_COL_DIM, int fdata_col_dim (Input)
The column dimension of fdata.
Default: fdata_col_dim = num_ydata
IMSL_WEIGHTS, float xweights[], float yweights[] (Input)
This option requires the user to provide the weights for the least-squares fit.
Default: all weights are equal to 1.
Description
The imsl_f_spline_2d_least_squares procedure computes a tensor-product spline least-squares approximation to weighted tensor-product data. The input for this function consists of data vectors to specify the tensor-product grid for the data, two vectors with the weights (optional, the default is 1), the values of the surface on the grid, and the specification for the tensor-product spline (optional, a default is chosen). The grid is specified by the two arrays x = xdata and y = ydata of length n = num_xdata and m = num_ydata, respectively. A two-dimensional array f = fdata contains the data values which are to be fit. The two vectors wx = xweights and wy = yweights contain the weights for the weighted least-squares problem. The information for the approximating tensor-product spline can be provided using the keywords IMSL_ORDER and IMSL_KNOTS. This information is contained in kx = xorder, tx = xknots, and N = xspline_space_dim for the spline in the first variable, and in ky = yorder, ty = yknots and M = y_spline_space_dim for the spline in the second variable.
This function computes coefficients for the tensor-product spline by solving the normal equations in tensor-product form as discussed in de Boor (1978, Chapter 17). Also see the paper by Grosse (1980).
As the computation proceeds, we obtain coefficients c minimizing
where the function Bkl is the tensor-product of two B-splines of order kx and ky. Specifically, we have
The spline
and its partial derivatives can be evaluated at any point (
x, y) using
imsl_f_spline_2d_value.
The return value for this function is a pointer to the structure Imsl_f_spline. The calling program must receive this in a pointer of type Imsl_f_spline. This structure contains all the information to determine the spline that is computed by this procedure. For example, the following code sequence evaluates this spline (stored in the structure sp at (x, y) and returns the value in v.
v = imsl_f_spline_2d_value (x, y, sp, 0)
Examples
Example 1
The data for this example comes from the function ex sin (x + y) on the rectangle [0, 3] × [0, 5]. This function is sampled on a 50 × 25 grid. Then recover or smooth the data by using tensor-product cubic splines and least-squares fitting. The values of the function ex sin (x + y) are printed on a 2 × 2 grid and compared with the values of the tensor-product spline least-squares fit.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define NXDATA 50
#define NYDATA 25
#define OUTDATA 2
/* Define function */
#define F(x,y) (float)(exp(x)*sin(x+y))
int main()
{
int i, j, num_xdata, num_ydata;
float fdata[NXDATA][NYDATA];
float xdata[NXDATA], ydata[NYDATA];
Imsl_f_spline *sp;
/* Set up grid */
for (i = 0; i < NXDATA; i++) {
xdata[i] = 3.*(float) i / ((float)(NXDATA-1));
}
for (i = 0; i < NYDATA; i++) {
ydata[i] = 5.*(float) i / ((float)(NYDATA-1));
}
/* Compute function values on grid */
for (i = 0; i < NXDATA; i++) {
for (j = 0; j < NYDATA; j++) {
fdata[i][j] = F(xdata[i], ydata[j]);
}
}
num_xdata = NXDATA;
num_ydata = NYDATA;
/* Compute tensor-product interpolant */
sp = imsl_f_spline_2d_least_squares(num_xdata, &xdata[0], num_ydata,
&ydata[0], &fdata[0][0], 5, 7, 0);
/* Print results */
printf(" x y F(x, y) Spline Fit Abs. 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) Spline Fit Abs. Error
0.000 0.000 0.000 -0.020 0.0204
0.000 0.500 0.479 0.500 0.0208
0.500 0.000 0.790 0.816 0.0253
0.500 0.500 1.387 1.384 0.0031
Example 2
The same data is used as in the previous example. Optional argument IMSL_SSE is used to return the error sum of squares.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define NXDATA 50
#define NYDATA 25
#define OUTDATA 2
/* Define function */
#define F(x,y) (float)(exp(x)*sin(x+y))
int main()
{
int i, j, num_xdata, num_ydata;
float fdata[NXDATA][NYDATA];
float xdata[NXDATA], ydata[NYDATA], sse, x, y, z;
Imsl_f_spline *sp;
/* Set up grid */
for (i = 0; i < NXDATA; i++) {
xdata[i] = 3.*(float) i / ((float) (NXDATA - 1));
}
for (i = 0; i < NYDATA; i++) {
ydata[i] = 5.*(float) i / ((float) (NYDATA - 1));
}
/* Compute function values on grid */
for (i = 0; i < NXDATA; i++) {
for (j = 0; j < NYDATA; j++) {
fdata[i][j] = F(xdata[i], ydata[j]);
}
}
num_xdata = NXDATA;
num_ydata = NYDATA;
/* Compute tensor-product interpolant */
sp = imsl_f_spline_2d_least_squares(num_xdata, &xdata[0], num_ydata,
&ydata[0], &fdata[0][0], 5, 7,
IMSL_SSE, &sse, 0);
/* Print results */
printf("The error sum of squares is %10.3f\n\n", sse);
printf(" x y F(x, y) Spline Fit Abs. 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 error sum of squares is 3.753
x y F(x, y) Spline Fit Abs. Error
0.000 0.000 0.000 -0.020 0.0204
0.000 0.500 0.479 0.500 0.0208
0.500 0.000 0.790 0.816 0.0253
0.500 0.500 1.387 1.384 0.0031
Warning Errors
IMSL_ILL_COND_LSQ_PROB | The least-squares matrix is ill-conditioned. The solution might not be accurate. |
IMSL_SPLINE_LOW_ACCURACY | There may be less than one digit of accuracy in the least-squares fit. Try using a higher precision if possible. |
Fatal Errors
IMSL_KNOT_MULTIPLICITY | Multiplicity of the knots cannot exceed the order of the spline. |
IMSL_KNOT_NOT_INCREASING | The knots must be nondecreasing. |
IMSL_SPLINE_LRGST_ELEMNT | The data arrays xdata and ydata must satisfy datai ≤ tspline_space_dim, for i = 1, …, num_data. |
IMSL_SPLINE_SMLST_ELEMNT | The data arrays xdata and ydata must satisfy datai ≥ torder-1, for i = 1, …, num_data. |
IMSL_NEGATIVE_WEIGHTS | All weights must be greater than or equal to zero. |
IMSL_DATA_DECREASING | The xdata values must be nondecreasing. |