naginterfaces.library.quad.md_numth¶

naginterfaces.library.quad.md_numth(f, region, npts, vk, nrand, itrans=0, data=None)[source]¶

md_numth calculates an approximation to a definite integral in up to $20$ dimensions, using the Korobov–Conroy number theoretic method.

For full information please refer to the NAG Library document for d01gc

https://www.nag.com/numeric/nl/nagdoc_29.3/flhtml/d01/d01gcf.html

Parameters

fcallable retval = f(x, data=None)

$f$ must return the value of the integrand $f$ at a given point.

Parameters

xfloat, ndarray, shape $(ndim)$: The coordinates of the point at which the integrand $f$ must be evaluated.
dataarbitrary, optional, modifiable in place: User-communication data for callback functions.

Returns

retvalfloat: The value of $f (x)$ evaluated at $x$ .

regioncallable (c, d) = region(x, j, data=None)

$r e g i o n$ must evaluate the limits of integration in any dimension.

Parameters

xfloat, ndarray, shape $(ndim)$: $x [0], \dots, x [j - 2]$ contain the current values of the first $(j - 1)$ variables, which may be used if necessary in calculating $c_{j}$ and $d_{j}$ .
jint: The index $j$ for which the limits of the range of integration are required.
dataarbitrary, optional, modifiable in place: User-communication data for callback functions.

Returns

cfloat: The lower limit $c_{j}$ of the range of $x_{j}$ .
dfloat: The upper limit $d_{j}$ of the range of $x_{j}$ .

nptsint

The Korobov rule to be used. There are two alternatives depending on the value of $n p t s$ .

$1 \leq n p t s \leq 6$ .

In this case one of six preset rules is chosen using $2129$ , $5003$ , $10007$ , $20011$ , $40009$ or $80021$ points depending on the respective value of $n p t s$ being $1$ , $2$ , $3$ , $4$ , $5$ or $6$ .
$n p t s > 6$ .

$n p t s$ is the number of actual points to be used with corresponding optimal coefficients supplied in the array $v k$ .

vkfloat, array-like, shape $(ndim)$

If $n p t s > 6$ , $v k$ must contain the $n$ optimal coefficients (which may be calculated using md_numth_coeff_prime() or md_numth_coeff_2prime()).

If $n p t s \leq 6$ , $v k$ need not be set.

nrandint

The number of random samples to be generated in the error estimation (generally a small value, say $3$ to $5$ , is sufficient). The total number of integrand evaluations will be $n r a n d \times n p t s$ .

itransint, optional

Indicates whether the periodising transformation is to be used.

$i t r a n s ='0'$

The transformation is to be used.

$i t r a n s \neq'0'$

The transformation is to be suppressed (to cover cases where the integrand may already be periodic or where you want to specify a particular transformation in the definition of $f$ ).

dataarbitrary, optional

User-communication data for callback functions.

Returns

vkfloat, ndarray, shape $(ndim)$

If $n p t s > 6$ , $v k$ is unchanged.

If $n p t s \leq 6$ , $v k$ contains the $n$ optimal coefficients used by the preset rule.

resfloat

The approximation to the integral $I$ .

errfloat

The standard error as computed from $n r a n d$ sample values. If $n r a n d = 1$ , $e r r$ contains zero.

Raises

NagValueError

(errno $1$ )

On entry, $ndim = ⟨ v a l u e ⟩$ .

Constraint: $1 \leq ndim \leq 20$ .

(errno $2$ )

On entry, $n p t s = ⟨ v a l u e ⟩$ .

Constraint: $n p t s \geq 1$ .

(errno $3$ )

On entry, $n r a n d = ⟨ v a l u e ⟩$ .

Constraint: $n r a n d \geq 1$ .

Notes

No equivalent traditional C interface for this routine exists in the NAG Library.

md_numth calculates an approximation to the integral

I = \int_{c_{1}}^{d_{1}} d x_{1}, \dots, \int_{c_{n}}^{d_{n}} d x_{n} f (x_{1}, x_{2}, \dots, x_{n})

using the Korobov–Conroy number theoretic method (see Korobov (1957), Korobov (1963) and Conroy (1967)). The region of integration defined in (1) is such that generally $c_{i}$ and $d_{i}$ may be functions of $x_{1}, x_{2}, \dots, x_{i - 1}$ , for $i = 2, 3, \dots, n$ , with $c_{1}$ and $d_{1}$ constants. The integral is first of all transformed to an integral over the $n$ -cube ${[0, 1]}^{n}$ by the change of variables

x_{i} = c_{i} + (d_{i} - c_{i}) y_{i}, i = 1, 2, \dots, n .

The method then uses as its basis the number theoretic formula for the $n$ -cube, ${[0, 1]}^{n}$ :

\int_{0}^{1} d x_{1} \dots \int_{0}^{1} d x_{n} g (x_{1}, x_{2}, \dots, x_{n}) = \frac{1}{p} p \sum k = 1 g ({k \frac{a_{1}}{p}}, \dots, {k \frac{a_{n}}{p}}) - E

where ${x}$ denotes the fractional part of $x$ , $a_{1}, a_{2}, \dots, a_{n}$ are the so-called optimal coefficients, $E$ is the error, and $p$ is a prime integer. (It is strictly only necessary that $p$ be relatively prime to all $a_{1}, a_{2}, \dots, a_{n}$ and is in fact chosen to be even for some cases in Conroy (1967).) The method makes use of properties of the Fourier expansion of $g (x_{1}, x_{2}, \dots, x_{n})$ which is assumed to have some degree of periodicity. Depending on the choice of $a_{1}, a_{2}, \dots, a_{n}$ the contributions from certain groups of Fourier coefficients are eliminated from the error, $E$ . Korobov shows that $a_{1}, a_{2}, \dots, a_{n}$ can be chosen so that the error satisfies

E \leq C K p^{- α} {l n}^{α β} (p)

where $α$ and $C$ are real numbers depending on the convergence rate of the Fourier series, $β$ is a constant depending on $n$ , and $K$ is a constant depending on $α$ and $n$ . There are a number of procedures for calculating these optimal coefficients. Korobov imposes the constraint that

a_{1} = 1 and a_{i} = a^{i - 1} (m o d p)

and gives a procedure for calculating the argument, $a$ , to satisfy the optimal conditions.

In this function the periodisation is achieved by the simple transformation

x_{i} = y_{i}^{2} (3 - 2 y_{i}), i = 1, 2, \dots, n .

More sophisticated periodisation procedures are available but in practice the degree of periodisation does not appear to be a critical requirement of the method.

An easily calculable error estimate is not available apart from repetition with an increasing sequence of values of $p$ which can yield erratic results. The difficulties have been studied by Cranley and Patterson (1976) who have proposed a Monte Carlo error estimate arising from converting (2) into a stochastic integration rule by the inclusion of a random origin shift which leaves the form of the error (3) unchanged; i.e., in the formula (2), ${k \frac{a_{i}}{p}}$ is replaced by ${α_{i} + k \frac{a_{i}}{p}}$ , for $i = 1, 2, \dots, n$ , where each $α_{i}$ , is uniformly distributed over $[0, 1]$ . Computing the integral for each of a sequence of random vectors $α$ allows a ‘standard error’ to be estimated.

This function provides built-in sets of optimal coefficients, corresponding to six different values of $p$ . Alternatively, the optimal coefficients may be supplied by you. Functions md_numth_coeff_prime() and md_numth_coeff_2prime() compute the optimal coefficients for the cases where $p$ is a prime number or $p$ is a product of two primes, respectively.

References

Conroy, H, 1967, Molecular Shroedinger equation VIII. A new method for evaluting multi-dimensional integrals, J. Chem. Phys. (47), 5307–5318

Cranley, R and Patterson, T N L, 1976, Randomisation of number theoretic methods for mulitple integration, SIAM J. Numer. Anal. (13), 904–914

Korobov, N M, 1957, The approximate calculation of multiple integrals using number theoretic methods, Dokl. Acad. Nauk SSSR (115), 1062–1065

Korobov, N M, 1963, Number Theoretic Methods in Approximate Analysis, Fizmatgiz, Moscow

NAG and Python

Return to Front

naginterfaces.library.quad.md_numth¶

naginterfaces.library.quad.md_​numth¶

naginterfaces.library.quad.md_numth¶