{ } Raw JSON

bundles / scipy latest / scipy / integrate / _rules / _base / Rule

class

scipy.integrate._rules._base:Rule

source: /scipy/integrate/_rules/_base.py :6

Signature

class   Rule ( )

Members

Summary

Base class for numerical integration algorithms (cubatures).

Extended Summary

Finds an estimate for the integral of f over the region described by two arrays a and b via estimate, and find an estimate for the error of this approximation via estimate_error.

If a subclass does not implement its own estimate_error, then it will use a default error estimate based on the difference between the estimate over the whole region and the sum of estimates over that region divided into 2^ndim subregions.

Examples

In the following, a custom rule is created which uses 3D Genz-Malik cubature for the estimate of the integral, and the difference between this estimate and a less accurate estimate using 5-node Gauss-Legendre quadrature as an estimate for the error.
import numpy as np
from scipy.integrate import cubature
from scipy.integrate._rules import (
    Rule, ProductNestedFixed, GenzMalikCubature, GaussLegendreQuadrature
)
def f(x, r, alphas):
    # f(x) = cos(2*pi*r + alpha @ x)
    # Need to allow r and alphas to be arbitrary shape
    npoints, ndim = x.shape[0], x.shape[-1]
    alphas_reshaped = alphas[np.newaxis, :]
    x_reshaped = x.reshape(npoints, *([1]*(len(alphas.shape) - 1)), ndim)
    return np.cos(2*np.pi*r + np.sum(alphas_reshaped * x_reshaped, axis=-1))
genz = GenzMalikCubature(ndim=3)
gauss = GaussKronrodQuadrature(npoints=21)
gauss_3d = ProductNestedFixed([gauss, gauss, gauss])
class CustomRule(Rule):
    def estimate(self, f, a, b, args=()):
        return genz.estimate(f, a, b, args)
    def estimate_error(self, f, a, b, args=()):
        return np.abs(
            genz.estimate(f, a, b, args)
            - gauss_3d.estimate(f, a, b, args)
        )
rng = np.random.default_rng()
res = cubature(
    f=f,
    a=np.array([0, 0, 0]),
    b=np.array([1, 1, 1]),
    rule=CustomRule(),
    args=(rng.random((2,)), rng.random((3, 2, 3)))
)
res.estimate

See also

FixedRule

Aliases

  • scipy.integrate._rules.Rule