intFcnQmc

Integrates a function on a hyper-rectangle using a quasi-Monte Carlo method.

Synopsis

intFcnQmc (fcn, a, b)

Required Arguments

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

User-supplied function to be integrated.

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: errAbs = 1.0e-4.

errRel, float (Input)

Relative accuracy desired.

Default: errAbs = 1.0e-4.

errEst (Output)

An estimate of the absolute value of the error.

maxEvals, int (Input)

Number of evaluations allowed.

Default: No limit.

base, int (Input)

The value of base used to compute the Faure sequence.

skip, int (Input)

The value of skip used to compute the Faure sequence.

Description

Integration of functions over hypercubes by direct methods, such as intFcnHyperRect, is practical only for fairly low dimensional hypercubes. This is because the amount of work required increases exponential as the dimension increases.

An alternative to direct methods is Monte Carlo, in which the integral is evaluated as the value of the function averaged over a sequence of randomly chosen points. Under mild assumptions on the function, this method will converge like \(1/n^{1/2}\), where n is the number of points at which the function is evaluated.

It is possible to improve on the performance of Monte Carlo by carefully choosing the points at which the function is to be evaluated. Randomly distributed points tend to be non-uniformly distributed. The alternative to at sequence of random points is a low-discrepancy sequence. A low-discrepancy sequence is one that is highly uniform.

This function is based on the low-discrepancy Faure sequence as computed by faureNextPoint (see Chapter 10, “Statistics and Random Number Generation”).

Example

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


def fcn(ndim, x):
    sign = -1.0
    sum = 0.0
    for i in range(0, ndim, 1):
        prod = 1.0
        for j in range(0, i + 1, 1):
            prod = prod * x[j]
        sum = sum + (sign * prod)
        sign = -sign
    return sum


a = zeros(10, dtype=double)
b = zeros(10, dtype=double)
for k in range(0, 10, 1):
    a[k] = 0.0
    b[k] = 1.0
q = intFcnQmc(fcn, a, b)
print("integral=%10.3f" % q)

Output

integral=    -0.333

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 = “#”.