Skip to content
/ quadpy Public
forked from sigma-py/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

macicco/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

Hundreds of 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).

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 vector-valued functions.

Line segments

val, error_estimate = quadpy.line_segment.integrate_adaptive(
        lambda x: x * sin(5 * x),
        [0.0, pi],
        1.0e-10
        )

tanh-sinh quadrature

The more modern tanh-sinh 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.0e-50  # !
        )

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.0e-10
        )

ProTip: You can provide many triangles that together form a domain to get an approximation of the integral over the domain.

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 (7 schemes up to degree 191)
  • 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:

val = quadpy.line_segment.integrate(
    lambda x: numpy.exp(x),
    [0.0, 1.0],
    quadpy.line_segment.GaussPatterson(5)
    )

1D half-space with weight function exp(-r)

  • Generalized Gauss-Laguerre

Example:

val = quadpy.e1r.integrate(
    lambda x: x**2,
    quadpy.e1r.GaussLaguerre(5, alpha=0)
    )

1D space with weight function exp(-r2)

  • Gauss-Hermite (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 seven-point schemes we have

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

Example:

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

Quadrilateral

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 7-2')
    )

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:

val = quadpy.e2r.integrate(
    lambda x: x[0]**2,
    quadpy.e2r.RabinowitzRichter(5)
    )

2D space with weight function exp(-r2)

Example:

val = quadpy.e2r2.integrate(
    lambda x: x[0]**2,
    quadpy.e2r2.RabinowitzRichter(3)
    )

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:

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

  • 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:

val = quadpy.ball.integrate(
    lambda x: numpy.exp(x[0]),
    [0.0, 0.0, 0.0], 1.0,
    quadpy.ball.HammerStroud('14-3a')
    )

Tetrahedron

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

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)

Example:

val = quadpy.e3r.integrate(
    lambda x: x[0]**2,
    quadpy.e3r.StroudSecrest('IX')
    )

3D space with weight function exp(-r2)

Example:

val = quadpy.e3r2.integrate(
    lambda x: x[0]**2,
    quadpy.e3r2.StroudSecrest("Xa")
    )

n-Simplex

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)
    )

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:

dim = 4
quadpy.nsphere.integrate(
    lambda x: numpy.exp(x[0]),
    numpy.zeros(dim), 1.0,
    quadpy.nsphere.Dobrodeev1978(dim)
    )

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:

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

n-Cube

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 3-3')
    )

nD space with weight function exp(-r)

Example:

dim = 4
val = quadpy.enr.integrate(
    lambda x: x[0]**2,
    quadpy.enr.Stroud(dim, '5-4')
    )

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:

dim = 4
val = quadpy.enr2.integrate(
    lambda x: x[0]**2,
    quadpy.enr2.Stroud(dim, '5-2')
    )

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:

  1. You need to compute the first 2*n moments of your measure

    integral(w(x) p_k(x) dx)

    with a particular set of polynomials p_k. A common choice are the monomials x^k. You can do that by hand or use

    moments = 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 ill-conditioned, meaning that small round-off 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!

  2. 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 unstable
    • stieltjes: moments not even needed here, but can also be numerically unstable
    • chebyshev: can be used if you chose monomials in the first step; again, potentially numerically unstable
    • chebyshev_modified: to be used if you chose something other than monomials in the first step; stable if the polynomial_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.)

  3. Lastly, we generate the Gauss points and weights from alpha and beta. Since symbolic computation can take very long even for small sizes, we convert alpha and beta 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 moment-based scheme with

    err = quadpy.tools.check_coefficients(moments, alpha, beta)

Relevant publications

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

  1. bump the __version__ number,

  2. publish to PyPi and GitHub:

    $ make publish

License

quadpy is published under the MIT license.

About

Numerical integration (quadrature, cubature) in Python

Resources

Readme

License

MIT license
Activity

Stars

0 stars

Watchers

2 watching

Forks

0 forks
Report repository

Releases

51 tags

Packages

No packages published

Languages

  • Python 99.9%
  • Makefile 0.1%