spline2dLeastSquares¶

Computes a two-dimensional, tensor-product spline approximate using least-squares.

Synopsis¶

spline2dLeastSquares (xdata, ydata, fdata, xSplineSpaceDim, ySplineSpaceDim)

Required Arguments¶

float xdata[] (Input): Array with numXdata components containing the data points in the x direction.
float ydata[] (Input): Array with numYdata components containing the data points in the y direction.
float fdata[[]] (Input): Array of size numXdata × numYdata containing the values to be approximated. fdata[i][j] is the value at (xdata[i], ydata[j]).
int xSplineSpaceDim (Input): The linear dimension of the spline subspace for the x variable. It should be smaller than numXdata and greater than or equal to xorder whose default value is 4.
int ySplineSpaceDim (Input): The linear dimension of the spline subspace for the y variable. It should be smaller than numYdata and greater than or equal to yorder whose default value is 4.

Return Value¶

This is the structure that represents the tensor-product spline interpolant. If an interpolant cannot be computed, then None is returned.

Optional Arguments¶

sse (Output)

This option places the weighted error sum of squares in the place pointed to by sse.

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)

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.

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 spline2dLeastSquares 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 = numXdata and m = numYdata, respectively. A two-dimensional array f = fdata contains the data values which are to be fit. The two vectors \(w_x\) = xweights and \(w_y\) = yweights contain the weights for the weighted least-squares problem. The information for the approximating tensor-product spline can be provided using the keywords order and knots. This information is contained in \(k_x\) = xorder, \(t_x\) = xknots, and N = xspline_space_dim for the spline in the first variable, and in \(k_y\) = yorder, \(t_y\) = yknots and M = ySplineSpaceDim 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

\[\sum_{i=0}^{n=1} \sum_{i=0}^{m-1} w_x(i)w_y(j) \left[ \sum_{k=0}^{N-1} \sum_{l=0}^{M-1} c_{kl} B_{kl} \left( x_i,y_i \right) - f_{ij} \right]^2\]

where the function \(B_{kl}\) is the tensor-product of two B-splines of order \(k_x\) and \(k_y\). Specifically, we have

\[B_{kl}(x,y) = B_{k,k_x,t_x}(x) B_{l,k_y,t_y}(y)\]

The spline

\[\sum_{k=0}^{N-1} \sum_{l=0}^{M-1} c_{kl}B_{kl}(x,y)\]

and its partial derivatives can be evaluated at any point (x, y) using spline2dValue.

The return value for this function is the structure Imsl_d_spline. The calling program must receive this structure. 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 ppoly at (x, y) and returns the value in v.

v = spline2dValue (x, y, ppoly)

Examples¶

Example 1¶

The data for this example comes from the function \(e^x \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 \(e^x \sin(x + y)\) are printed on a 2 × 2 grid and compared with the values of the tensor-product spline least-squares fit.

from __future__ import print_function
from numpy import *
from pyimsl.math.spline2dLeastSquares import spline2dLeastSquares
from pyimsl.math.spline2dValue import spline2dValue

# Define function


def F(x, y):
    return exp(x) * sin(x + y)


# Set up grid
nxdata = 50
nydata = 25
outdata = 2
xdata = empty(nxdata)
ydata = empty(nydata)
fdata = empty((nxdata, nydata))
for i in range(0, nxdata):
    xdata[i] = 3. * i / (nxdata - 1)
for i in range(0, nydata):
    ydata[i] = 5. * i / (nydata - 1)

#  Compute function values on grid
for i in range(0, nxdata):
    for j in range(0, nydata):
        fdata[i, j] = F(xdata[i], ydata[j])

# Compute tensor-product interpolant
sp = spline2dLeastSquares(xdata, ydata, fdata, 5, 7)

# Print results
print("     x       y        F(x, y)   Fitted Values    Error")
for i in range(0, outdata):
    x = float(i) / outdata
    for j in range(0, outdata):
        y = float(j) / outdata
        z = spline2dValue(x, y, sp)
        print("  %6.3f  %6.3f  %10.3f   %10.3f   %10.4f"
              % (x, y, F(x, y), z, abs(F(x, y) - z)))

Output¶

     x       y        F(x, y)   Fitted Values    Error
000   0.000       0.000       -0.020       0.0204
000   0.500       0.479        0.500       0.0208
500   0.000       0.790        0.816       0.0253
500   0.500       1.387        1.384       0.0031

Example 2¶

The same data is used as in the previous example. Optional argument sse is used to return the error sum of squares.

from __future__ import print_function
from numpy import *
from pyimsl.math.spline2dLeastSquares import spline2dLeastSquares
from pyimsl.math.spline2dValue import spline2dValue

# Define function


def F(x, y):
    return exp(x) * sin(x + y)


# Set up grid
nxdata = 50
nydata = 25
outdata = 2
xdata = empty(nxdata)
ydata = empty(nydata)
fdata = empty((nxdata, nydata))
for i in range(0, nxdata):
    xdata[i] = 3. * i / (nxdata - 1)
for i in range(0, nydata):
    ydata[i] = 5. * i / (nydata - 1)

#  Compute function values on grid
for i in range(0, nxdata):
    for j in range(0, nydata):
        fdata[i, j] = F(xdata[i], ydata[j])

# Compute tensor-product interpolant
x = []
sp = spline2dLeastSquares(xdata, ydata, fdata, 5, 7, sse=x)

# Print results
print("The error sum of squares is %10.3f\n" % (x[0]))
print("       x       y      F(x, y)   Fitted Values    Error")
for i in range(0, outdata):
    x = float(i) / outdata
    for j in range(0, outdata):
        y = float(j) / outdata
        z = spline2dValue(x, y, sp)
        print("  %6.3f  %6.3f  %10.3f   %10.3f   %10.4f"
              % (x, y, F(x, y), z, abs(F(x, y) - z)))

Output¶

The error sum of squares is      3.753

       x       y      F(x, y)   Fitted Values    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 \(\text{data}_i \leq t_{splineSpaceDim}\), for i = 1, …, `numData`.
`IMSL_SPLINE_SMLST_ELEMNT`	The data arrays `xdata` and `ydata` must satisfy \(\text{data}_i \geq t_{order-1}\), for i = 1, …, `numData`.
`IMSL_NEGATIVE_WEIGHTS`	All weights must be greater than or equal to zero.
`IMSL_DATA_DECREASING`	The `xdata` values must be nondecreasing.