splineNdInterp¶

Performs multidimensional interpolation and differentiation for up to 7 dimensions.

Synopsis¶

float splineNdInterp (d, x, xdata, fdata)

Required Arguments¶

int d[] (Input)

Array of length n. d[i] contains the number of gridpoints in the i-th direction.

float x[] (Input)

Array of length n containing the point at which interpolation is to be done. An interpolant is to be calculated at the point:

\[\left(X_1,X_2, \ldots, X_n\right)\]

where

\[X_1 = x[0], X_2 = x[1], \ldots\]

float xdata[[]] (Input)

Array of size n × max(d[0], …, d[n-1]) containing the gridpoint values for the grid.

float fdata[] (Input)

Array of size d[0]× d[1]× …× d[n-1] containing the values of the function to be interpolated at the gridpoints.

fdata(i, j, k, …) is the value of the function at

\[\left(Z_{1,i}, Z_{2,j}, Z_{3,k}, \ldots\right)\]

where

\[Z_{1,i} = \mathit{xdata}[0][i-1]\]

\[Z_{2,j} = \mathit{xdata}[1][j-1]\]

\[Z_{3,k} = \mathit{xdata}[2][k-1]\]

\(\text{for } i = 1, \ldots, d[0], j = 1, \ldots, d[1], k = 1, \ldots, d[2], \ldots\)

Return Value¶

Interpolated value of the function. If no value can be computed, NaN is returned.

Optional Arguments¶

nDegree, int (Input)

Array of length n containing the degree of polynomial interpolation to be used in each dimension. nDegree[i] must be less than or equal to 15.

Default: nDegree[i] = 5, for i = 0, ..., n‑1.

order, int (Input)

Maximum order of derivatives to be computed with respect to each variable. order cannot be larger than max (7‑ n, 2). All partial derivatives up to and including order order are returned in each of the n dimensions. See deriv for more details.

Default: order = 0.

deriv, float (Output)

An n dimensional array, dimensioned (order +1) × (order +1) × …, containing derivative estimates at the interpolation point. deriv [i] [j] … will hold an estimate of the derivative with respect to \(\text{x}_1i\) times, and with respect to \(\text{x}_2j\) times, etc. where i = 0, …, order, j = 0, …, order, …. The 0‑th derivative means the function value, thus, deriv[0][0] ... = splineNdInterp.

errEst (Output)

Estimate of the error.

Description¶

The function splineNdInterp interpolates a function of up to 7 variables, defined on a (possibly nonuniform) grid. It fits a polynomial of up to degree 15 in each variable through the grid points nearest the interpolation point. Approximations of partial derivatives are calculated, if requested. If derivatives are desired, high precision is strongly recommended. For more details, see Krogh (1970).

Example¶

The 3D function \(f(x, y, z) = \exp(x + 2y + 3z)\), defined on a 20 by 30 by 40 uniform grid, is interpolated together with several partial derivatives.

from __future__ import print_function
from numpy import *
import math
from pyimsl.math.splineNdInterp import splineNdInterp


N = 3
ND1 = 20
ND2 = 30
ND3 = 40
NDERS = 1

nders = NDERS
d = [ND1, ND2, ND3]
ndeg = [8, 7, 9]

# 20 by 30 by 40 uniform grid used for
# interpolation of F(x,y,z) = exp(x+2*y+3*z)

xdata = zeros((N, ND3), dtype=double)
fdata = zeros((ND1, ND2, ND3), dtype=double)

xdata[0, 0:ND1] = [i * 0.05 for i in range(ND1)]
xdata[1, 0:ND2] = [i * 0.03 for i in range(ND2)]
xdata[2, 0:ND3] = [i * 0.025 for i in range(ND3)]

for i in range(ND1):
    for j in range(ND2):
        for k in range(ND3):
            fdata[i, j, k] = math.e**(xdata[0, i]
                                      + 2 * xdata[1, j] + 3 * xdata[2, k])

#  Interpolate at (0.18,0.43,0.35)
x = [0.18, 0.43, 0.35]
derout = []
error = []

yout = splineNdInterp(d, x, xdata, fdata,
                      nDegree=ndeg,
                      order=nders,
                      deriv=derout,
                      errorEst=error)

# Output F,Fx,Fy,Fz,Fxy,Fxz,Fyz,Fxyz at
# interpolation point

xx = x[0]
yy = x[1]
zz = x[2]
order = [' ', ' ', ' ']
print("EST. VALUE = %g, EST. ERROR = %g" % (yout, error[0]))
print("")
print("            Computed Der.  True Der.     Rel. Err")
for k in range(NDERS + 1):
    for j in range(NDERS + 1):
        for i in range(NDERS + 1):
            order[0] = 'x' if i == 1 else ''
            order[1] = 'y' if j == 1 else ''
            order[2] = 'z' if k == 1 else ''
            tr = (2**j) * (3**k) * math.e**(xx + 2 * yy + 3 * zz)
            relerr = (derout[i][j][k] - tr) / tr
            print("F%-3s   %14.6f %14.6f %14.3e" %
                  (''.join(order), derout[i][j][k], tr, relerr))

Output¶

EST. VALUE = 8.08492, EST. ERROR = 5.93863e-12

            Computed Der.  True Der.     Rel. Err
F            8.084915       8.084915     -1.569e-13
Fx           8.084915       8.084915     -1.799e-13
Fy          16.169830      16.169830     -4.571e-12
Fxy         16.169830      16.169830     -4.400e-12
Fz          24.254745      24.254745     -2.153e-13
Fxz         24.254745      24.254745     -3.329e-13
Fyz         48.509491      48.509491     -4.693e-12
Fxyz        48.509491      48.509491     -5.920e-12

Warning Errors¶

`IMSL_ARG_TOO_BIG`	“`order`” is too large, it has been reset to max(7‑`n`,2).
`IMSL_INTERP_OUTSIDE_DOMAIN`	The interpolation point is outside the domain of the table, so extrapolation is used.

Fatal Errors¶

`IMSL_TOO_MANY_DERIVATIVES`	Too many derivatives requested for the polynomial degree used.
`IMSL_POLY_DEGREE_TOO_LARGE`	One of the polynomial degrees requested is too large for the number of gridlines in that direction.