Skip to content
Numerical integration (quadrature, cubature) in Python
Branch: master
Clone or download
Latest commit b24b944 Jul 2, 2019
Permalink
Type Name Latest commit message Commit time
Failed to load latest commit information.
.circleci circleci fix Jun 25, 2019
quadpy version bump Jul 2, 2019
test rename secrest to roman numerals Jul 1, 2019
tools remove -*- coding: utf-8 -*- Jul 1, 2019
.bandit add bandit config Feb 14, 2018
.flake8 flake9: ignore too-complex (C901) Jun 10, 2019
.gitignore more gitignore Jun 4, 2018
.travis.yml travis: install pylint Jul 13, 2017
LICENSE fix e3r2 Jun 20, 2019
Makefile small fixes for scheme_from_rc Jun 27, 2019
README.md add on logo Jul 2, 2019
codecov.yml codecov config Jun 4, 2019
setup.py remove -*- coding: utf-8 -*- Jul 1, 2019
test_requirements.txt install accupy for tests Jun 6, 2018

README.md

quadpy

Your one-stop shop for numerical integration in Python.

CircleCI codecov Code style: black awesome PyPi Version DOI GitHub stars PyPi downloads

Maxwell's quad scheme from 1890

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
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.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[0]), numpy.sin(x[1])]

More examples under test/examples_test.py.

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

Advanced topics:

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
import quadpy

scheme = quadpy.line_segment.gauss_patterson(5)
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:

import quadpy

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:

import quadpy

scheme = quadpy.e1r2.gauss_hermite(5)
scheme.show()
val = scheme.integrate(lambda x: x**2)

Circle

  • Krylov (1959, arbitrary degree)

Example:

import quadpy

scheme = quadpy.circle.krylov(7)
scheme.show()
val = scheme.integrate(lambda x: numpy.exp(x[0]), [0.0, 0.0], 1.0)

Triangle

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

Example:

import quadpy

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

Disk

Example:

import numpy
import quadpy

scheme = quadpy.disk.lether(6)
scheme.show()
val = scheme.integrate(lambda x: numpy.exp(x[0]), [0.0, 0.0], 1.0)

Quadrilateral

Example:

import numpy
import quadpy

scheme = quadpy.quadrilateral.stroud_c2_7_2()
val = scheme.integrate(
    lambda x: numpy.exp(x[0]),
    [[[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[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)

Example:

import quadpy

scheme = quadpy.e2r.rabinowitz_richter_5()
scheme.show()
val = scheme.integrate(lambda x: x[0]**2)

2D space with weight function exp(-r2)

Example:

import quadpy

scheme = quadpy.e2r2.rabinowitz_richter_3()
scheme.show()
val = scheme.integrate(lambda x: x[0]**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
import quadpy

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

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

import numpy
import quadpy

scheme = quadpy.sphere.lebedev_019()
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
import quadpy

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

Tetrahedron

Example:

import numpy
import quadpy

scheme = quadpy.tetrahedron.keast_10()
scheme.show()
val = scheme.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]],
    )

Hexahedron

Example:

import numpy
import quadpy

scheme = quadpy.hexahedron.product(quadpy.line_segment.newton_cotes_closed(3))
scheme.show()
val = scheme.integrate(
    lambda x: numpy.exp(x[0]),
    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
import quadpy

scheme = quadpy.pyramid.felippa_5()

val = scheme.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],
    ]
    )

Wedge

Example:

import numpy
import quadpy

scheme = quadpy.wedge.felippa_3
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]],
    ]
    )

3D space with weight function exp(-r)

Example:

import quadpy

scheme = quadpy.e3r.stroud_secrest_09()
scheme.show()
val = scheme.integrate(lambda x: x[0]**2)

3D space with weight function exp(-r2)

Example:

import quadpy

scheme = quadpy.e3r2.stroud_secrest_10a
scheme.show()
val = scheme.integrate(lambda x: x[0]**2)

n-Simplex

Example:

import numpy
import quadpy

scheme = quadpy.nsimplex.GrundmannMoeller(dim, 3)
dim = 4
val = scheme.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],
        ])
    )

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
import quadpy

scheme = quadpy.nsphere.dobrodeev_1978(dim)
dim = 4
val = scheme.integrate(lambda x: numpy.exp(x[0]), 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
import quadpy

scheme = quadpy.nball.dobrodeev_1970(dim)
dim = 4
val = scheme.integrate(lambda x: numpy.exp(x[0]), numpy.zeros(dim), 1.0)

n-Cube

Example:

import numpy
import quadpy

dim = 4
scheme = quadpy.ncube.stroud_cn_3_3(dim)
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]
        )
    )

nD space with weight function exp(-r)

Example:

import quadpy

dim = 4
scheme = quadpy.enr.stroud_5_4(dim)
val = scheme.integrate(lambda x: x[0]**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:

import quadpy

dim = 4
scheme = quadpy.enr2.stroud_5_2(dim)
val = scheme.integrate(lambda x: x[0]**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

License

quadpy is published under the MIT license.

You can’t perform that action at this time.