Skip to content
/ quadpy Public
forked from Programing4Mathematics/quadpy

Numerical integration (quadrature, cubature) in Python

License

GPL-3.0 license
0 stars 84 forks Branches Tags Activity
Star
Notifications You must be signed in to change notification settings

BonneelP/quadpy

BranchesTags
 
 

Repository files navigation

quadpy

Your one-stop shop for numerical integration in Python.

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

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 interval, install quadpy from the Python Package Index with

pip install quadpy

and do

import numpy
import quadpy

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

val, err = quadpy.quad(f, 0.0, 6.0)

This is just like scipy with the addition that quadpy handles complex-, vector-, matrix-valued integrands, and lines in spaces of arbitrary dimension.

To integrate over a triangle with Strang's rule of degree 6, 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)

All schemes 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_3()
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

pip install quadpy

you can install.

Testing

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

MPLBACKEND=Agg pytest

License

This software is published under the GPLv3 license.

About

Numerical integration (quadrature, cubature) in Python

Resources

Readme

License

GPL-3.0 license

Code of conduct

Code of conduct
Activity

Stars

0 stars

Watchers

0 watching

Forks

0 forks
Report repository

Releases

No releases published

Packages

No packages published

Languages

  • Python 99.9%
  • Makefile 0.1%