Skip to content
/ quadpy Public
forked from LJPapenfort/quadpy

Numerical integration (quadrature, cubature) in Python

License

MIT license
0 stars 82 forks Branches Tags Activity
Star
Notifications You must be signed in to change notification settings

ermost/quadpy

BranchesTags
 
 

Repository files navigation

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.

About

Numerical integration (quadrature, cubature) in Python

Resources

Readme

License

MIT license

Code of conduct

Code of conduct
Activity

Stars

0 stars

Watchers

0 watching

Forks

0 forks
Report repository

Releases

55 tags

Packages

No packages published

Languages

  • Python 99.9%
  • Makefile 0.1%