intFcnHyperRect

Integrate a function on a hyper-rectangle,

\[\int_{a_0}^{b_0} \ldots \int_{a_{n-1}}^{b_{n-1}} f\left(x_0, \ldots, x_{n-1}\right) dx_{n-1} \ldots dx_0\]

Synopsis

intFcnHyperRect (fcn, ndim, a, b)

Required Arguments

float fcn (ndim, x[]) (Input)

User-supplied function to be integrated.

int ndim (Input)

The dimension of the hyper-rectangle.

float a[] (Input)

Lower limits of integration.

float b[] (Input)

Upper limits of integration.

Return Value

The value of

\[\int_{a_0}^{b_0} \ldots \int_{a_{n-1}}^{b_{n-1}} f\left(x_0, \ldots, x_{n-1}\right) dx_{n-1} \ldots dx_0\]

is returned. If no value can be computed, then NaN is returned.

Optional Arguments

errAbs, float (Input)

Absolute accuracy desired.

Default: \(\mathit{errAbs} = \sqrt{\varepsilon}\) where ɛ is the machine precision

errRel, float (Input)

Relative accuracy desired.

Default: \(\mathit{errRel} = \sqrt{\varepsilon}\) where ɛ is the machine precision

errEst (Output)

An estimate of the absolute value of the error.

maxEvals, int (Input)

Number of evaluations allowed.

Default: maxEvals = \(32^n\).

Description

The function intFcnHyperRect approximates the n-dimensional iterated integral

\[\int_{a_0}^{b} \ldots \int_{a_{n-1}}^{b_{n-1}} f\left(x_0, \ldots, x_{n-1}\right) dx_{n-1} \ldots dx_0\]

An estimate of the error is returned in the optional argument errEst. The approximation is achieved by iterated applications of product Gauss formulas. The integral is first estimated by a two-point tensor product formula in each direction. Then for \(i = 1, ..., n\), the function calculates a new estimate by doubling the number of points in the i-th direction, then halving the number immediately afterwards if the new estimate does not change appreciably. This process is repeated until either one complete sweep results in no increase in the number of sample points in any dimension; the number of Gauss points in one direction exceeds 256; or the number of function evaluations needed to complete a sweep exceeds maxEvals.

Example

In this example, we compute the integral of

\[e^{-\left(x_1^2 + x_2^2 + x_3^2\right)}\]

on an expanding cube. The values of the error estimates are machine dependent. The exact integral over \(R^3\) is \(\pi^{3/2}\).

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


def fcn(n, x):
    s = x[0] * x[0] + x[1] * x[1] + x[2] * x[2]
    return exp(-s)


print("       integral       limit")
a = zeros(3, dtype=double)
b = zeros(3, dtype=double)
limit = pow(constant("pi"), 1.5)

# Evaluate the integral
for i in range(0, 6):
    for j in range(0, 3):
        a[j] = (-(i + 1) / 2.)
        b[j] = (i + 1) / 2.
        q = intFcnHyperRect(fcn, a, b)
    # Print the result and the limiting answer
    print("   %10.3f    %10.3f" % (q, limit))

Output

       integral       limit
        0.785         5.568
        3.332         5.568
        5.021         5.568
        5.491         5.568
        5.562         5.568
        5.568         5.568

Warning Errors

IMSL_MAX_EVALS_TOO_LARGE

The argument maxEvals was set greater than \(2^{8n}\).

Fatal Errors

IMSL_NOT_CONVERGENT	The maximum number of function evaluations has been reached, and convergence has not been attained
IMSL_STOP_USER_FCN	Request from user supplied function to stop algorithm. User flag = “#”.