Skip to content

Creating your own Gauss quadrature in two simple steps

Jump to bottom
Nico Schlömer edited this page Feb 14, 2020 · 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].

TLDR:

import quadpy
import sympy


N = 10
moments = quadpy.tools.integrate(lambda x: [x ** (2 + k) for k in range(2 * N)], -1, +1)
alpha, beta = quadpy.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)

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 quadpy

N = 10

alpha, beta = quadpy.tools.stieltjes(lambda x: x ** 2, -1, 1, N)
print(alpha)
print(beta)
      [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[2/3, 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 quadpy

N = 10

moments = quadpy.tools.integrate(lambda x: [x ** (2 + k) for k in range(2 * N + 1)], -1, +1)
moments = [float(m) for m in moments]  # convert to floats
print(moments)
print()

alpha, beta = quadpy.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.]
[0.66666667 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 quadpy

N = 10

moments = quadpy.tools.integrate(lambda x: [x ** (2 + k) for k in range(2 * N)], -1, +1)
print(moments)
print()

alpha, beta = quadpy.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]
[2/3, 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, here taken from orthopy. 
      import orthopy
import quadpy

N = 10


def leg_polys(x):
    return orthopy.line_segment.tree_legendre(x, 2 * N - 1, "monic", symbolic=True)


moments = quadpy.tools.integrate(
    lambda x: [x ** 2 * leg_poly for leg_poly in leg_polys(x)], -1, +1
)
print(moments)
print()

_, _, a, b = orthopy.line_segment.recurrence_coefficients.legendre(
    2 * N, "monic", symbolic=True
)
alpha, beta = quadpy.tools.chebyshev_modified(moments, a, b)
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]
[2/3 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 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]

points, weights = quadpy.tools.scheme_from_rc(alpha, beta)
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.05327099 0.09881669 0.0993154  0.06283658 0.01909367 0.01909367
 0.06283658 0.0993154  0.09881669 0.05327099]
    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 = quadpy.tools.check_coefficients(moments, alpha, beta)

Relevant publications

Clone this wiki locally