mmf.math.integrate.multidimen

dcuhre(f, a, b[, epsabs, epsrel, ...])

Description:

Multi-dimensional integration.

Integrates a function over a hypercube in multiple dimensions.

mmf.math.integrate.multidimen.dcuhre(f, a, b, epsabs=1.49e-08, epsrel=1.49e-08, full_output=False, minpts=0, maxpts=100000, key=0, nproc=1, nw=None)

Description:

Return the integral of f over the hypercube a <= x <= b.

Inputs:

f – f(x) should return a one dimensional vector (ndarray for

performance) of the function(s) evaluated at the point x.

a – Lower limits of integration (2 <= ndim = len(a) <= 15) b – Upper limits of integration (2 <= ndim = len(b) <= 15) epsabs – Desired absolute tolerance epsabs – Desired relative tolerance minpts – Minimum number of points to sample function at maxpts – Maximum number of function evaluations key – Key to selected local integration rule.

key = 0 is the default value.

For ndim = 2 the degree 13 rule is selected. For ndim = 3 the degree 11 rule is selected. For ndim > 3 the degree 9 rule is selected.

key = 1 gives the user the 2 dimensional degree 13

integration rule that uses 65 evaluation points.

key = 2 gives the user the 3 dimensional degree 11

integration rule that uses 127 evaluation points.

key = 3 gives the user the degree 9 integration rule. key = 4 gives the user the degree 7 integration rule.

This is the recommended rule for problems that require great adaptivity.

nproc – Number of processors (if DCUHRE compiled with

multiple processor support)

nw – Defines the length of the working array. Should be

automatically determined (see below).

full_output – non-zero to return a dictionary of integration

information. If non-zero, warning messages are also suppressed and the message is appended to the output tuple.

The result hopefully satisfies for each component k of f the following claim for accuracy:

abs(ans - exact) <= max(epsabs,epsrel*abs(ans))

Based on the FORTRAN code DCUHRE which is a driver for the integration routine DADHRE, which repeatedly subdivides the region of integration and estimates the integrals and the errors over the subregions with greatest estimated errors until the error request is met or maxpts function evaluations have been used.

Let ndim = len(a) == len(b) be the dimension of the domain hypercube:

• For ndim = 2 the default integration rule is of degree 13 and uses 65 evaluation points.

• For ndim = 3 the default integration rule is of degree 11 and uses 127 evaluation points.

• For ndim greater then 3 the default integration rule is of degree 9 and uses num evaluation points where num = 1 + 4*2*ndim + 2*ndim*(ndim-1) + 4*ndim*(ndim-1) +

  4*ndim*(ndim-1)*(ndim-2)/3 + 2**ndim

The degree 9 rule may also be applied for ndim = 2 and ndim = 3. A rule of degree 7 is available in all dimensions. The number of evaluation points used by the degree 7 rule is

num = 1 + 3*2*ndim + 2*ndim*(ndim-1) + 2**ndim

The number of function evaluations over each subregion is num: if (key == 0 or key == 1) and (ndim == 2):

num = 65

elif (key == 0 or key == 2) and (ndim == 3):

num = 127

elif (key == 0 and ndim > 3) or (key == 3):

num = 1 + 4*2*ndim + 2*ndim*(ndim-1) + 4*ndim*(ndim-1) + 4*ndim*(ndim-1)*(ndim-2)/3 + 2**ndim

elif (key == 4):

num = 1 + 3*2*ndim + 2*ndim*(ndim-1) + 2**ndim

maxpts >= 3*num and maxpts >= minpts

for 3 < ndim < 13 the minimum values for maxpts are: ndim = 4 5 6 7 8 9 10 11 12 key = 3: 459 819 1359 2151 3315 5067 7815 12351 20235 key = 4: 195 309 483 765 1251 2133 3795 7005 13299

Let maxsub denote the maximum allowed number of subregions for the given values of maxpts, key and ndim.

maxsub = (maxpts-num)/(2*num) + 1 nw >= maxsub*(2*ndim + 2*numfun + 2) + 17*numfun + 1

For efficient execution on parallel computers nw should be chosen greater or equal to maxsub*(2*ndim+2*numfun+2)+17*numfun*mdiv+1 where mdiv is the number of subregions that are divided in each subdivision step. mdiv is default set internally in DCUHRE equal to 1. For efficient execution on parallel computers with nproc processors mdiv should be set equal to the smallest integer such that (2*mdiv)%nproc == 0.

When DCUHRE computes estimates to a vector of integrals, all components of the vector are given the same treatment. That is, I(f[j]) and I(f[k]) for j != k, are estimated with the same subdivision of the region of integration. For integrals with enough similarity, we may save time by applying DCUHRE to all integrands in one call. For integrals that vary continuously as functions of some parameter, the estimates produced by DCUHRE will also vary continuously when the same subdivision is applied to all components. This will generally not be the case when the different components are given separate treatment.

On the other hand this feature should be used with caution when the different components of the integrals require clearly different subdivisions.

Credit: Jarle Berntsen, The Computing Centre, University of Bergen, Thormohlens gt. 55, N-5008 Bergen, Norway Phone.. 47-5-544055 Email.. jarle@eik.ii.uib.no Terje O. Espelid, Department of Informatics, University of Bergen, Thormohlens gt. 55, N-5008 Bergen, Norway Phone.. 47-5-544180 Email.. terje@eik.ii.uib.no Alan Genz, Computer Science Department, Washington State University, Pullman, WA 99163-1210, USA Phone.. 509-335-2131 Email.. acg@cs2.cs.wsu.edu

>>> def f(x):
...     ndim = len(x)
...     nf = ndim + 1
...     s = 0
...     for n in xrange(ndim):
...         s = s + (n+1)*x[n]**2
...     fv = numpy.zeros(nf)
...     fv[0] = numpy.exp(-s/2.0)
...     for n in xrange(ndim):
...         fv[n+1] = x[n]*fv[0]
...     return fv
>>> ndim = 5
>>> a = numpy.zeros(ndim,float)
>>> b = numpy.ones(ndim,float)
>>> abseps = 0
>>> releps = 0.001
>>> ans = dcuhre(f=f,a=a,b=b,epsabs=abseps,epsrel=releps)
Success!
>>> def f(x):
...     return sum(x*x)
>>> Ndim = 4.0
>>> abs_tol = 0.0001
>>> rel_tol = 0.001
>>> ans = dcuhre(f,a=numpy.zeros(Ndim),b=numpy.ones(Ndim),                         epsabs=abs_tol,epsrel=rel_tol)
Success!

>>> abs(ans[0] - Ndim/3.0) < abs_tol
True
>>> abs(ans[0] - Ndim/3.0)/Ndim*3.0 < rel_tol
True