Numerical integration, quadrature for various domains

Project description

Your one-stop shop for numerical integration in Python.

More than 1500 numerical integration schemes for line segments, circles, disks, triangles, quadrilaterals, spheres, balls, tetrahedra, hexahedra, wedges, pyramids, n-spheres, n-balls, n-cubes, n-simplices, and the 1D/2D/3D/nD spaces with weight functions exp(-r) and exp(-r2) for fast integration of real-, complex-, and vector-valued functions.

For example, to numerically integrate any function over any given triangle, install quadpy from the Python Package Index with

pip3 install quadpy --user

and do

import numpy

def f(x):
return numpy.sin(x) * numpy.sin(x)

triangle = numpy.array([[0.0, 0.0], [1.0, 0.0], [0.7, 0.5]])

val = quadpy.triangle.strang_fix_cowper_09().integrate(f, triangle)

This uses Strang's rule of degree 6.

All schemese have

scheme.points
scheme.weights
scheme.degree
scheme.citation

scheme.show()
scheme.integrate(
# ...
)

and many have

scheme.points_symbolic
scheme.weights_symbolic

quadpy is fully vectorized, so if you like to compute the integral of a function on many domains at once, you can provide them all in one integrate() call, e.g.,

# shape (3, 5, 2), i.e., (corners, num_triangles, xy_coords)
triangles = numpy.stack([
[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
[[1.2, 0.6], [1.3, 0.7], [1.4, 0.8]],
[[26.0, 31.0], [24.0, 27.0], [33.0, 28]],
[[0.1, 0.3], [0.4, 0.4], [0.7, 0.1]],
[[8.6, 6.0], [9.4, 5.6], [7.5, 7.4]]
], axis=-2)

The same goes for functions with vectorized output, e.g.,

def f(x):
return [numpy.sin(x), numpy.sin(x)]

More examples under test/examples_test.py.

Read more about the dimensionality of the input/output arrays in the wiki.

Schemes

Line segment • Chebyshev-Gauss (both variants, arbitrary degree)
• Clenshaw-Curtis (after Waldvogel, arbitrary degree)
• Fejér-type-1 (after Waldvogel, arbitrary degree)
• Fejér-type-2 (after Waldvogel, arbitrary degree)
• Gauss-Jacobi
• Gauss-Legendre (via NumPy, arbitrary degree)
• Gauss-Lobatto (arbitrary degree)
• Gauss-Kronrod (after Laurie, arbitrary degree)
• Gauss-Patterson (9 nested schemes up to degree 767)
• Gauss-Radau (arbitrary degree)
• closed Newton-Cotes (arbitrary degree)
• open Newton-Cotes (arbitrary degree)
• tanh-sinh quadrature (see above)

See below for how to generate Gauss formulas for your own weight functions.

Example:

import numpy

scheme.show()
val = scheme.integrate(lambda x: numpy.exp(x), [0.0, 1.0])

1D half-space with weight function exp(-r) • Generalized Gauss-Laguerre

Example:

scheme = quadpy.e1r.gauss_laguerre(5, alpha=0)
scheme.show()
val = scheme.integrate(lambda x: x**2)

1D space with weight function exp(-r2) Example:

scheme.show()
val = scheme.integrate(lambda x: x**2)

Circle • Krylov (1959, arbitrary degree)

Example:

scheme.show()
val = scheme.integrate(lambda x: numpy.exp(x), [0.0, 0.0], 1.0)

Triangle Apart from the classical centroid, vertex, and seven-point schemes we have

Example:

scheme.show()
val = scheme.integrate(lambda x: numpy.exp(x), [[0.0, 0.0], [1.0, 0.0], [0.5, 0.7]])

Disk Example:

import numpy

scheme.show()
val = scheme.integrate(lambda x: numpy.exp(x), [0.0, 0.0], 1.0) Example:

import numpy

val = scheme.integrate(
lambda x: numpy.exp(x),
[[[0.0, 0.0], [1.0, 0.0]], [[0.0, 1.0], [1.0, 1.0]]],
)

The points are specified in an array of shape (2, 2, ...) such that arr is the lower left corner, arr the upper right. If your quadrilateral has its sides aligned with the coordinate axes, you can use the convenience function

to generate the array.

2D space with weight function exp(-r) Example:

scheme.show()
val = scheme.integrate(lambda x: x**2)

2D space with weight function exp(-r2) Example:

scheme.show()
val = scheme.integrate(lambda x: x**2)

Sphere • via Stroud (1971):
• Lebedev (1976, 34 schemes up to degree 131)
• Bažant-Oh (1986, 3 schemes up to degree 11)
• Heo-Xu (2001, 27 schemes up to degree 39, single-precision)
• Fliege-Maier (2007, 4 schemes up to degree 4, single-precision)

Example:

import numpy

scheme.show()
val = scheme.integrate(lambda x: numpy.exp(x), [0.0, 0.0, 0.0], 1.0)

Integration on the sphere can also be done for function defined in spherical coordinates:

import numpy

val = scheme.integrate_spherical(
lambda azimuthal, polar: numpy.sin(azimuthal)**2 * numpy.sin(polar),
)

Ball • Hammer-Stroud (1958, 6 schemes up to degree 7)
• via: Stroud (1971):
• Ditkin (1948, 3 schemes up to degree 7)
• Mysovskih (1964, degree 7)
• 2 schemes up to degree 14

Example:

import numpy

scheme.show()
val = scheme.integrate(
lambda x: numpy.exp(x),
[0.0, 0.0, 0.0], 1.0,
)

Tetrahedron Example:

import numpy

scheme.show()
val = scheme.integrate(
lambda x: numpy.exp(x),
[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 0.7, 0.0], [0.3, 0.9, 1.0]],
)

Hexahedron Example:

import numpy

scheme.show()
val = scheme.integrate(
lambda x: numpy.exp(x),
quadpy.hexahedron.cube_points([0.0, 1.0], [-0.3, 0.4], [1.0, 2.1]),
)

Pyramid • Felippa (9 schemes up to degree 5)

Example:

import numpy

val = scheme.integrate(
lambda x: numpy.exp(x),
[
[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 0.7, 0.0], [0.3, 0.9, 0.0],
[0.0, 0.1, 1.0],
]
)

Wedge Example:

import numpy

lambda x: numpy.exp(x),
[
[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 0.7, 0.0]],
[[0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.5, 0.7, 1.0]],
]
)

3D space with weight function exp(-r) Example:

scheme.show()
val = scheme.integrate(lambda x: x**2)

3D space with weight function exp(-r2) Example:

scheme.show()
val = scheme.integrate(lambda x: x**2)

n-Simplex

Example:

import numpy

scheme = quadpy.nsimplex.GrundmannMoeller(dim, 3)
dim = 4
val = scheme.integrate(
lambda x: numpy.exp(x),
numpy.array([
[0.0, 0.0, 0.0, 0.0],
[1.0, 2.0, 0.0, 0.0],
[0.0, 1.0, 0.0, 0.0],
[0.0, 3.0, 1.0, 0.0],
[0.0, 0.0, 4.0, 1.0],
])
)

n-Sphere

• via Stroud (1971):
• Stroud (1967, degree 7)
• Stroud (1969, 3 <= n <= 16, degree 11)
• 6 schemes up to degree 5
• Dobrodeev (1978, n >= 2, degree 5)

Example:

import numpy

dim = 4
val = scheme.integrate(lambda x: numpy.exp(x), numpy.zeros(dim), 1.0)

n-Ball

• Dobrodeev (1970, n >= 3, degree 7)
• via Stroud (1971):
• Stroud (1957, degree 2)
• Hammer-Stroud (1958, 2 schemes up to degree 5)
• Stroud (1966, 4 schemes of degree 5)
• Stroud (1967, 4 <= n <= 7, 2 schemes of degree 5)
• Stroud (1967, n >= 3, 3 schemes of degree 7)
• Stenger (1967, 6 schemes up to degree 11)
• Dobrodeev (1978, 2 <= n <= 20, degree 5)

Example:

import numpy

dim = 4
val = scheme.integrate(lambda x: numpy.exp(x), numpy.zeros(dim), 1.0)

n-Cube

Example:

import numpy

dim = 4
lambda x: numpy.exp(x),
[0.0, 1.0], [0.1, 0.9], [-1.0, 1.0], [-1.0, -0.5]
)
)

nD space with weight function exp(-r)

Example:

dim = 4
val = scheme.integrate(lambda x: x**2)

nD space with weight function exp(-r2)

• via Stroud (1971):
• Stroud-Secrest (1963, 4 schemes up to degree 5)
• Stroud (1967, 2 schemes of degree 5)
• Stroud (1967, 3 schemes of degree 7)
• Stenger (1971, 6 schemes up to degree 11, varying dimensionality restrictions)
• 5 schemes up to degree 5

Example:

dim = 4
val = scheme.integrate(lambda x: x**2)

Installation

quadpy is available from the Python Package Index, so with

pip3 install quadpy --user

you can install.

Testing

To run the tests, just check out this repository and type

MPLBACKEND=Agg pytest

quadpy is published under the MIT license.

Project details

This version 0.13.2 0.13.1 0.13.0 0.12.11 0.12.10 0.12.9 0.12.8 0.12.7 0.12.6 0.12.5 0.12.4 0.12.3 0.12.2 0.12.1 0.12.0 0.11.4 0.11.3 0.11.2 0.11.1 0.11.0 0.10.7 0.10.6 0.10.5 0.10.4 0.10.3 0.10.2 0.10.1 0.10.0 0.9.4 0.9.3 0.9.2 0.9.0 0.8.3 0.8.2 0.8.1 0.8.0 0.7.1 0.7.0 0.6.1 0.6.0 0.5.5 0.5.4 0.5.3 0.5.2 0.5.1 0.5.0 0.4.2 0.4.1 0.4.0