quadpy

Quadpy is fully vectorized. You can integrate a function over multiple domains simultaneously by passing an array of domains to the scheme.integrate() method. Similarly, you can integrate functions that return vectorized outputs (e.g., a list or array of values).

# shape (3, 5, 2), i.e., (corners, num_triangles, xy_coords)
triangles = np.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,
)

# Example of a vectorized function
def f(x):
    return [np.sin(x[0]), np.sin(x[1])]

# Integration call would use the 'triangles' array
# val = scheme.integrate(f, triangles)

For basic numerical integration over a 1D interval, use the quadpy.quad function. This is similar to scipy.integrate.quad but supports complex-, vector-, and matrix-valued integrands, as well as arbitrary-dimensional intervals.

import numpy as np
import quadpy


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

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

To integrate over specific geometric domains (e.g., a triangle), you must first obtain a quadrature scheme and then call its .integrate() method. Most domains provide a get_good_scheme(degree) method to retrieve a scheme optimized for a specific algebraic degree. Alternatively, you can select a specific scheme from the domain's scheme dictionary (e.g., quadpy.t2.schemes).

import numpy as np
import quadpy


def f(x):
    return np.sin(x[0]) * np.sin(x[1])


triangle = np.array([[0.0, 0.0], [1.0, 0.0], [0.7, 0.5]])

# get a "good" scheme of degree 10
scheme = quadpy.t2.get_good_scheme(10)
val = scheme.integrate(f, triangle)

Install the quadpy package from PyPI using pip.

pip install quadpy

Use the quadpy.c1 module to perform integration on 1D line segments. Supported schemes include Chebyshev-Gauss, Clenshaw-Curtis, Gauss-Legendre, Gauss-Lobatto, and Gauss-Patterson (up to degree 767).

import numpy as np
import quadpy

scheme = quadpy.c1.gauss_patterson(5)
scheme.show()
val = scheme.integrate(lambda x: np.exp(x), [0.0, 1.0])

Use the quadpy.p3 module for integration over a 3D pyramid. The felippa_5() scheme is available. Integration limits are a list of the five vertices.

import numpy as np
import quadpy

scheme = quadpy.p3.felippa_5()
val = scheme.integrate(
    lambda x: np.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],
    ],
)

Use the quadpy.t3 module for integration over a 3D tetrahedron. Use get_good_scheme(degree) to select a scheme. Integration limits are a list of the four vertices: [[x1, y1, z1], [x2, y2, z2], [x3, y3, z3], [x4, y4, z4]].

import numpy as np
import quadpy

scheme = quadpy.t3.get_good_scheme(5)
val = scheme.integrate(
    lambda x: np.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]],
)

Use the quadpy.cn module for integration over an n-dimensional cube. The stroud_cn_3_3(dim) scheme is available. Integration limits are provided via quadpy.cn.ncube_points(...) which takes ranges for each dimension.

import numpy as np
import quadpy

dim = 4
scheme = quadpy.cn.stroud_cn_3_3(dim)
val = scheme.integrate(
    lambda x: np.exp(x[0]),
    quadpy.cn.ncube_points([0.0, 1.0], [0.1, 0.9], [-1.0, 1.0], [-1.0, -0.5]),
)

Once you have a scheme object, you can access several attributes and methods to understand its composition and behavior:

• scheme.points: The quadrature points.
• scheme.weights: The quadrature weights.
• scheme.degree: The algebraic degree of the scheme.
• scheme.source: The bibliographic source of the scheme.
• scheme.test_tolerance: The tolerance used for testing.
• scheme.show(): Displays information about the scheme.
• scheme.integrate(f, domain): Performs the integration.
• scheme.points_symbolic: Symbolic representation of points.
• scheme.weights_symbolic: Symbolic representation of weights.

Use the quadpy.t2 module for integration over 2D triangles. Use get_good_scheme(degree) to select a scheme. The integration limits should be a list of the three vertices: [[x1, y1], [x2, y2], [x3, y3]].

import numpy as np
import quadpy

scheme = quadpy.t2.get_good_scheme(12)
scheme.show()
val = scheme.integrate(lambda x: np.exp(x[0]), [[0.0, 0.0], [1.0, 0.0], [0.5, 0.7]])

Use the quadpy.c2 module for integration over 2D quadrilaterals. Limits are specified as an array of shape (2, 2, ...) where arr[0][0] is the lower-left corner and arr[1][1] is the upper-right corner. For axis-aligned rectangles, use quadpy.c2.rectangle_points([x0, x1], [y0, y1]) to generate the points array.

import numpy as np
import quadpy

scheme = quadpy.c2.get_good_scheme(7)
val = scheme.integrate(
    lambda x: np.exp(x[0]),
    [[[0.0, 0.0], [1.0, 0.0]], [[0.0, 1.0], [1.0, 1.0]]],
)

Use the quadpy.e1r2 module for integration on a 1D space with the weight function $e^{-r^2}$. Supported schemes include Gauss-Hermite and Genz-Keister.

import quadpy

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