Your one-stop shop for numerical integration in Python.
Hundreds of 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 2D/3D/nD spaces with weight functions exp(-r) and exp(-r2).
To numerically integrate any function over any given triangle, do
import numpy
import quadpy
def f(x):
return numpy.sin(x[0]) * numpy.sin(x[1])
triangle = numpy.array([[0.0, 0.0], [1.0, 0.0], [0.7, 0.5]])
val = quadpy.triangle.integrate(f, triangle, quadpy.triangle.Strang(9))
This uses Strang's rule of degree 6.
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.,
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[0]), numpy.sin(x[1])]
More examples under test/examples_test.py.
quadpy can do adaptive quadrature for certain domains. Again, everything is fully vectorized, so you can provide multiple intervals and vector-valued functions.
val, error_estimate = quadpy.line_segment.adaptive_integrate(
lambda x: x * sin(5 * x),
[0.0, pi],
1.0e-10
)
val, error_estimate = quadpy.triangle.adaptive_integrate(
lambda x: x[0] * sin(5 * x[1]),
[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]],
1.0e-10
)
ProTip: You can provide many triangles that together form a domain to get an approximation of the integral over the domain.
- 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-Hermite (via NumPy, arbitrary degree)
- Gauss-Laguerre (via NumPy, arbitrary degree)
- Gauss-Legendre (via NumPy, arbitrary degree)
- Gauss-Lobatto (arbitrary degree)
- Gauss-Kronrod (after Laurie, arbitrary degree)
- Gauss-Patterson (7 schemes up to degree 191)
- Gauss-Radau (arbitrary degree)
- closed Newton-Cotes (arbitrary degree)
- open Newton-Cotes (arbitrary degree)
You can use orthopy to generate Gauss formulas for your own weight functions.
Example:
val = quadpy.line_segment.integrate(
lambda x: numpy.exp(x),
[0.0, 1.0],
quadpy.line_segment.GaussPatterson(5)
)
- Krylov (1959, arbitrary degree)
Example:
val = quadpy.circle.integrate(
lambda x: numpy.exp(x[0]),
[0.0, 0.0], 1.0,
quadpy.circle.Krylov(7)
)
Apart from the classical centroid, vertex, and seven-point schemes we have
- Hammer-Marlowe-Stroud (1956, 5 schemes up to degree 5),
- Hammer-Stroud (1958, 2 schemes up to degree 3)
- open and closed Newton-Cotes schemes (1970, after Silvester, arbitrary degree),
- via Stroud (1971):
- Albrecht-Collatz (1958, degree 3)
- conical product scheme (degree 7)
- Strang/Cowper (1973, 10 schemes up to degree 7),
- Lyness-Jespersen (1975, 21 schemes up to degree 11),
- Lether (1976, degree 2n-2, arbitrary n, not symmetric)
- Hillion (1977, 8 schemes up to degree 3),
- Grundmann-Möller (1978, arbitrary degree),
- Laursen-Gellert (1978, 17 schemes up to degree 10),
- CUBTRI (Laurie, 1982, degree 8),
- TRIEX (de Doncker-Robinson, 1984, degrees 9 and 11),
- Dunavant (1985, 20 schemes up to degree 20),
- Cools-Haegemans (1987, degrees 8 and 11),
- Gatermann (1988, degree 7)
- Berntsen-Espelid (1990, 4 schemes of degree 13, the first one being DCUTRI),
- Liu-Vinokur (1998, 13 schemes up to degree 5),
- Walkington (2000, 5 schemes up to degree 5),
- Wandzura-Xiao (2003, 6 schemes up to degree 30),
- Taylor-Wingate-Bos (2005, 5 schemes up to degree 14),
- Zhang-Cui-Liu (2009, 3 schemes up to degree 20),
- Xiao-Gimbutas (2010, 50 schemes up to degree 50),
- Vioreanu-Rokhlin (2014, 20 schemes up to degree 62),
- Williams-Shunn-Jameson (2014, 8 schemes up to degree 12),
- Witherden-Vincent (2015, 19 schemes up to degree 20),
- Papanicolopulos (2016, 27 schemes up to degree 25).
Example:
val = quadpy.triangle.integrate(
lambda x: numpy.exp(x[0]),
[[0.0, 0.0], [1.0, 0.0], [0.5, 0.7]],
quadpy.triangle.XiaoGimbutas(5)
)
- Peirce (1957, arbitrary degree)
- via Stroud:
- Radon (1948, degree 5)
- Peirce (1956, 3 schemes up to degree 11)
- Albrecht-Collatz (1958, degree 3)
- Hammer-Stroud (1958, 8 schemes up to degree 15)
- Albrecht (1960, 8 schemes up to degree 17)
- Mysovskih (1964, 3 schemes up to degree 15)
- Rabinowitz-Richter (1969, 6 schemes up to degree 15)
- Lether (1971, arbitrary degree)
- Cools-Haegemans (1985, 3 schemes up to degree 9)
- Wissmann-Becker (1986, 3 schemes up to degree 8)
- Cools-Kim (2000, 3 schemes up to degree 21)
Example:
val = quadpy.disk.integrate(
lambda x: numpy.exp(x[0]),
[0.0, 0.0], 1.0,
quadpy.disk.Lether(6)
)
- Hammer-Stroud (1958, 3 schemes up to degree 7)
- via Stroud (1971, 15 schemes up to degree 15):
- Maxwell (1890, degree 7)
- Burnside (1908, degree 5)
- Irwin (1923, 3 schemes up to degree 5)
- Tyler (1953, 3 schemes up to degree 7)
- Albrecht-Collatz (1958, 4 schemes up to degree 5)
- Miller (1960, degree 1)
- Meister (1966, degree 7)
- Phillips (1967, degree 7)
- Rabinowitz-Richter (1969, 6 schemes up to degree 15)
- Cools-Haegemans (1985, 3 schemes up to degree 13)
- Dunavant (1985, 11 schemes up to degree 19)
- Morrow-Patterson (1985, 2 schemes up to degree 20, single precision)
- Wissmann-Becker (1986, 6 schemes up to degree 8)
- Cools-Haegemans (1988, 2 schemes up to degree 13)
- products of line segment schemes
- all formulas from the n-cube
Example:
val = quadpy.quadrilateral.integrate(
lambda x: numpy.exp(x[0]),
[[[0.0, 0.0], [1.0, 0.0]], [[0.0, 1.0], [1.0, 1.0]]],
quadpy.quadrilateral.Stroud('C2 7-2')
)
The points are specified in an array of shape (2, 2, ...) such that arr[0][0]
is the lower left corner, arr[1][1]
the upper right. If your quadrilateral
has its sides aligned with the coordinate axes, you can use the convenience
function
quadpy.quadrilateral.rectangle_points([x0, x1], [y0, y1])
to generate the array.
- via Stroud (1971):
- Stroud-Secrest (1963, 2 schemes up to degree 7)
- Rabinowitz-Richter (1969, 4 schemes up to degree 15)
- a scheme of degree 4
Example:
val = quadpy.e2r.integrate(
lambda x: numpy.exp(x[0]),
quadpy.e2r.RabinowitzRichter(5)
)
- via Stroud (1971):
- Stroud-Secrest (1963, 2 schemes up to degree 7)
- Rabinowitz-Richter (1969, 5 schemes up to degree 15)
- 3 schemes up to degree 7
Example:
val = quadpy.e2r2.integrate(
lambda x: numpy.exp(x[0]),
quadpy.e2r2.RabinowitzRichter(3)
)
- via Stroud (1971):
- Albrecht-Collatz (1958, 5 schemes up to degree 7)
- McLaren (1963, 10 schemes up to degree 14)
- Lebedev (1976, 32 schemes up to degree 131)
- Heo-Xu (2001, 27 schemes up to degree 39, single-precision)
Example:
val = quadpy.sphere.integrate(
lambda x: numpy.exp(x[0]),
[0.0, 0.0, 0.0], 1.0,
quadpy.sphere.Lebedev(19)
)
Integration on the sphere can also be done for function defined in spherical coordinates:
val = quadpy.sphere.integrate_spherical(
lambda phi_theta: numpy.sin(phi_theta[0])**2 * numpy.sin(phi_theta[1]),
radius=1.0,
rule=quadpy.sphere.Lebedev(19)
)
Note that phi_theta[0]
is the azimuthal, phi_theta[1]
the polar angle here.
- 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:
val = quadpy.ball.integrate(
lambda x: numpy.exp(x[0]),
[0.0, 0.0, 0.0], 1.0,
quadpy.ball.HammerStroud('14-3a')
)
- Hammer-Marlowe-Stroud (1956, 3 schemes up to degree 3)
- Hammer-Stroud (1958, 2 schemes up to degree 3)
- open and closed Newton-Cotes (1970, after Silvester) (arbitrary degree)
- Stroud (1971, degree 7)
- Grundmann-Möller (1978, arbitrary degree),
- Yu (1984, 5 schemes up to degree 6)
- Keast (1986, 11 schemes up to degree 8)
- Beckers-Haegemans (1990, degrees 8 and 9)
- Gatermann (1992, degree 5)
- Liu-Vinokur (1998, 14 schemes up to degree 5)
- Walkington (2000, 6 schemes up to degree 7)
- Zienkiewicz (2005, 2 schemes up to degree 3)
- Zhang-Cui-Liu (2009, 2 schemes up to degree 14)
- Xiao-Gimbutas (2010, 15 schemes up to degree 15)
- Shunn-Ham (2012, 6 schemes up to degree 7)
- Vioreanu-Rokhlin (2014, 10 schemes up to degree 13)
- Williams-Shunn-Jameson (2014, 1 scheme with degree 9)
- Witherden-Vincent (2015, 9 schemes up to degree 10)
Example:
val = quadpy.tetrahedron.integrate(
lambda x: numpy.exp(x[0]),
[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.5, 0.7, 0.0], [0.3, 0.9, 1.0]],
quadpy.tetrahedron.Keast(10)
)
- Product schemes derived from line segment schemes
- via: Stroud (1971):
- Sadowsky (1940, degree 5)
- Tyler (1953, 2 schemes up to degree 5)
- Hammer-Wymore (1957, degree 7)
- Albrecht-Collatz (1958, degree 3)
- Hammer-Stroud (1958, 6 schemes up to degree 7)
- Mustard-Lyness-Blatt (1963, 6 schemes up to degree 5)
- Stroud (1967, degree 5)
- Sarma-Stroud (1969, degree 7)
- all formulas from the n-cube
Example:
val = quadpy.hexahedron.integrate(
lambda x: numpy.exp(x[0]),
quadpy.hexahedron.cube_points([0.0, 1.0], [-0.3, 0.4], [1.0, 2.1]),
quadpy.hexahedron.Product(quadpy.line_segment.NewtonCotesClosed(3))
)
- Felippa (9 schemes up to degree 5)
Example:
val = quadpy.pyramid.integrate(
lambda x: numpy.exp(x[0]),
[
[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],
],
quadpy.pyramid.Felippa(5)
)
- Felippa (6 schemes up to degree 6)
Example:
val = quadpy.wedge.integrate(
lambda x: numpy.exp(x[0]),
[
[[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]],
],
quadpy.wedge.Felippa(3)
)
- via Stroud (1971):
- Stroud-Secrest (1963, 5 schemes up to degree 7)
Example:
val = quadpy.e2r.integrate(
lambda x: numpy.exp(x[0]),
quadpy.e2r.StroudSecrest('IX')
)
- via Stroud (1971):
- Stroud-Secrest (1963, 7 schemes up to degree 7)
- scheme of degree 14
Example:
val = quadpy.e2r2.integrate(
lambda x: numpy.exp(x[0]),
quadpy.e2r2.RabinowitzRichter(3)
)
- via Stroud:
- Lauffer (1955, 5 schemes up to degree 5)
- Hammer-Stroud (1956, 3 schemes up to degree 3)
- Stroud (1966, 7 schemes of degree 3)
- Stroud (1969, degree 5)
- Grundmann-Möller (1978, arbitrary degree)
- Walkington (2000, 5 schemes up to degree 7)
Example:
dim = 4
val = quadpy.simplex.integrate(
lambda x: numpy.exp(x[0]),
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],
]),
quadpy.simplex.GrundmannMoeller(dim, 3)
)
Example:
dim = 4
quadpy.nsphere.integrate(
lambda x: numpy.exp(x[0]),
numpy.zeros(dim), 1.0,
quadpy.nsphere.Dobrodeev1978(dim)
)
Example:
dim = 4
quadpy.nball.integrate(
lambda x: numpy.exp(x[0]),
numpy.zeros(dim), 1.0,
quadpy.nball.Dobrodeev1970(dim)
)
- Dobrodeev (1970, n >= 5, degree 7)
- via Stroud (1971):
- Ewing (1941, degree 3)
- Tyler (1953, degree 3)
- Stroud (1957, 2 schemes up to degree 3)
- Hammer-Stroud (1958, degree 5)
- Mustard-Lyness-Blatt (1963, degree 5)
- Thacher (1964, degree 2)
- Stroud (1966, 4 schemes of degree 5)
- Phillips (1967, degree 7, single precision)
- Stroud (1968, degree 5)
- Dobrodeev (1978, n >= 2, degree 5)
Example:
dim = 4
quadpy.ncube.integrate(
lambda x: numpy.exp(x[0]),
quadpy.ncube.ncube_points(
[0.0, 1.0], [0.1, 0.9], [-1.0, 1.0], [-1.0, -0.5]
),
quadpy.ncube.Stroud(dim, 'Cn 3-3')
)
- via Stroud (1971):
- Stroud-Secrest (1963, 4 schemes up to degree 5)
- 2 schemes up to degree 5
Example:
dim = 4
val = quadpy.enr.integrate(
lambda x: numpy.exp(x[0]),
quadpy.enr.Stroud(dim, '5-4')
)
- 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 = quadpy.enr2.integrate(
lambda x: numpy.exp(x[0]),
quadpy.enr2.Stroud(dim, '5-2')
)
quadpy is available from the Python Package Index, so with
pip install -U quadpy
you can install/upgrade.
To run the tests, just check out this repository and type
MPLBACKEND=Agg pytest
To create a new release
-
bump the
__version__
number, -
publish to PyPi and GitHub:
$ make publish
quadpy is published under the MIT license.