Generating 1D recurrence coefficients for a given weight

orthopy contains four different algorithms that (in exact arithmetic) all do the same thing: They generate the recurrence coefficients for a 1D interval with a given weight function. Let's run this example with the weight function x ** 2 over the interval (-1, 1).

Stieltjes

import orthopy
import sympy

N = 10

alpha, beta, int_1 = orthopy.tools.stieltjes(
    # You need to provide a function that integrates a given function over the domain
    lambda x, fx: sympy.integrate(fx * x ** 2, (x, -1, 1)),
    N
)
print(alpha)
print(beta)

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[None, 3/5, 4/35, 25/63, 16/99, 49/143, 12/65, 27/85, 64/323, 121/399]

Golub-Welsch

(Only works in floating-point arithmetic.) For this, you first need to compute the first 2*N moments of your measure

integral(w(x) x^k dx)

import orthopy
import sympy

N = 10

x = sympy.Symbol("x")
moments = [sympy.integrate(x ** (2 + k), (x, -1, +1)) for k in range(2 * N + 1)]
moments = [float(m) for m in moments]  # convert to floats
print(moments)
print()

alpha, beta, int_1 = orthopy.tools.golub_welsch(moments)
print(alpha)
print(beta)

[0.6666666666666666, 0.0, 0.4, 0.0, 0.2857142857142857, 0.0, 0.2222222222222222, 0.0, 0.18181818181818182, 0.0, 0.15384615384615385, 0.0, 
0.13333333333333333, 0.0, 0.11764705882352941, 0.0, 0.10526315789473684, 0.0, 0.09523809523809523, 0.0, 0.08695652173913043]

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[nan 0.6        0.11428571 0.3968254  0.16161616 0.34265734
 0.18461538 0.31764706 0.19814241 0.30325815]

Chebyshev

The moments can be used in exact arithmetic here.

import orthopy
import sympy

N = 10

x = sympy.Symbol("x")
moments = [sympy.integrate(x ** (2 + k), (x, -1, +1)) for k in range(2 * N)]
print(moments)
print()

alpha, beta, int_1 = orthopy.tools.chebyshev(moments)
print(alpha)
print(beta)

[2/3, 0, 2/5, 0, 2/7, 0, 2/9, 0, 2/11, 0, 2/13, 0, 2/15, 0, 2/17, 0, 2/19, 0, 2/21, 0]

[0 0 0 0 0 0 0 0 0 0]
[nan 3/5 4/35 25/63 16/99 49/143 12/65 27/85 64/323 121/399]

Modified Chebyshev

For this, you first need to compute the first 2*N modified moments of your measure

integral(w(x) p_k(x) dx)

with a particular set of orthogonal polynomials p_k. A common choice are the Legendre polynomials.

import orthopy
import sympy

N = 10


x = sympy.Symbol("x")
evaluator = orthopy.c1.legendre.Eval(x, "monic")
legendre_polys = [sympy.expand(next(evaluator)) for _ in range(2 * N)]

mod_moments = [sympy.integrate(x ** 2 * poly, (x, -1, +1)) for poly in legendre_polys]
print(mod_moments)
print()

rc = orthopy.c1.legendre.RecurrenceCoefficients("monic", symbolic=True)
alpha, beta, int_1 = orthopy.tools.chebyshev_modified(mod_moments, rc)
print(alpha)
print(beta)

[2/3, 0, 8/45, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[nan, 3/5, 4/35, 25/63, 16/99, 49/143, 12/65, 27/85, 64/323, 121/399]

Note: If you perform the computation in floating-point arithmetic for example because you need to compute the scheme fast), it makes sense to think about the numerical implications here. That's because the map from the moments to the recurrence coefficients can be very ill-conditioned, meaning that small round-off errors can lead to very wrong coefficients. For further computation, it's numerically beneficial if the moments are either 0 or in the same order of magnitude. That is why the flexibility provided by the modified Chebyshev scheme can be important. (Note how almost all moments are 0 there.)