Numerical integration, quadrature for various domains
Project description
quadpy
Your onestop shop for numerical integration in Python.
Hundreds of numerical integration schemes for line segments, circles, disks, triangles, quadrilaterals, spheres, balls, tetrahedra, hexahedra, wedges, pyramids, nspheres, nballs, ncubes, nsimplices, and the 1D/2D/3D/nD spaces with weight functions exp(r) and exp(r^{2}).
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.,
# 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[0]), numpy.sin(x[1])]
More examples under test/examples_test.py.
Adaptive quadrature
quadpy can do adaptive quadrature for certain domains. Again, everything is fully vectorized, so you can provide multiple intervals and vectorvalued functions.
Line segments
val, error_estimate = quadpy.line_segment.integrate_adaptive( lambda x: x * sin(5 * x), [0.0, pi], 1.0e10 )
tanhsinh quadrature
The more modern tanhsinh quadrature is different from all other methods in quadpy in that it doesn't exactly integrate any function exactly, not even polynomials of low degree. Its tremendous usefulness rather comes from the fact that a wide variety of function, even seemingly difficult ones with (integrable) singularities at the end points, can be integrated with arbitrary precision.
from mpmath import mp import sympy mp.dps = 50 val, error_estimate = quadpy.line_segment.tanh_sinh( lambda x: mp.exp(x) * sympy.cos(x), 0, mp.pi/2, 1.0e50 # ! )
Note the usage of mpmath
here for arbirtrary precision arithmetics.
If the function has a singularity at a boundary, it needs to be shifted such
that the singularity is at 0. If there are singularities at both ends, the
function can be shifted both ways and be handed off to tanh_sinh_lr
:
tanh_sinh_lr(f_left, f_right, interval_length, tol)
Triangles
val, error_estimate = quadpy.triangle.integrate_adaptive( lambda x: x[0] * sin(5 * x[1]), [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]], 1.0e10 )
ProTip: You can provide many triangles that together form a domain to get an approximation of the integral over the domain.
Schemes
Line segment
 ChebyshevGauss (both variants, arbitrary degree)
 ClenshawCurtis (after Waldvogel, arbitrary degree)
 Fejértype1 (after Waldvogel, arbitrary degree)
 Fejértype2 (after Waldvogel, arbitrary degree)
 GaussJacobi
 GaussLegendre (via NumPy, arbitrary degree)
 GaussLobatto (arbitrary degree)
 GaussKronrod (after Laurie, arbitrary degree)
 GaussPatterson (7 schemes up to degree 191)
 GaussRadau (arbitrary degree)
 closed NewtonCotes (arbitrary degree)
 open NewtonCotes (arbitrary degree)
 tanhsinh quadrature (see above)
See below for how 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) )
1D halfspace with weight function exp(r)
 Generalized GaussLaguerre
Example:
val = quadpy.e1r.integrate( lambda x: x**2, quadpy.e1r.GaussLaguerre(5, alpha=0) )
1D space with weight function exp(r^{2})
 GaussHermite (via NumPy, arbitrary degree)
Example:
val = quadpy.e1r2.integrate( lambda x: x**2, quadpy.e1r2.GaussHermite(5) )
Circle
 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) )
Triangle
Apart from the classical centroid, vertex, and sevenpoint schemes we have
 HammerMarloweStroud (1956, 5 schemes up to degree 5),
 HammerStroud (1958, 2 schemes up to degree 3)
 open and closed NewtonCotes schemes (1970, after Silvester, arbitrary degree),
 via Stroud (1971):
 AlbrechtCollatz (1958, degree 3)
 conical product scheme (degree 7)
 Strang/Cowper (1973, 10 schemes up to degree 7),
 LynessJespersen (1975, 21 schemes up to degree 11),
 Lether (1976, degree 2n2, arbitrary n, not symmetric)
 Hillion (1977, 10 schemes up to degree 3),
 GrundmannMöller (1978, arbitrary degree),
 LaursenGellert (1978, 17 schemes up to degree 10),
 CUBTRI (Laurie, 1982, degree 8),
 TRIEX (de DonckerRobinson, 1984, degrees 9 and 11),
 Dunavant (1985, 20 schemes up to degree 20),
 CoolsHaegemans (1987, degrees 8 and 11),
 Gatermann (1988, degree 7)
 BerntsenEspelid (1990, 4 schemes of degree 13, the first one being DCUTRI),
 LiuVinokur (1998, 13 schemes up to degree 5),
 Walkington (2000, 5 schemes up to degree 5),
 WandzuraXiao (2003, 6 schemes up to degree 30),
 TaylorWingateBos (2005, 5 schemes up to degree 14),
 ZhangCuiLiu (2009, 3 schemes up to degree 20),
 XiaoGimbutas (2010, 50 schemes up to degree 50),
 VioreanuRokhlin (2014, 20 schemes up to degree 62),
 WilliamsShunnJameson (2014, 8 schemes up to degree 12),
 WitherdenVincent (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) )
Disk
 Peirce (1957, arbitrary degree)
 via Stroud:
 Radon (1948, degree 5)
 Peirce (1956, 3 schemes up to degree 11)
 AlbrechtCollatz (1958, degree 3)
 HammerStroud (1958, 8 schemes up to degree 15)
 Albrecht (1960, 8 schemes up to degree 17)
 Mysovskih (1964, 3 schemes up to degree 15)
 RabinowitzRichter (1969, 6 schemes up to degree 15)
 Lether (1971, arbitrary degree)
 CoolsHaegemans (1985, 3 schemes up to degree 9)
 WissmannBecker (1986, 3 schemes up to degree 8)
 CoolsKim (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) )
Quadrilateral
 HammerStroud (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)
 AlbrechtCollatz (1958, 4 schemes up to degree 5)
 Miller (1960, degree 1)
 Meister (1966, degree 7)
 Phillips (1967, degree 7)
 RabinowitzRichter (1969, 6 schemes up to degree 15)
 CoolsHaegemans (1985, 3 schemes up to degree 13)
 Dunavant (1985, 11 schemes up to degree 19)
 MorrowPatterson (1985, 2 schemes up to degree 20, single precision)
 WissmannBecker (1986, 6 schemes up to degree 8)
 CoolsHaegemans (1988, 2 schemes up to degree 13)
 products of line segment schemes
 all formulas from the ncube
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 72') )
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.
2D space with weight function exp(r)
 via Stroud (1971):
 StroudSecrest (1963, 2 schemes up to degree 7)
 RabinowitzRichter (1969, 4 schemes up to degree 15)
 a scheme of degree 4
Example:
val = quadpy.e2r.integrate( lambda x: x[0]**2, quadpy.e2r.RabinowitzRichter(5) )
2D space with weight function exp(r^{2})
 via Stroud (1971):
 StroudSecrest (1963, 2 schemes up to degree 7)
 RabinowitzRichter (1969, 5 schemes up to degree 15)
 3 schemes up to degree 7
Example:
val = quadpy.e2r2.integrate( lambda x: x[0]**2, quadpy.e2r2.RabinowitzRichter(3) )
Sphere
 via Stroud (1971):
 AlbrechtCollatz (1958, 5 schemes up to degree 7)
 McLaren (1963, 10 schemes up to degree 14)
 Lebedev (1976, 34 schemes up to degree 131)
 HeoXu (2001, 27 schemes up to degree 39, singleprecision)
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 azimuthal, polar: numpy.sin(azimuthal)**2 * numpy.sin(polar), rule=quadpy.sphere.Lebedev("19") )
Ball
 HammerStroud (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('143a') )
Tetrahedron
 HammerMarloweStroud (1956, 3 schemes up to degree 3)
 HammerStroud (1958, 2 schemes up to degree 3)
 open and closed NewtonCotes (1970, after Silvester) (arbitrary degree)
 Stroud (1971, degree 7)
 GrundmannMöller (1978, arbitrary degree),
 Yu (1984, 5 schemes up to degree 6)
 Keast (1986, 11 schemes up to degree 8)
 BeckersHaegemans (1990, degrees 8 and 9)
 Gatermann (1992, degree 5)
 LiuVinokur (1998, 14 schemes up to degree 5)
 Walkington (2000, 6 schemes up to degree 7)
 Zienkiewicz (2005, 2 schemes up to degree 3)
 ZhangCuiLiu (2009, 2 schemes up to degree 14)
 XiaoGimbutas (2010, 15 schemes up to degree 15)
 ShunnHam (2012, 6 schemes up to degree 7)
 VioreanuRokhlin (2014, 10 schemes up to degree 13)
 WilliamsShunnJameson (2014, 1 scheme with degree 9)
 WitherdenVincent (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) )
Hexahedron
 Product schemes derived from line segment schemes
 via: Stroud (1971):
 Sadowsky (1940, degree 5)
 Tyler (1953, 2 schemes up to degree 5)
 HammerWymore (1957, degree 7)
 AlbrechtCollatz (1958, degree 3)
 HammerStroud (1958, 6 schemes up to degree 7)
 MustardLynessBlatt (1963, 6 schemes up to degree 5)
 Stroud (1967, degree 5)
 SarmaStroud (1969, degree 7)
 all formulas from the ncube
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)) )
Pyramid
 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) )
Wedge
 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) )
3D space with weight function exp(r)
 via Stroud (1971):
 StroudSecrest (1963, 5 schemes up to degree 7)
Example:
val = quadpy.e3r.integrate( lambda x: x[0]**2, quadpy.e3r.StroudSecrest('IX') )
3D space with weight function exp(r^{2})
 via Stroud (1971):
 StroudSecrest (1963, 7 schemes up to degree 7)
 scheme of degree 14
Example:
val = quadpy.e3r2.integrate( lambda x: x[0]**2, quadpy.e3r2.StroudSecrest("Xa") )
nSimplex
 via Stroud:
 Lauffer (1955, 5 schemes up to degree 5)
 HammerStroud (1956, 3 schemes up to degree 3)
 Stroud (1964, degree 3)
 Stroud (1966, 7 schemes of degree 3)
 Stroud (1969, degree 5)
 GrundmannMöller (1978, arbitrary degree)
 Walkington (2000, 5 schemes up to degree 7)
Example:
dim = 4 val = quadpy.nsimplex.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.nsimplex.GrundmannMoeller(dim, 3) )
nSphere
Example:
dim = 4 quadpy.nsphere.integrate( lambda x: numpy.exp(x[0]), numpy.zeros(dim), 1.0, quadpy.nsphere.Dobrodeev1978(dim) )
nBall
Example:
dim = 4 quadpy.nball.integrate( lambda x: numpy.exp(x[0]), numpy.zeros(dim), 1.0, quadpy.nball.Dobrodeev1970(dim) )
nCube
 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)
 HammerStroud (1958, degree 5)
 MustardLynessBlatt (1963, degree 5)
 Thacher (1964, degree 2)
 Stroud (1966, 4 schemes of degree 5)
 Phillips (1967, degree 7)
 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 33') )
nD space with weight function exp(r)
 via Stroud (1971):
 StroudSecrest (1963, 4 schemes up to degree 5)
 2 schemes up to degree 5
Example:
dim = 4 val = quadpy.enr.integrate( lambda x: x[0]**2, quadpy.enr.Stroud(dim, '54') )
nD space with weight function exp(r^{2})
 via Stroud (1971):
 StroudSecrest (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: x[0]**2, quadpy.enr2.Stroud(dim, '52') )
Extras
Classical schemes
With quadpy, it's easy to regenerate classical Gauss quadrature schemes are listed in, e.g., Stroud & Secrest.
Some examples:
scheme = quadpy.line_segment.GaussLegendre(96, mode='mpmath', decimal_places=30) scheme = quadpy.e1r2.GaussHermite(14, mode='mpmath', decimal_places=20) scheme = quadpy.e1r.GaussLaguerre(13, mode='mpmath', decimal_places=50)
Generating your own Gauss quadrature in three simple steps
You have a measure (or, more colloquially speaking, a domain and a nonnegative weight function) and would like to generate the matching Gauss quadrature? Great, here's how to do it.
As an example, let's try and generate the Gauss quadrature with 10 points for
the weight function x^2
on the interval [1, +1]
.
TLDR:
moments = quadpy.tools.integrate( lambda x: [x**(2+k) for k in range(20)], 1, +1 ) alpha, beta = quadpy.tools.chebyshev(moments) points, weights = quadpy.tools.scheme_from_rc(alpha, beta, decimal_places=20)
Some explanations:

You need to compute the first
2*n
moments of your measureintegral(w(x) p_k(x) dx)
with a particular set of polynomials
p_k
. A common choice are the monomialsx^k
. You can do that by hand or usemoments = quadpy.tools.integrate(lambda x: [x**(2+k) for k in range(20)], 1, +1)
[2/3, 0, 2/5, 0, 2/7, 0, 2/9, 0, 2/11, 0, 2/13, 0, 2/15, 0, 2/17, 0, 2/19, 0, 2/21, 0]
Note that the moments have all been computed symbolically here.
If you have the moments in floating point (for example because you need to compute the scheme fast), it makes sense to think about the numerical implications here. That's because the map to the recurrence coefficients (step 2) can be very illconditioned, meaning that small roundoff errors can lead to an unusable scheme. For further computation, it's numerically beneficial if the moments are either 0 or in the same order of magnitude. The above numbers are alright, but if you want to max it out, you could try Legendre polynomials from orthopy for
p_k
:import orthopy def leg_polys(x): return orthopy.line_segment.tree_legendre(x, 20, "monic", symbolic=True) moments = quadpy.tools.integrate( lambda x: [x**2 * leg_poly for leg_poly in leg_polys(x)], 1, +1 )
[2/3, 0, 8/45, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
Better!

From the moments, we generate the recurrence coefficients of our custom orthogonal polynomials. There are a few choices to accomplish this:
golub_welsch
: uses Cholesky at its core; can be numerically unstablestieltjes
: moments not even needed here, but can also be numerically unstablechebyshev
: can be used if you chose monomials in the first step; again, potentially numerically unstablechebyshev_modified
: to be used if you chose something other than monomials in the first step; stable if thepolynomial_class
was chosen wisely
Since we have computed modified moments in step one, let's use the latter method:
_, _, a, b = \ orthopy.line_segment.recurrence_coefficients.legendre(20, "monic", symbolic=True) alpha, beta = quadpy.tools.chebyshev_modified(moments, a, b)
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0] [2/3, 3/5, 4/35, 25/63, 16/99, 49/143, 12/65, 27/85, 64/323, 121/399]
(Note that, since everything is done symbolically in this example, we could have used Stieltjes's or Chebyshev's unmodified method; the results are the same.)

Lastly, we generate the Gauss points and weights from
alpha
andbeta
. Since symbolic computation can take very long even for small sizes, we convertalpha
andbeta
to numpy arrays first. (If you need more digits, look at mpmath arrays.)points, weights = quadpy.tools.scheme_from_rc( numpy.array([sympy.N(a) for a in alpha], dtype=float), numpy.array([sympy.N(b) for b in beta], dtype=float), mode='numpy' )
[0.97822866 0.8870626 0.73015201 0.51909613 0.26954316 0.26954316 0.51909613 0.73015201 0.8870626 0.97822866]
[0.05327099 0.09881669 0.0993154 0.06283658 0.01909367 0.01909367 0.06283658 0.0993154 0.09881669 0.05327099]
Congratulations! Your Gaussian quadrature rule.
Other tools

Transforming Gaussian points and weights back to recurrence coefficients:
alpha, beta = quadpy.tools.coefficients_from_gauss(points, weights)

The Gautschi test: As recommended by Gautschi, you can test your momentbased scheme with
err = quadpy.tools.check_coefficients(moments, alpha, beta)
Relevant publications
 A.H. Stroud and D. Secrest, Gaussian Quadrature Formulas, 1966, Prentice Hall, Series in Automatic Computation
 Gene H. Golub and John H. Welsch, Calculation of Gauss Quadrature Rules, Mathematics of Computation, Vol. 23, No. 106 (Apr., 1969), pp. 221230+s1s10
 W. Gautschi, On Generating Orthogonal Polynomials, SIAM J. Sci. and Stat. Comput., 3(3), 289–317
 W. Gautschi, How and how not to check Gaussian quadrature formulae, BIT Numerical Mathematics, June 1983, Volume 23, Issue 2, pp 209–216
 D. Boley and G.H. Golub, A survey of matrix inverse eigenvalue problems, Inverse Problems, 1987, Volume 3, Number 4
 W. Gautschi, Algorithm 726: ORTHPOL–a package of routines for generating orthogonal polynomials and Gausstype quadrature rules, ACM Transactions on Mathematical Software (TOMS), Volume 20, Issue 1, March 1994, Pages 2162
Installation
quadpy is available from the Python Package Index, so with
pip install U quadpy
you can install/upgrade.
Testing
To run the tests, just check out this repository and type
MPLBACKEND=Agg pytest
Distribution
To create a new release

bump the
__version__
number, 
publish to PyPi and GitHub:
$ make publish
License
quadpy is published under the MIT license.
Project details
Release history Release notifications
Download files
Download the file for your platform. If you're not sure which to choose, learn more about installing packages.
Filename, size & hash SHA256 hash help  File type  Python version  Upload date 

quadpy0.12.6py2.py3noneany.whl (462.3 kB) Copy SHA256 hash SHA256  Wheel  py2.py3  Oct 10, 2018 
quadpy0.12.6.tar.gz (358.3 kB) Copy SHA256 hash SHA256  Source  None  Oct 10, 2018 