Support for construction of new quadrature rules

Added in version 2025.2.

modepy.quadrature.construction.Integrand

alias of Callable[[ndarray], ndarray]

class modepy.quadrature.construction.LinearCombinationIntegrand(coefficients: ndarray, functions: Sequence[Callable[[ndarray], ndarray]])

Note

New linear combinations should be created by linearly_combine(), which flattens nested combinations by taking advantage of associativity, substantially decreasing evaluation cost.

coefficients: ndarray

functions: Sequence[Callable[[ndarray], ndarray]]

__call__(points: ndarray) → ndarray

Evaluate the linear combination of functions with coefficients at points.

modepy.quadrature.construction.linearly_combine(coefficients: ndarray, functions: Sequence[Callable[[ndarray], ndarray]]) → LinearCombinationIntegrand

Returns a linear combination of functions with coefficients. If functions are themselves LinearCombinationIntegrands, and all functions use a consistent basis, then the resulting LinearCombinationIntegrand will be flattened via associativity, leading to faster execution time.

modepy.quadrature.construction.orthogonalize_basis(integrate: Callable[[Callable[[ndarray], ndarray]], inexact], basis: Sequence[Callable[[ndarray], ndarray]]) → Sequence[LinearCombinationIntegrand]

Let \(\Omega\subset\mathbb C^n\) be a domain. (Domains over the reals are allowable as well.) Returns linear combinations of functions in basis that is orthogonal under the (complex-valued) \(L^2\) inner product induced by integrate.

Parameters:

• integrate – Computes an integral of the passed integrand over \(\Omega\). Must use integration nodes compatible with basis.

• basis – A sequence of functions that accept an array of nodes of shape either (n, nnodes) or (nnodes,) and return an array of shape (nnodes,) with the value of the basis function at the node.

modepy.quadrature.construction.adapt_2d_integrands_to_complex_arg(functions: Sequence[Callable[[ndarray], ndarray]]) → Sequence[Callable[[ndarray], ndarray]]

Converts a list of Integrands taking nodes in real-valued arrays of shape (2, nnodes) to ones accepting a single complex-valued array of shape (nnodes,). See guess_nodes_vioreanu_rokhlin() for the main intended use case.

modepy.quadrature.construction.guess_nodes_vioreanu_rokhlin(integrate: Callable[[Callable[[ndarray], ndarray]], inexact], onb: Sequence[Callable[[ndarray], ndarray]]) → ndarray

Let \(\Omega\subset\mathbb C\) be a convex domain. Finds interpolation nodes for \(\Omega\) based on the multiplication-operator technique in [Vioreanu2011]. May be useful as an initial guess for a Gauss-Newton process to find new quadrature rules, as driven by quad_residual_and_jacobian().

Parameters:

• integrate – Must accurately integrate a product of two functions from onb and a degree-1 monomial over \(\Omega\).

• onb –

  An orthonormal basis of functions. Each function takes in an array of (complex) nodes and returns an array of the same shape containing the function values. Functions are expected to be orthogonal with respect to the complex-valued inner product.

  Integrands accepting real-valued node coordinates in two dimensions (shaped (2, nnodes) may be converted to functions acceptable by this function via adapt_2d_integrands_to_complex_arg().)

Returns:

An array of shape (len(onb),) containing nodes, complex-valued.

Note

This function is based on complex arithmetic and thus only usable for domains $Omega$ with one and two (real-valued) dimensions.

Note

(Empirically) it seems acceptable for the functions in onb to be exclusively real-valued, i.e., particularly, no assumptions on complex differentiability are needed.

Note

As noted in [Vioreanu2011], the returned nodes should obey the symmetry of the domain.

modepy.quadrature.construction.find_weights_undetermined_coefficients(integrands: Sequence[Callable[[ndarray], ndarray]], nodes: ndarray, reference_integrals: ndarray) → ndarray

Uses the method of undetermined coefficients to find weights for a quadrature rule using nodes, for the provided integrands.

Parameters:

• integrands – A sequence of functions that accept an array of nodes of shape either (ndim, nnodes) or (nnodes,) and return an array of shape (nnodes,) with the value of the basis function at the node.

• nodes – An array with shape (ndim, nnodes) or (nnodes,), real-valued. Must be compatible with integrands.

• reference_integrals – An array with shape (len(integrands),)

Note

This function also supports overdetermined systems. In that case, it will provide the least squares solution.

Note

Unlike guess_nodes_vioreanu_rokhlin(), domains of any dimensionality are allowed.

class modepy.quadrature.construction.QuadratureResidualJacobian(residual: ndarray, dresid_dweights: ndarray, dresid_dnodes: ndarray)

Contains residual value and Jacobian for a quadrature residual, see quad_residual_and_jacobian(). May be reused for residuals that include additional terms, such as those that penalize asymmetry of the noddes.

residual: ndarray

Shaped (nintegrands,).

dresid_dweights: ndarray

Shaped (nintegrands, nnodes).

dresid_dnodes: ndarray

Shaped (nintegrands, ndim*nnodes).

modepy.quadrature.construction.quad_residual_and_jacobian(nodes: ndarray, weights: ndarray, integrands: Sequence[Callable[[ndarray], ndarray]], integrand_derivatives: Sequence[Sequence[Callable[[ndarray], ndarray]]], reference_integrals: ndarray) → QuadratureResidualJacobian

Computes the residual and Jacobian of the objective function

\[\begin{split}\begin{bmatrix} \sum_{j=0}^{\text{nnodes-1}} \psi_0(\boldsymbol x_j) w_j - I_0\\ \vdots\\ \sum_{j=0}^{\text{nnodes-1}} \psi_{\text{nintegrands-1}} (\boldsymbol x_j) w_j - I_{\text{nintegrands}-1} \end{bmatrix}.\end{split}\]

Typically used with quad_gauss_newton_increment() to drive an iterative process for finding nodes and weights of a quadrature rule.

An initial guess for the nodes may be obtained via guess_nodes_vioreanu_rokhlin(), and an initial guess for the weights (given nodes) may be found via find_weights_undetermined_coefficients().

Note

Unlike guess_nodes_vioreanu_rokhlin(), domains of any dimensionality are allowed.

Parameters:

• nodes – An array with shape (ndim, nnodes) or (nnodes,), real-valued. Must be compatible with integrands.

• weights – An array with shape (nnodes,), real-valued

• integrands – A sequence of functions that accept an array of nodes of shape either (ndim, nnodes) or (nnodes,) and return an array of shape (nnodes,) with the value of the basis function at the node.

• integrand_derivatives – Derivatives of integrands along the coordinate axes, one list of functions per axis, with each sub-list matching integrands in structure.

• reference_integrals – Integrals of integrands, shaped (len(integrands),).

modepy.quadrature.construction.quad_gauss_newton_increment(qrj: QuadratureResidualJacobian) → tuple[ndarray, ndarray]

Return the Gauss-Newton increment based on the residual and Jacobian (see quad_residual_and_jacobian()), separated into the weight increment and the nodes increment, in the two entries of the returned tuple. The nodes increment has shape (ndim, nnodes).

Note

[Vioreanu2011] suggests that step size control should be used in the Gauss-Newton process. No notion of step size control is included in the returned increments.