scattered2dInterp

Computes a smooth bivariate interpolant to scattered data that is locally a quintic polynomial in two variables.

Synopsis

scattered2dInterp (xydata, fdata, xOut, yOut)

Required Arguments

float xydata[[]] (Input)

Array with ndata\*2 components containing the data points for the interpolation problem. The i-th data point \((x_i, y_i)\) is stored consecutively in the \(2i\) and \(2i + 1\) positions of xydata.

float fdata[] (Input)

Array of size ndata containing the values to be interpolated.

float xOut[] (Input)

Array of length nxOut specifying the x values for the output grid. It must be strictly increasing.

float yOut[] (Input)

Array of length nyOut specifying the y values for the output grid. It must be strictly increasing.s

Return Value

The nxOut × nyOut grid of values of the interpolant. If no answer can be computed, then None is returned.

Description

The function scattered2dInterp computes a \(C^1\) interpolant to scattered data in the plane. Given the data points

\[\left\{\left(x_i,y_i,f_i\right)\right\}_{i=0}^{n-1}\]

in \(R^3\) where n = ndata, scattered2dInterp returns the values of the interpolant s on the user-specified grid. The computation of s is as follows: First the Delaunay triangulation of the points

\[\left\{\left(x_i,y_i\right)\right\}_{i=0}^{n-1}\]

is computed. On each triangle T in this triangulation, s has the form

\[s(x,y) = \sum_{m+n \leq 5} c_{mn}^{T} x^m y^n \phantom{...} \forall x,y \in T\]

Thus, s is a bivariate quintic polynomial on each triangle of the triangulation. In addition, we have

\[s\left(x_i,y_i\right) = f_i \text{ for } i = 0, \ldots, n-1\]

and s is continuously differentiable across the boundaries of neighboring triangles. These conditions do not exhaust the freedom implied by the above representation. This additional freedom is exploited in an attempt to produce an interpolant that is faithful to the global shape properties implied by the data. For more information on this procedure, refer to the article by Akima (1978). The output grid is specified by the two integer variables nxOut and nyOut that represent the number of grid points in the first (second) variable and by two real vectors that represent the first (second) coordinates of the grid.

Examples

Example 1

In this example, the interpolant to the linear function (3 + 7x + 2y) is computed from 20 data points equally spaced on the circle of radius 3. Then the values are printed on a 3 ×3 grid.

from __future__ import print_function
from numpy import *
from pyimsl.math.scattered2dInterp import scattered2dInterp
from pyimsl.math.constant import constant

# Define function


def F(x, y):
    return 3. + 7. * x + 2. * y


ndata = 20
outdata = 3
fdata = empty(ndata)
xydata = empty((ndata, 2))
x_out = empty(outdata)
y_out = empty(outdata)
pi = constant("pi")

# Set up output grid
for i in range(0, outdata):
    x_out[i] = y_out[i] = float(i) / (outdata - 1)
for i in range(0, ndata):
    xydata[i, 0] = 3. * cos(pi * i * 2 / ndata)
    xydata[i, 1] = 3. * sin(pi * i * 2 / ndata)
    fdata[i] = F(xydata[i, 0], xydata[i, 1])

# Compute scattered data interpolant
surf = scattered2dInterp(xydata, fdata, x_out, y_out)

# Print results
print("       x       y      F(x, y)    Interpolant    Error")
for i in range(0, outdata):
    for j in range(0, outdata):
        x = x_out[i]
        y = y_out[j]
        z = surf[i, j]
        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
   0.000   0.000       3.000        3.000       0.0000
   0.000   0.500       4.000        4.000       0.0000
   0.000   1.000       5.000        5.000       0.0000
   0.500   0.000       6.500        6.500       0.0000
   0.500   0.500       7.500        7.500       0.0000
   0.500   1.000       8.500        8.500       0.0000
   1.000   0.000      10.000       10.000       0.0000
   1.000   0.500      11.000       11.000       0.0000
   1.000   1.000      12.000       12.000       0.0000

Example 2

Recall that in the first example, the interpolant to the linear function 3 + 7x + 2y is computed from 20 data points equally spaced on the circle of radius 3. We then print the values on a 3 × 3 grid. This example used the optional arguments to indicate that the answer is stored noncontiguously in a two-dimensional array surf with column dimension equal to 11.

from __future__ import print_function
from numpy import *
from pyimsl.math.constant import constant
from pyimsl.math.scattered2dInterp import scattered2dInterp

# Define function


def F(x, y):
    return 3. + 7. * x + 2. * y


ndata = 20
outdata = 3
fdata = empty(ndata)
xydata = empty((ndata, 2))
x_out = empty(outdata)
y_out = empty(outdata)
pi = constant("pi")

# Set up output grid
for i in range(0, outdata):
    x_out[i] = y_out[i] = float(i) / (outdata - 1)
for i in range(0, ndata):
    xydata[i, 0] = 3. * cos(pi * i * 2 / ndata)
    xydata[i, 1] = 3. * sin(pi * i * 2 / ndata)
    fdata[i] = F(xydata[i, 0], xydata[i, 1])

# Compute scattered data interpolant
surf = scattered2dInterp(xydata, fdata, x_out, y_out)

# Print results
print("       x       y      F(x, y)    Interpolant    Error")
for i in range(0, outdata):
    for j in range(0, outdata):
        x = x_out[i]
        y = y_out[j]
        z = surf[i, j]
        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
   0.000   0.000       3.000        3.000       0.0000
   0.000   0.500       4.000        4.000       0.0000
   0.000   1.000       5.000        5.000       0.0000
   0.500   0.000       6.500        6.500       0.0000
   0.500   0.500       7.500        7.500       0.0000
   0.500   1.000       8.500        8.500       0.0000
   1.000   0.000      10.000       10.000       0.0000
   1.000   0.500      11.000       11.000       0.0000
   1.000   1.000      12.000       12.000       0.0000

Fatal Errors

IMSL_DUPLICATE_XYDATA_VALUES	The two-dimensional data values must be distinct.
IMSL_XOUT_NOT_STRICTLY_INCRSING	The vector xOut must be strictly increasing.
IMSL_YOUT_NOT_STRICTLY_INCRSING	The vector yOut must be strictly increasing.