Cubature (cubature)

There is one function that embodies all functionalities offered by Cubature.

One of the nice things about Cubature is that you can perform multi-dimensional integrations at once, which can be achieved by defining:

• ndim

• xmin

• xmax

There are two types of refinement in the integration methods, both necessary to achieve a desired integration error threshold:

• h additional integration points are added

• p the order of the integration polynomials is increased

It allows the evaluation of vectorized functions, making it convenient to take advantage of NumPy’s speed (see examples with vectorized=True).

See the detailed description below on how to use a vector valued fuction (function that returs an array) or a scalar function.

cubature.cubature(func, ndim, fdim, xmin, xmax, args=(), kwargs={}, abserr=1e-08, relerr=1e-08, norm=0, maxEval=0, adaptive='h', vectorized=False)

Numerical-integration using the cubature method.

Parameters:

funccallable

If vectorized=False the callable must have the form:

f(x_array, *args, **kwargs)

where:

• x_array in an array containing the ndim variables

  being integrated

• args is a tuple containing any other arguments required by the function

• kwargs is a dict containing any keyword arguments required by the function

• the function must return a 1-D np.ndarray object

• example 1, vector valued function:

  fdim = 3
  def func(x_array, *args, **kwargs):
      # note that here ndim=2 (2 variables)
      x, y = x_array
      return np.array([x**2-y**2, x*y, x*y**2])
  

• example 2, scalar function:

  fdim = 1
  def func(x_array, *args, **kwargs):
      # note that here ndim=2 (2 variables)
      x, y = x_array
      return x**2-y**2
  

If vectorized=True the function must have the form:

f(x_array, *args, **kwargs)

where:

• x_array has shape=(npt, ndim)

• args is a tuple containing any other arguments required by the function

• kwargs is a dict containing any keyword arguments required by the function

• example 1, vector valued function:

  # function that returns a vector with 3 values
  fdim = 3
  def func(x_array, *args, **kwargs):
      # note that here ndim=2 (2 variables)
      x = x_array[:, 0]
      y = x_array[:, 1]
      npt = x_array.shape[0]
      out = np.zeros((npt, fdim))
      out[:, 0] = x**2 - y**2
      out[:, 1] = x*y
      out[:, 2] = x*y**2
      return out
  

• example 2, scalar function:

  fdim = 1
  def func(x_array, *args, **kwargs):
      # note that here ndim=3 (3 variables)
      x = x_array[:, 0]
      y = x_array[:, 1]
      z = x_array[:, 2]
      return x**2 + y**2 + z**2
  

The results from both vectorized and non-vectorized examples above should be the same, but the vectorized implementation is much faster since it will take advantage of NumPy’s vectorization capabilities.

ndiminteger

Number dimensions or number of variables being integrated.

fdiminteger

Length of the output vector given by func. It should be 1 if the function returns a scalar.

xmin, xmaxarray-like

These are 1-D arrays with the minimum and maximum integration limits for each variable being integrated. In cubature, if you are using vectorized=True, it is also possible to have one integration limit per value in the vector-valued function.

If vectorized=False:

• xmin, xmax are vectors with the integration limits for each of the ndim dimensions, thus we must assert that xmin.shape[0] == ndim.

If vectorized=True:

• xmin, xmax are vectors with the integration limits for each of the ndim dimensions. With vectorized=True it is possible to use one integration limit per value in the vector-valued function. In this case, we must assert that xmin.shape[0] == ndim*fdim. Here, the limits are read in the order xmin[i*ndim + j], meaning the j-th dimension of the i-th element of the vector-valued function.

argstuple or list, optional

Contains the extra arguments required by func.

kwargsdict-like, optional

Contains the extra keyword arguments passed to func.

adaptivestring, optional

The adaptive scheme used along the adaptive integration:

• ‘h’ means ‘h-adaptive’, where the domain is partitioned

• ‘p’ means ‘p-adaptive’, where the order of the integration rule is increased

The ‘p-adaptive’ scheme is often better for smoth functions in low dimensions.

abserrdouble, optional

Integration stops when estimated absolute error is below this threshold

relerrdouble, optional

Integration stops when estimated error in integral value is below this threshold

norminteger, optional

Specifies the norm that is used to measure the error and determine convergence properties (irrelevant for single-valued functions). The norm argument takes one of the values:

• 0 for ERROR_INDIVIDUAL: convergence is achieved only when each

  integrand individually satisfies the requested error tolerances;

• 1 for ERROR_PAIRED: like ERROR_INDIVIDUAL, except that the

  integrands are grouped into consecutive pairs, with the error tolerance applied in a L2 sense to each pair. This option is mainly useful for integrating vectors of complex numbers, where each consecutive pair of real integrans is the real and imaginary parts of a single complex integrand, and you only care about the error in the complex plane rather than the error in the real and imaginary parts separately;

• 2 for ERROR_L2;

• 3 for ERROR_L1;

• 4 for ERROR_LINF,

  the absolute error is measured as |e| and the relative error as |err|/|val|, where |...| is the L1, L2, or L-infinity norm, respectively. (|x| in the L1 norm is the sum of the absolute values of the components, in the L2 norm is the root mean square of the components, and in the L-infinity norm is the maximum absolute value of the components).

maxEvalinteger, optional

The maximum number of function evaluations.

vectorizedboolean, optional

If vectorized=True the integration points are passed to the integrand function as an array of points, allowing parallel evaluation of different points.

Returns:

valnumpy.ndarray

The 1-D array of length fdim with the computed integral values

errnumpy.ndarray

The 1-D array of length fdim with the estimated errors. For smooth functions this estimate is usually conservative (see the results from the test_cubature.py script.

Notes

• The supplied function must return a 1-D np.ndarray object, even for single-valued functions

References

[1]

Cubature (Multi-dimensional integration).

Examples

>>> # Volume of a sphere:
>>> import numpy as np
>>> from cubature import cubature

>>> def integrand_sphere(x_array, *args):
>>>     r, theta, phi = x_array
>>>     return np.array([r**2*sin(phi)])

>>> ndim = 3
>>> fdim = 1
>>> radius = 1.
>>> xmin = np.array([0, 0, 0])
>>> xmax = np.array([radius, 2*pi, pi])
>>> val, err = cubature(integrand_sphere, ndim, fdim, xmin, xmax)

More Examples

import numpy as np
from numpy import pi, sin

from cubature import cubature

def integrand_brick(x_array):
    return 1.

def integrand_sphere(x_array):
    r, theta, phi = x_array
    return r**2*sin(phi)

def integrand_ellipsoid(x_array, a, b, c):
    rho, phi, theta = x_array
    return a*b*c*rho**2*sin(theta)

def integrand_ellipsoid_v(x_array, a, b, c):
    rho = np.array(x_array[:, 0])
    phi = np.array(x_array[:, 1])
    theta = np.array(x_array[:, 2])
    return a*b*c*rho**2*sin(theta)

def exact_brick(a, b, c):
    return a*b*c

def exact_sphere(r):
    return 4./3*pi*r**3

def exact_ellipsoid(a,b,c):
    return 4./3*pi*a*b*c

if __name__ == '__main__':
    # brick
    print('_________________')
    print('')
    print('Brick')
    a, b, c = 1., 2., 3.
    xmin = np.zeros((3,), dtype=float)
    xmax = np.array([a, b, c], dtype=float)
    val, err = cubature(integrand_brick, 3, 1, xmin, xmax)
    print('Approximated: {0}'.format(val))
    print('Exact: {0}'.format(exact_brick(a,b,c)))
    # sphere
    print('_________________')
    print('')
    print('Sphere')
    radius = 1.
    xmin = np.array([0, 0, 0], np.float64)
    xmax = np.array([radius, 2*pi, pi], np.float64)
    val, err = cubature(integrand_sphere, 3, 1, xmin, xmax)
    print('Approximated: {0}'.format(val))
    print('Exact: {0}'.format(exact_sphere(radius)))
    # ellipsoid
    print('_________________')
    print('')
    print('Ellipsoid')
    a, b, c = 1., 2., 3.
    xmin = np.array([0, 0, 0], np.float64)
    xmax = np.array([1., 2*pi, pi], np.float64)
    val, err = cubature(integrand_ellipsoid, 3, 1, xmin, xmax, args=(a,b,c))
    print('Approximated: {0}'.format(val))
    print('Exact: {0}'.format(exact_ellipsoid(a,b,c)))
    print('_________________')
    # ellipsoid vectorized
    print('_________________')
    print('')
    print('Ellipsoid Vectorized')
    a, b, c = 1., 2., 3.
    xmin = np.array([0, 0, 0], np.float64)
    xmax = np.array([1., 2*pi, pi], np.float64)
    val, err = cubature(integrand_ellipsoid_v, 3, 1, xmin, xmax, args=(a,b,c),
        vectorized=True)
    print('Approximated: {0}'.format(val))
    print('Exact: {0}'.format(exact_ellipsoid(a,b,c)))
    print('_________________')

import numpy as np
from numpy import pi, sin

from cubature import cubature

def integrand_rectangle(x_array):
    return 1.

def integrand_rectangle_v(x_array):
    return np.ones_like(x_array[:, 0])

def integrand_circle(x_array):
    return x_array[0]

def integrand_circle_v(x_array):
    return x_array[:, 0]

def exact_rectangle(a, b):
    return a*b

def exact_circle(r):
    return pi*r**2

if __name__ == '__main__':
    # rectangle
    print('_________________')
    print('')
    print('Rectangle')
    a, b = 3, 5
    xmin = [0, 0]
    xmax = [a, b]
    val, err = cubature(integrand_rectangle, 2, 1, xmin, xmax)
    print('Approximated: {0}'.format(val))
    print('Exact: {0}'.format(exact_rectangle(a, b)))
    # rectangle (vectorized)
    print('_________________')
    print('')
    print('Rectangle (vectorized)')
    a, b = 3, 5
    xmin = [0, 0]
    xmax = [a, b]
    val, err = cubature(integrand_rectangle_v, 2, 1, xmin, xmax, vectorized=True)
    print('Approximated: {0}'.format(val))
    print('Exact: {0}'.format(exact_rectangle(a, b)))
    # circle
    print('_________________')
    print('')
    print('Circle')
    r = 3.
    xmin = [0, 0]
    xmax = [r, 2*pi]
    val, err = cubature(integrand_circle, 2, 1, xmin, xmax)
    print('Approximated: {0}'.format(val))
    print('Exact: {0}'.format(exact_circle(r)))
    print('_________________')
    # circle (vectorized)
    print('_________________')
    print('')
    print('Circle (vectorized)')
    r = 3.
    xmin = [0, 0]
    xmax = [r, 2*pi]
    val, err = cubature(integrand_circle_v, 2, 1, xmin, xmax, vectorized=True)
    print('Approximated: {0}'.format(val))
    print('Exact: {0}'.format(exact_circle(r)))
    print('_________________')

Indices and tables

• Index

• Module Index

• Search Page