Skip to main content

Numerical integration, quadrature for various domains

Project description

quadpy

Your one-stop shop for numerical integration in Python.

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

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.

Project details


Download files

Download the file for your platform. If you're not sure which to choose, learn more about installing packages.

Files for quadpy, version 0.13.2
Filename, size File type Python version Upload date Hashes
Filename, size quadpy-0.13.2-py3-none-any.whl (761.7 kB) File type Wheel Python version py3 Upload date Hashes View hashes
Filename, size quadpy-0.13.2.tar.gz (628.8 kB) File type Source Python version None Upload date Hashes View hashes

Supported by

Elastic Elastic Search Pingdom Pingdom Monitoring Google Google BigQuery Sentry Sentry Error logging AWS AWS Cloud computing DataDog DataDog Monitoring Fastly Fastly CDN SignalFx SignalFx Supporter DigiCert DigiCert EV certificate StatusPage StatusPage Status page