spline2dInterp¶

Computes a two-dimensional, tensor-product spline interpolant from two-dimensional, tensor-product data.

Synopsis¶

spline2dInterp (xdata, ydata, fdata)

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 interpolated. fdata[i][j] is the value at (xdata[i], ydata[j]).

Return Value¶

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

Optional Arguments¶

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. The default knots are selected by the function splineKnots using its defaults.

Description¶

The function spline2dInterp computes a tensor-product spline interpolant. The tensor-product spline interpolant to data \(\{(x_j, y_j, f_{jj})\}\), where \(0 \leq i \leq n_x − 1\) and \(0 \leq j \leq n_y − 1\) has the form

\[\sum_{m=0}^{n_y-1} \sum_{n=0}^{n_x-1} c_{nm}B_{n,k_x,t_x}(x)B_{m,k_y,t_y}(y)\]

where \(k_x\) and \(k_y\) are the orders of the splines. These numbers are defaulted to be 4, but can be set to any positive integer using the keyword, order. Likewise, \(t_x\) and \(t_y\) are the corresponding knot sequences (xknots and yknots). These values are defaulted to the knots returned by splineKnots. The algorithm requires that

\[\begin{split}\begin{array}{l} t_x \left(k_x-1\right) \leq x_i \leq t_x\left(n_x\right)\phantom{...}0 \leq i \leq n_x-1 \\ t_y\left(k_y-1\right) \leq y_i \leq t_y\left(n_y-1\right)\phantom{...}0 \leq j \leq n_y-1 \end{array}\end{split}\]

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

\[\sum_{m=0}^{n_y-1} \sum_{n=0}^{n_x-1} c_{nm}B_{n,k_x,t_x}(x_i)B_{m,k_y,t_y}(y_j) = f_{ij}\]

Setting

\[h_{mi} = \textstyle\sum_{n=0}^{n_x-1} c_{nm}B_{n,k_x,t_x}\left(x_i\right)\]

note that for each fixed i from 1 to \(n_x − 1\), we have \(n_y − 1\) linear equations in the same number of unknowns as can be seen below:

\[\sum_{m=0}^{n_y-1} h_{mi}B_{m,k_y,t_y}\left(y_i\right) = f_{ij}\]

The same matrix appears in all of the equations above:

\[\left[B_{m,k_y,t_y}\left(y_j\right)\right] \phantom{...}1 \leq m,j \leq n_y-1\]

Thus, only factor this matrix once and then apply this factorization to the \(n_x\) right-hand sides. Once this is done and \(h_{mi}\) is computed, then solve for the coefficients \(c_{nm}\) using the relation

\[\sum_{n=0}^{n_x-1} c_{nm}B_{n,k_x,t_x}\left(x_i\right) = h_{mi}\]

for m from 0 to \(n_y − 1\), which again involves one factorization and \(n_y\) solutions to the different right-hand sides. The function spline2dInterp is based on the routine SPLI2D by de Boor (1978, p. 347).

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 (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 = spline2dValue (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.

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

# Define function


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


# Set up grid
ndata = 11
outdata = 2
xdata = empty(ndata)
ydata = empty(ndata)
fdata = empty((ndata, ndata))
for i in range(0, ndata):
    xdata[i] = float(i) / (ndata - 1)
    ydata[i] = float(i) / (ndata - 1)
for i in range(0, ndata):
    for j in range(0, ndata):
        fdata[i, j] = F(xdata[i], ydata[j])

# Compute tensor-product interpolant
sp = spline2dInterp(xdata, ydata, fdata)

# Print results
print("     x       y        F(x, y)    Interpolant    Error")
for i in range(0, outdata):
    x = float(i) / float(outdata)
    for j in range(0, outdata):
        y = float(j) / 2.0
        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)    Interpolant    Error
000   0.000       0.000        0.000       0.0000
000   0.500       0.250        0.250       0.0000
500   0.000       0.125        0.125       0.0000
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.

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

# Define function


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


# Set up grid
ndata = 11
outdata = 4
xdata = empty(ndata)
ydata = empty(ndata)
fdata = empty((ndata, ndata))
for i in range(0, ndata):
    xdata[i] = float(i) / (ndata - 1)
    ydata[i] = float(i) / (ndata - 1)
for i in range(0, ndata):
    for j in range(0, ndata):
        fdata[i, j] = F(xdata[i], ydata[j])

for order in range(3, 7, 3):
    # Compute tensor-product interpolant
    sp = spline2dInterp(xdata, ydata, fdata,
                        order={"xorder": order, "yorder": order})
    # Print results
    print("\nThe order of the spline is: ", order)
    print("     x       y        F(x, y)    Interpolant    Error")
    for i in range(0, outdata):
        x = float(i) / float(outdata)
        for j in range(0, outdata):
            y = float(j) / float(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 order of the spline is:  3
     x       y        F(x, y)    Interpolant    Error
000   0.000       0.000        0.000       0.0000
000   0.250       0.062        0.062       0.0000
000   0.500       0.250        0.250       0.0000
000   0.750       0.562        0.562       0.0000
250   0.000       0.016        0.016       0.0000
250   0.250       0.078        0.078       0.0000
250   0.500       0.266        0.266       0.0000
250   0.750       0.578        0.578       0.0000
500   0.000       0.125        0.125       0.0000
500   0.250       0.188        0.188       0.0000
500   0.500       0.375        0.375       0.0000
500   0.750       0.688        0.688       0.0000
750   0.000       0.422        0.422       0.0000
750   0.250       0.484        0.484       0.0000
750   0.500       0.672        0.672       0.0000
750   0.750       0.984        0.984       0.0000

The order of the spline is:  6
     x       y        F(x, y)    Interpolant    Error
000   0.000       0.000        0.000       0.0000
000   0.250       0.062        0.063       0.0000
000   0.500       0.250        0.250       0.0000
000   0.750       0.562        0.562       0.0000
250   0.000       0.016        0.016       0.0000
250   0.250       0.078        0.078       0.0000
250   0.500       0.266        0.266       0.0000
250   0.750       0.578        0.578       0.0000
500   0.000       0.125        0.125       0.0000
500   0.250       0.188        0.188       0.0000
500   0.500       0.375        0.375       0.0000
500   0.750       0.688        0.688       0.0000
750   0.000       0.422        0.422       0.0000
750   0.250       0.484        0.484       0.0000
750   0.500       0.672        0.672       0.0000
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 \(t_i \leq \text{data}_i < t_{i+order}\), where t is the knot sequence.
`IMSL_DATA_TOO_LARGE`	The data arrays `xdata` and `ydata` must satisfy \(\text{data}_i \leq t_{numData}\), for i = 1, …, `numData`.
`IMSL_DATA_TOO_SMALL`	The data arrays `xdata` and `ydata` must satisfy \(\text{data}_i \geq t_{order-1}\), for i = 1, …, `numData`.