Skip to content

Creating your own Gauss quadrature in two simple steps

Jump to bottom
Nico Schlömer edited this page Sep 29, 2021 · 12 revisions

You have a measure (or, more colloquially speaking, a domain and a nonnegative weight function) and would like to generate the matching Gauss quadrature? Great, here's how to do it.

As an example, let's try and generate the Gauss quadrature with 10 points for the weight function x^2 on the interval [-1, +1].

The package orthopy (which computes orthogonal polynomials) helps us here.

TLDR:

import orthopy
import quadpy
import sympy


N = 10

# int_{-1}^{+1} x^2 x^k dx
moments = [sympy.S(1 + (-1) ** kk) / (kk + 3) for kk in range(2 * N)]

alpha, beta, int_1 = orthopy.tools.chebyshev(moments)
alpha = [float(sympy.N(a)) for a in alpha]
beta = [float(sympy.N(b)) for b in beta]
points, weights = quadpy.tools.scheme_from_rc(alpha, beta, int_1)

Some explanations:

  1. You first need to compute the recurrence coefficients of the corresponding orthogonal polynomials. There are four ways to do that. They all yield the same result when computing in exact arithmetic (which is the default).

    • Stieltjes.

      import orthopy
import sympy

N = 10

alpha, beta, int_1 = orthopy.tools.stieltjes(
    lambda t, ft: sympy.integrate(t ** 2 * ft, (t, -1, 1)), N
)
print(alpha)
print(beta)
print(int_1)

      [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]
2/3

    • Golub-Welsch. (Only works in floating-point arithmetic.) For this, you first need to compute the first 2*N + 1 moments of your measure $$ I_k = \int_{-1}^{1} x^2 x^k\, dx = \frac{1 + (-1)^k}{k + 3} $$ (Rendered with Purple Pi.)

      import orthopy

N = 10

moments = [(1 + (-1) ** kk) / (kk + 3) for kk in range(2 * N + 1)]
print(moments)
print()

alpha, beta = 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]
0.816496580927726

    • Chebyshev. The moments can be used in exact arithmetic here.

      import orthopy
import sympy

N = 10

moments = [sympy.S(1 + (-1) ** kk) / (kk + 3) for kk in range(2 * N)]
print(moments)
print()

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

      [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]
[None, 3/5, 4/35, 25/63, 16/99, 49/143, 12/65, 27/85, 64/323, 121/399]
2/3

    • Modified Chebyshev. For this, you first need to compute the first 2*N modified moments of your measure $$ I_k = \int_{-1}^{1} x^2 p_k(x)\, dx $$ with a particular set of orthogonal polynomials $p_k$. A common choice are the Legendre polynomials.

      import orthopy
import quadpy
import sympy

N = 10


x = sympy.Symbol("x")
evaluator = orthopy.c1.legendre.Eval(x, "monic")
moments = []
for _ in range(2 * N):
    leg = next(evaluator)
    val = sympy.integrate(x ** 2 * leg, (x, -1, 1))
    moments.append(val)

print(moments)
print()

rc = orthopy.c1.legendre.RecurrenceCoefficients("monic", symbolic=True)
rc = [rc[k] for k in range(2 * N)]

alpha, beta, int_1 = orthopy.tools.chebyshev_modified(moments, rc)

print(alpha)
print(beta)
print(int_1)

      [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]
[None 3/5 4/35 25/63 16/99 49/143 12/65 27/85 64/323 121/399]
2/3

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 an unusable scheme. 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.)

  1. From the recurrence coefficients, we generate the Gauss points and weights. Since symbolic computation can take very long even for small sizes, we convert alpha and beta to numpy float arrays first. (If you need more digits, look at mpmath arrays.) 
    alpha = [float(sympy.N(a)) for a in alpha]
beta = [float(sympy.N(b)) for b in beta]
int_1 = float(int_1)

points, weights = quadpy.tools.scheme_from_rc(alpha, beta, int_1)
print(points)
print(weights)
    [-0.97822866 -0.8870626  -0.73015201 -0.51909613 -0.26954316  0.26954316 0.51909613  0.73015201  0.8870626   0.97822866]
[0.0532709947237133 0.0988166881454095 0.0993154007474127 0.0628365763465905 0.0190936733702069 0.0190936733702070 0.0628365763465907 0.0993154007474143 0.0988166881454080 0.0532709947237134]
    Congratulations! Your Gaussian quadrature rule.

Related tools

  • Transforming Gaussian points and weights back to recurrence coefficients:

    alpha, beta = quadpy.tools.coefficients_from_gauss(points, weights)

  • The Gautschi test: As recommended by Gautschi, you can test your moment-based scheme with

    err = orthopy.tools.gautschi_test_3(moments, alpha, beta)

Relevant publications

Clone this wiki locally