archimedes.measure.stieltjes_recurrence¶

archimedes.measure.stieltjes_recurrence(

weight: Callable[[ndarray], ndarray],

support: tuple[float, float],

n: int,

**quad_kwargs,

) → tuple[ndarray, ndarray]¶

Monic recurrence coefficients for an arbitrary weight function.

Computes the coefficients of the monic three-term recurrence

\[\pi_{k+1}(x) = (x - \alpha_k) \, \pi_k(x) - \beta_k \, \pi_{k-1}(x), \qquad \pi_{-1} = 0, \ \pi_0 = 1\]

satisfied by the polynomials orthogonal with respect to weight on support, via the (discretized) Stieltjes procedure [1]:

\[\alpha_k = \frac{\langle x \pi_k, \pi_k \rangle} {\langle \pi_k, \pi_k \rangle}, \qquad \beta_k = \frac{\langle \pi_k, \pi_k \rangle} {\langle \pi_{k-1}, \pi_{k-1} \rangle}, \quad k \geq 1\]

where \(\langle f, g \rangle = \int_{\text{support}} f(x) \, g(x) \, w(x) \, dx\). The inner products are evaluated directly by adaptive quadrature (scipy.integrate.quad()).

This is the default implementation of recurrence_coeffs(); any Measure subclass that supplies weight and support (i.e. domain) gets a working recurrence_coeffs – and hence, via Golub-Welsch (golub_welsch_rule()), a Gauss quadrature rule – without deriving a closed-form recursion.

Parameters:

weight (callable) – Weight function \(w(x)\), evaluated at x.
support (tuple of float) – Integration bounds (a, b); either may be infinite.
n (int) – Number of coefficients to compute, i.e. degrees 0, ..., n - 1. Must be >= 1; not validated here.
**quad_kwargs – Forwarded to every internal scipy.integrate.quad() call (e.g. epsabs, epsrel, limit) – tightening these can extend the usable range of n for weights singular at an endpoint (see Notes).

Returns:

alpha, beta – Monic recurrence coefficients. beta[0] is the zeroth moment (total mass) of weight rather than a recursion coefficient – see recurrence_coeffs().

Return type:

ndarray, shape (n,)

Notes

Accuracy is limited by the compounding of each step’s adaptive quadrature error into the next, and degrades with n – how quickly depends on the weight. For smooth, bounded weights (e.g. Legendre-like), coefficients are accurate to near machine precision through about \(n \sim 15\). For weights singular at an endpoint (e.g. Jacobi-like), accuracy degrades starting around \(n \sim 8\). For larger n or better accuracy, either supply a closed-form recurrence_coeffs override, or pass tighter epsabs/epsrel/limit via **quad_kwargs as a partial mitigation.

References