spline2dIntegral¶

Evaluates the integral of a tensor-product spline on a rectangular domain.

Synopsis¶

spline2dIntegral (a, b, c, d, sp)

Required Arguments¶

float a (Input)

float b (Input): The integration limits for the first variable of the tensor-product spline.

float c (Input)

float d (Input): The integration limits for the second variable of the tensor-product spline.
spline sp (Input): The structure that represents the spline.

Return Value¶

The value of the integral of the tensor-product spline over the rectangle [a, b] × [c, d]. If no value can be computed, NaN is returned.

Description¶

The function spline2dIntegral computes the integral of a tensor-product spline. If s is the spline, then this function returns

\[\int_a^b \int_c^d s(x,y)dydx\]

This function uses the (univariate integration) identity (22) in de Boor (1978, p. 151)

\[\int_{t_0}^{x} \sum_{i=0}^{n-1} \alpha_i B_{i,k}(\tau)d\tau = \sum_{i=0}^{r-1} \left[ \sum_{j=0}^{i} \alpha_j \frac{t_{j+k} - t_j}{k} \right] B_{i,k+1}(x)\]

where \(t_0 \leq x \leq t_r\).

It assumes (for all knot sequences) that the first and last k knots are stacked, that is, \(t_0 = ... = t_{k-1}\) and \(t_n = ... = t_{n+k-1}\), where k is the order of the spline in the x or y direction.

Example¶

This example integrates a two-dimensional, tensor-product spline over the rectangle [0, x] × [0, y].

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

# Define function


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

# The integral of F from 0 to x and 0 to y


def FI(x, y):
    return y * x * x * x * x / 4. + x * y * y * y / 3.


# Set up grid
ndata = 11
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     FI(x, y)     Integral       Error")
for i in range(0, 2):
    x = float(i + 1) / 3.0
    for j in range(0, 2):
        y = float(j + 1) / 3.0
        z = spline2dIntegral(0.0, x, 0.0, y, sp)
        print("  %6.3f  %6.3f  %10.3f   %10.3f   %10.4f"
              % (x, y, FI(x, y), z, abs(FI(x, y) - z)))

Output¶

       x       y     FI(x, y)     Integral       Error
333   0.333       0.005        0.005       0.0000
333   0.667       0.035        0.035       0.0000
667   0.333       0.025        0.025       0.0000
667   0.667       0.099        0.099       0.0000

Warning Errors¶

`IMSL_SPLINE_LEFT_ENDPT`	The left endpoint of X integration is not within the knot sequence. Integration occurs only from \(t_{order-1}\) to `b`.
`IMSL_SPLINE_RIGHT_ENDPT`	The right endpoint of X integration is not within the knot sequence. Integration occurs only from \(t_{order-1}\) to `a`.
`IMSL_SPLINE_LEFT_ENDPT_1`	The left endpoint of X integration is not within the knot sequence. Integration occurs only from `b` to \(t_{splineSpaceDim‑1}\).
`IMSL_SPLINE_RIGHT_ENDPT_1`	The right endpoint of X integration is not within the knot sequence. Integration occurs only from `a` to \(t_{splineSpaceDim‑1}\).
`IMSL_SPLINE_LEFT_ENDPT_2`	The left endpoint of Y integration is not within the knot sequence. Integration occurs only from \(t_{order-1}\) to `d`.
`IMSL_SPLINE_RIGHT_ENDPT_2`	The right endpoint of Y integration is not within the knot sequence. Integration occurs only from \(t_{order-1}\) to `c`.
`IMSL_SPLINE_LEFT_ENDPT_3`	The left endpoint of Y integration is not within the knot sequence. Integration occurs only from `d` to \(t_{splineSpaceDim‑1}\).
`IMSL_SPLINE_RIGHT_ENDPT_3`	The right endpoint of Y integration is not within the knot sequence. Integration occurs only from `c` to \(t_{splineSpaceDim‑1}\).

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.