Quadrature

Base classes

class modepy.Quadrature(nodes: ArrayF, weights: ArrayF, exact_to: int | _Inf | None = None)

The basic interface for a quadrature rule.

nodes: ArrayF

An array of shape (dim, nnodes), where dim is the dimension of the quadrature rule. In 1D, the shape is just (nnodes,).

weights: ArrayF

A vector of length nnodes that contains the quadrature weights.

property dim: int

Dimension of the space on which the quadrature rule applies.

property exact_to: int | _Inf

Summed polynomial degree up to which the quadrature is exact.

In higher-dimensions, the quadrature is supposed to be exact on (at least) \(P^N\), where \(N\) = exact_to. If the quadrature accuracy is not known, exact_to will not be set, and an AttributeError will be raised when attempting to access this information.

__call__(f: Callable[[ArrayF], ArrayF]) → ArrayF | floating

Evaluate the callable f at the quadrature nodes and return its integral.

f is assumed to accept arrays of shape (dims, npts), or of shape (npts,) for 1D quadrature.

class modepy.ZeroDimensionalQuadrature

Bases: Quadrature

A quadrature rule that should be used for 0d domains (i.e. points).

modepy.quadrature_for_space(space: FunctionSpace, shape: Shape) → Quadrature

modepy.quadrature_for_space(space: PN, shape: Simplex) → Quadrature

modepy.quadrature_for_space(space: TensorProductSpace, shape: TensorProductShape) → Quadrature

Returns:

a Quadrature that exactly integrates the functions in space over shape.

class modepy.quadrature._Inf

A sentinel type for infinite results. Do not use directly. Use isinf() instead.

modepy.quadrature.isinf(obj: object) → bool

Jacobi-Gauss quadrature in one dimension

class modepy.JacobiGaussQuadrature(alpha: float, beta: float, N: int, backend: str | None = None, force_dim_axis: bool = False)

Bases: Quadrature

A Gauss quadrature of order N associated with the Jacobi polynomials.

The quadrature rule can be used for weighted integrals of the form

\[I[f] = \int_{-1}^1 f(x) (1 - x)^\alpha (1 + x)^\beta\, \mathrm{d}x,\]

where \(\alpha, \beta > -1\). The quadrature rule is exact up to degree \(2N + 1\). The quadrature has N+1 nodes.

__init__(alpha: float, beta: float, N: int, backend: str | None = None, force_dim_axis: bool = False) → None

Parameters:

backend – Either "builtin" or "scipy". When the "builtin" backend is in use, there is an additional requirement that \(\alpha + \beta \ne -1\), with the exception of the Chebyshev quadrature \(\alpha = \beta = -1/2\). The "scipy" backend has no such restriction.

alpha: float

Power of \((1 - x)\) term in Jacobi quadrature.

beta: float

Power of \((1 + x)\) term in Jacobi quadrature.

static compute_weights_and_nodes(N: int, alpha: float, beta: float) → tuple[ArrayF, ArrayF]

Parameters:

• N – order of the Gauss quadrature (the order of exactly integrated polynomials is \(2 N + 1\)).

• alpha – power of \(1 - x\) in the Jacobi polynomial weight.

• beta – power of \(1 + x\) in the Jacobi polynomial weight.

Returns:

a tuple (nodes, weights) of quadrature notes and weights.

class modepy.LegendreGaussQuadrature(N: int, backend: str | None = None, force_dim_axis: bool = False)

Bases: JacobiGaussQuadrature

A Gauss quadrature rule with weight \(1\) and N+1 nodes.

Corresponds to a Gauss-Jacobi quadrature rule with \(\alpha = \beta = 0\).

class modepy.ChebyshevGaussQuadrature(N: int, kind: int = 1, backend: str | None = None, force_dim_axis: bool = False)

Bases: JacobiGaussQuadrature

A Gauss quadrature rule with weight \(\sqrt{1-x^2}^{\mp 1}\) and N+1 nodes.

The Chebyshev-Gauss quadrature rules of the first kind and second kind correspond to Gauss-Jacobi quadrature rules with \(\alpha = \beta = -0.5\) and \(\alpha = \beta = 0.5\), respectively.

Added in version 2019.1.

class modepy.GaussGegenbauerQuadrature(alpha: float, N: int, backend: str | None = None, force_dim_axis: bool = False)

Bases: JacobiGaussQuadrature

Gauss-Gegenbauer quadrature is a special case of Gauss-Jacobi quadrature with \(\alpha = \beta\) and N+1 nodes.

Added in version 2019.1.

modepy.quadrature.jacobi_gauss.jacobi_gauss_lobatto_nodes(alpha: float, beta: float, N: int, backend: str | None = None, force_dim_axis: bool = False) → ArrayF

Compute Gauss-Lobatto quadrature nodes associated with JacobiGaussQuadrature with the same parameters.

This helper returns only the N+1 nodes; the corresponding Gauss-Lobatto quadrature rule (using these nodes with appropriate weights) is exact for polynomials up to degree \(2N - 1\) when \(N \ge 1\).

For \(N = 0\), this function returns a single node at 0. This degenerate case does not correspond to a Gauss-Lobatto node set with endpoints \(\pm 1\), but is provided for convenience.

modepy.quadrature.jacobi_gauss.legendre_gauss_lobatto_nodes(N: int, backend: str | None = None, force_dim_axis: bool = False) → ArrayF

Compute the Legendre-Gauss-Lobatto quadrature nodes. N+1 is the number of nodes.

Exact to degree \(2N - 1\).

class modepy.JacobiGaussLobattoQuadrature(alpha: float, beta: float, N: int, *, backend: str | None = None, force_dim_axis: bool = False)

Compute the Jacobi-Gauss-Lobatto quadrature with respect to the Jacobi weight function

\[I[f] = \int_{-1}^1 f(x) (1 - x)^\alpha (1 + x)^\beta\, \mathrm{d}x,\]

There will be N+1 nodes. Exact to degree \(2N - 1\).

Added in version 2024.2.

class modepy.LegendreGaussLobattoQuadrature(N, *, backend: str | None = None, force_dim_axis: bool = False)

Clenshaw-Curtis and Fejér quadrature in one dimension

class modepy.ClenshawCurtisQuadrature(N: int, force_dim_axis: bool = False)

Bases: Quadrature

Clenshaw-Curtis quadrature of order N with N + 1 points.

The quadrature rule is exact up to degree \(N\) and can be nested. Its performance for differentiable functions is comparable with the classic Gauss-Legendre quadrature, which is exact for polynomials of degree up to \(2N + 1\). Implementation is based on [Waldvogel2003].

Integrates on the interval \((-1, 1)\).

[Waldvogel2003] (1,2)

Jörg Waldvogel, Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules, BIT Numerical Mathematics, 2003, Vol. 43, No. 1, pp. 001-018. DOI

class modepy.FejerQuadrature(N: int, kind: int = 1, force_dim_axis: bool = False)

Bases: Quadrature

Fejér quadrature rules of order N.

• Fejér quadrature of the first kind has N points and uses Chebyshev nodes, i.e. the roots of Chebyshev polynomials.

• Fejér quadrature of the second kind has N - 1 points and uses only the interior extrema of the Chebyshev nodes, i.e. the true stationary points. This rule is almost identical to Clenshaw-Curtis and can be nested.

Integrates on the interval \((-1, 1)\). Implementation is based on [Waldvogel2003].

kind: int

Kind of the Fejér quadrature, either first-kind or second-kind.

Gauss-Kronrod quadrature

class modepy.quadrature.kronrod.KronrodGaussQuadrature(nodes: ArrayF, weights: ArrayF, exact_to: int | _Inf | None = None)

class modepy.quadrature.kronrod.KronrodQuadrature(nodes: ArrayF, weights: ArrayF, weights_at_gauss_nodes: ArrayF, exact_to: int, gauss_quadrature: KronrodGaussQuadrature)

Uses a computation following [Laurie1997].

nodes: ArrayF

an array of shape (d, nnodes), where d is the dimension of the quadrature rule.

weights: ArrayF

an array of length nnodes.

weights_at_gauss_nodes: ArrayF

gauss_quadrature: KronrodGaussQuadrature

exact_to: int

__call__(f: Callable[[ArrayF], ArrayF]) → ArrayF | floating

Evaluate the Kronrod estimate for the integral.

Note

This will (re-)evaluate f at both the Gauss and the Kronrod nodes, likely negating any cost advantage of Kronrod. As a result, it is recommended to reimplement the evaluation logic, following the code of this function, while reusing the function values from the Gauss evaluation.

modepy.quadrature.kronrod.make_kronrod_quadrature(n: int) → KronrodQuadrature

References

[Laurie1997]

Dirk P. Laurie, “Calculation of Gauss-Kronrod quadrature rules,” Mathematics of Computation, vol. 66, no. 219, pp. 1133-1145 (1997). DOI

Transplanted quadrature in one dimension

The transplanted maps implemented here include the conformal maps from Hale-Trefethen (the sausage polynomial family and the strip map) as well as the earlier Kosloff-Tal-Ezer \(\arcsin\) map.

Note

In using the term ‘transplanted’, we are following the terminology from [HaleTrefethen2008]. In other nomenclature, this is also referred to as a change of variables transformation using a conformal mapping.

Given a base rule \((s_i, w_i^{(s)})\) on \([-1,1]\), transplanted quadrature uses a map \(g(s): [-1, 1] \to [-1, 1]\) to build

\[x_i = g(s_i), \qquad w_i = w_i^{(s)} g'(s_i),\]

so that

\[\int_{-1}^1 f(x)\,\mathrm{d}x = \int_{-1}^1 f(g(s))\,g'(s)\,\mathrm{d}s \approx \sum_i w_i f(x_i).\]

Map functions

Identity map

modepy.quadrature.transplanted.map_identity(s: ArrayF) → tuple[ArrayF, ArrayF]

Identity transplant map on \([-1, 1]\).

Returns:

(nodes, jacobian) where nodes is a copy of s and jacobian is an array of ones with the same shape.

Sausage polynomial maps

modepy.quadrature.transplanted.map_sausage(s: ArrayF, degree: int) → tuple[ArrayF, ArrayF]

Odd-degree polynomial sausage map from Hale-Trefethen [HaleTrefethen2008].

This is the normalized odd Taylor truncation of \(\arcsin(s)\) through the monomial of degree degree.

Parameters:

• s – quadrature nodes on \([-1, 1]\).

• degree – positive odd degree in {1, 3, 5, ...}.

Returns:

(nodes, jacobian).

Kosloff-Tal-Ezer map

modepy.quadrature.transplanted.map_kosloff_tal_ezer(s: ArrayF, *, rho: float = 1.4, alpha: float | None = None) → tuple[ArrayF, ArrayF]

Kosloff-Tal-Ezer map from [KosloffTalEzer1993].

The map is

\[g(s) = \frac{\arcsin(\alpha s)}{\arcsin(\alpha)},\]

where \(0 < \alpha < 1\). If alpha is not provided, it is chosen from rho using

\[\alpha = \frac{2}{\rho + \rho^{-1}},\]

matching the parameter choice discussed by Hale-Trefethen [HaleTrefethen2008] for a \(\rho\)-ellipse analyticity model.

Parameters:

• s – quadrature nodes on \([-1, 1]\).

• rho – ellipse parameter, must satisfy rho > 1. Ignored if alpha is given explicitly.

• alpha – map parameter satisfying 0 < alpha < 1. If None, computed from rho.

Returns:

(nodes, jacobian).

Strip conformal map

modepy.quadrature.transplanted.map_strip(s: ArrayF, *, rho: float = 1.4) → tuple[ArrayF, ArrayF]

Strip map from Hale-Trefethen [HaleTrefethen2008] transplanted quadrature.

The map is based on the Schwarz-Christoffel transformation that maps the unit disk to a rectangle, composed with an \(\arcsin\) to pull back to \([-1, 1]\). The parameter rho controls the half-width of the analyticity strip: a larger rho concentrates nodes near the endpoints.

Parameters:

• s – quadrature nodes on \((-1, 1)\) (strict interior).

• rho – strip parameter, must satisfy rho > 1.

Returns:

(nodes, jacobian).

Important

This map requires interior nodes (abs(s) < 1), so it is intended for base rules such as Legendre-Gauss that do not include endpoints. Other common rules, such as Gauss-Lobatto or Clenshaw-Curtis, cannot be used.

Quadrature wrappers

modepy.transplanted_1d_quadrature(quadrature: Quadrature, map_fn: Callable[[ArrayF], tuple[ArrayF, ArrayF]]) → Quadrature

Map an existing 1D quadrature rule using a transplant map.

The transformed rule approximates

\[\int_{-1}^1 f(x)\,\mathrm{d}x = \int_{-1}^1 f(g(s)) g'(s)\,\mathrm{d}s,\]

by mapping existing nodes \(s_i\) and scaling existing weights \(w_i\) with \(g'(s_i)\).

Parameters:

• quadrature – a one-dimensional Quadrature whose nodes lie in \([-1, 1]\).

• map_fn – a callable (s: ArrayF) -> (nodes, jacobian), such as map_identity(), map_sausage(), map_kosloff_tal_ezer(), or map_strip(), with parameters (if any) bound via functools.partial().

Returns:

a new Quadrature with mapped nodes and adjusted weights.

modepy.transplanted_legendre_gauss_quadrature(n: int, map_fn: Callable[[ArrayF], tuple[ArrayF, ArrayF]], *, backend: str | None = None, force_dim_axis: bool = False) → Quadrature

Legendre-Gauss quadrature transplanted by a Trefethen map.

Example

from functools import partial

import modepy as mp
from modepy.quadrature.transplanted import map_kosloff_tal_ezer, map_sausage

q_kte = mp.transplanted_legendre_gauss_quadrature(
    20,
    partial(map_kosloff_tal_ezer, rho=1.4),
    force_dim_axis=True,
)

q_sausage = mp.transplanted_legendre_gauss_quadrature(
    20,
    partial(map_sausage, degree=9),
    force_dim_axis=True,
)

References

[HaleTrefethen2008] (1,2,3,4)

N. Hale and L. N. Trefethen, New Quadrature Formulas from Conformal Maps, SIAM Journal on Numerical Analysis 46(2), 930-948 (2008), doi:10.1137/07068607X.

[KosloffTalEzer1993]

D. Kosloff and H. Tal-Ezer, A Modified Chebyshev Pseudospectral Method with an \(O(N^{-1})\) Time Step Restriction, Journal of Computational Physics 104(2), 457-469 (1993), doi:10.1006/jcph.1993.1044.

Quadratures on the simplex

exception modepy.QuadratureRuleUnavailable

Added in version 2013.3.

class modepy.GrundmannMoellerSimplexQuadrature(order: int, dimension: int)

Bases: Quadrature

Cubature on an n-simplex from [Grundmann1978].

This cubature rule has both negative and positive weights. It is exact for polynomials up to order \(2s + 1\), where \(s\) is given as order. The integration domain is the unit simplex (see Coordinates on the triangle and Coordinates on the tetrahedron).

[Grundmann1978]

A. Grundmann and H.M. Moeller, Invariant integration formulas for the n-simplex by combinatorial methods, SIAM Journal on Numerical Analysis, vol. 15, pp. 282–290, 1978. DOI

__init__(order: int, dimension: int) → None

Parameters:

• order – a parameter correlated with the total degree of polynomials that are integrated exactly (see also exact_to).

• dimension – the number of dimensions for the quadrature rule.

__call__(f: Callable[[ArrayF], ArrayF]) → ArrayF | floating

Evaluate the callable f at the quadrature nodes and return its integral.

f is assumed to accept arrays of shape (dims, npts), or of shape (npts,) for 1D quadrature.

class modepy.XiaoGimbutasSimplexQuadrature(order: int, dims: int)

Bases: Quadrature

A (nearly) Gaussian simplicial quadrature with very few quadrature nodes from [Xiao2010].

This rule is available for low-to-moderate orders. The integration domain is the unit simplex (see Coordinates on the triangle and Coordinates on the tetrahedron).

[Xiao2010]

H. Xiao and Z. Gimbutas, A numerical algorithm for the construction of efficient quadrature rules in two and higher dimensions, Computers & Mathematics with Applications, vol. 59, no. 2, pp. 663-676, 2010. DOI

__init__(order: int, dims: int) → None

Parameters:

• order – the total degree to which the quadrature rule is exact.

• dims – the number of dimensions for the quadrature rule. 2 for quadrature on triangles and 3 for tetrahedra.

Raises:

modepy.QuadratureRuleUnavailable if no quadrature rule for therequested parameters is available.

__call__(f: Callable[[ArrayF], ArrayF]) → ArrayF | floating

Evaluate the callable f at the quadrature nodes and return its integral.

f is assumed to accept arrays of shape (dims, npts), or of shape (npts,) for 1D quadrature.

class modepy.VioreanuRokhlinSimplexQuadrature(order: int, dims: int)

Bases: Quadrature

Simplicial quadratures with symmetric node sets and positive weights suitable for well-conditioned interpolation.

The integration domain is the unit simplex (see Coordinates on the triangle and Coordinates on the tetrahedron). When using these nodes, please acknowledge Zydrunas Gimbutas, who generated them as follows:

• The 2D nodes are based on the interpolation node set derived in the article [Vioreanu2011].

  Note that in Vioreanu’s tables, only orders 5, 6, 9, and 12 are rotationally symmetric, which gives one extra order for integration and better interpolation conditioning. Also note that since the tables have been re-generated independently, the nodes and weights may be different.

• The 3D nodes were derived from the modepy.warp_and_blend_nodes().

• A tightening algorithm was then applied, as described in [Vioreanu2012].

[Vioreanu2011]

B. Vioreanu and V. Rokhlin, Spectra of Multiplication Operators as a Numerical Tool, Yale CS Tech Report 1443. PDF

[Vioreanu2012]

B. Vioreanu, Spectra of Multiplication Operators as a Numerical Tool, Yale University, 2012. PDF

Added in version 2013.3.

__init__(order: int, dims: int) → None

Parameters:

• order – The total degree to which the quadrature rule is exact for interpolation.

• dims – The number of dimensions for the quadrature rule. 2 for quadrature on triangles and 3 for tetrahedra.

Raises:

modepy.QuadratureRuleUnavailable if no quadrature rule for the requested parameters is available.

__call__(f: Callable[[ArrayF], ArrayF]) → ArrayF | floating

Evaluate the callable f at the quadrature nodes and return its integral.

f is assumed to accept arrays of shape (dims, npts), or of shape (npts,) for 1D quadrature.

class modepy.JaskowiecSukumarQuadrature(order: int, dims: int)

Bases: Quadrature

Symmetric quadrature rules with positive weights for tetrahedra from [Jaskowiec2021].

The integration domain is the unit tetrahedron.

[Jaskowiec2021]

J. Jaśkowiec, N. Sukumar, High‐order Symmetric Cubature Rules for Tetrahedra and Pyramids, International Journal for Numerical Methods in Engineering, Vol. 122, pp. 148–171, 2021, DOI

__init__(order: int, dims: int) → None

Parameters:

• nodes – an array of shape (d, nnodes), where d is the dimension of the quadrature rule.

• weights – an array of length nnodes.

• exact_to – an optional argument denoting the summed polynomial degree to which the quadrature is exact. By default, exact_to is None and will not be set as an attribute.

__call__(f: Callable[[ArrayF], ArrayF]) → ArrayF | floating

Evaluate the callable f at the quadrature nodes and return its integral.

f is assumed to accept arrays of shape (dims, npts), or of shape (npts,) for 1D quadrature.

Quadratures on the hypercube

class modepy.WitherdenVincentQuadrature(order: int, dims: int)

Bases: Quadrature

Symmetric quadrature rules with positive weights for rectangles and hexahedra from [Witherden2015].

The integration domain is the unit hypercube \([-1, 1]^d\), where \(d\) is the dimension.

[Witherden2015]

F. D. Witherden, P. E. Vincent, On the Identification of Symmetric Quadrature Rules for Finite Element Methods, Computers & Mathematics with Applications, Vol. 69, pp. 1232–1241, 2015. DOI

__init__(order: int, dims: int) → None

Parameters:

• nodes – an array of shape (d, nnodes), where d is the dimension of the quadrature rule.

• weights – an array of length nnodes.

• exact_to – an optional argument denoting the summed polynomial degree to which the quadrature is exact. By default, exact_to is None and will not be set as an attribute.

__call__(f: Callable[[ArrayF], ArrayF]) → ArrayF | floating

Evaluate the callable f at the quadrature nodes and return its integral.

f is assumed to accept arrays of shape (dims, npts), or of shape (npts,) for 1D quadrature.

class modepy.TensorProductQuadrature(quads: Iterable[Quadrature])

Bases: Quadrature

A tensor product quadrature of one-dimensional Quadratures.

quadratures: Sequence[Quadrature]

The lower-dimensional quadratures from which the tensor product quadrature is composed.

__init__(quads: Iterable[Quadrature]) → None

Parameters:

quad – a iterable of Quadrature objects for one-dimensional intervals, one for each dimension of the tensor product.

class modepy.LegendreGaussTensorProductQuadrature(N: int, dims: int, backend: str | None = None)

Bases: TensorProductQuadrature

A tensor product using only LegendreGaussQuadrature one-dimensional rules.