Gaussian Quadrature

Posted on by ComPhys

Gaussian quadrature is a technique to evaluate integrals by summing a weighted integrand at a finite number of points. It can correctly evaluate some forms of integral to machine precision, rather than the approximations produced by techniques like the trapezium rule. When I first looked into using it, I found a huge number of explanations but it still took a while for me to put all the pieces together. In this post I will explain the background of the method in the way I found most helpful, and give a few demonstrative examples. I’ll show how to integrate polynomial functions between any finite limits, and give an example with an infinite upper limit. Before that though, I will introduce the more general quadrature method.

Integration by Quadrature

Suppose we try to evaluate an integral numerically by expanding it as a summation

\int_{a}^{b} f(x) \,\textrm{d}x = w_0 f(x_0) + w_1 f(x_1) + w_2 f(x_2) + ... + w_n f(x_n)

For this to work the coefficients (“weights”) w_n and sample points x_n have to be found. Weights can be found for any set of sample points; to illustrate this suppose there is an integral to evaluate between limits [-1, 1] and we are happy to approximate f(x) as an n-term polynomial series, with n=3 so that f(x) \approx a_0 + a_1 x + a_2 x^2. I might chose the sample points (x_0, x_1, x_2) = (-1, 0, 1) then we can write

\begin{aligned} \int_{-1}^{1} f(x) \,\textrm{d} x =& w_0 \left( a_0 + a_1 x_0 + a_2 x_0^2 \right) +\\ & w_1 \left( a_0 + a_1 x_1 + a_2 x_1^2 \right) +\\ & w_2 \left( a_0 + a_1 x_2 + a_2 x_2^2 \right) \\ = & \int_{-1}^{1} a_0 + a_1 x + a_2 x^2 \,\textrm{d} x  = \left. a_0 x \right|_{-1}^{1} \left. a_1 \frac{1}{2}x^2 \right|_{-1}^{1} \left. a_2 \frac{1}{3}x^3 \right|_{-1}^{1} = 2 a_0 + 0 a_1 + \frac{2}{3} a_2, \\ \end{aligned}

where the last line was just finding the integral of the polynomial. This would certainly be satisfied if the coefficients of a_i were equated, so \left( w_0 a_0 + w_1 a_0 + w_2 a_0 \right) = 2 a_0, \left( w_0 a_1 x_0 + w_1 a_1 x_1 + w_2 a_1 x_2 \right) = 0 and \left( w_0 a_2 x_0^2 + w_1 a_2 x_1^2 + w_2 a_2 x_2^2 \right) = \frac{2}{3}a_2, so we can cancel the a‘s and write the matrix equation

\begin{pmatrix} 1 & 1 & 1 \\ x_0 & x_1 & x_2 \\ x_0^2 & x_1^2 & x_2^2 \end{pmatrix} \begin{pmatrix} w_0 \\ w_1 \\ w_2 \end{pmatrix} = \begin{pmatrix} 2 \\ 0 \\ \frac{2}{3} \end{pmatrix}.

After substituting for the sample points x_i, this linear set of equations can be solved using linear algebra techniques to get w_0 = \frac{1}{3}, w_1 = \frac{4}{3} and w_2 = \frac{1}{3}.

Using these weights let’s take a specific example of f(x) = 5x^2 + 3x + 2. It is easy to calculate f(-1) = 4, f(0) = 2 and f(1) = 10, and so the result of the integral is \frac{1}{3} 4 + \frac{4}{3} 2 + \frac{1}{3} 10 = \frac{22}{3}. Evaluating the integral directly gives \left. \frac{5 x^3}{3} + \frac{3 x^2}{2} + 2x \right|_{-1}^{1} = \frac{22}{3}, so we had *exactly* the right answer with only three sample points. For comparison using the trapezium rule with 1000 sample points gives 17.333360 which is only accurate to four decimal places.

With this technique polynomials up to degree n-1 can be integrated exactly between given limits from n sample points. If a smaller number of sample points are used, then the result will be an approximation from the series expansion between the limits; this could be a good approximation but there is no guarantee. Functions which do not have finite polynomial expansions such as cos(x) cannot be integrated exactly, but a high-degree polynomial can be an excellent approximation in a range.

Gaussian Quadrature – Choosing Sample Points

I said earlier that any sample points could be used if appropriate weights were calculated. In the previous example, if we had chosen (x_0, x_1, x_2) = (0, 0.5, 1) then the weights would have been (w_0, w_1, w_2) = (10/3, -8/3, 4/3) but the answer would have been identical. One could ask the question: “is there an optimum set of sample locations?” It turns out that for certain classes of functions there is, and in these cases n sample points can produce exact answers of integrands with degree 2n - 1. This relies on some of the fundamental properties of orthogonal polynomials, and is what distinguishes “quadrature rules” from Gaussian-quadrature.

Suppose the function to be integrated can be written as

\begin{aligned} \int_{a}^{b} w(x) f(x) \,\textrm{d} x, \end{aligned}

where w(x) is referred to as a weight function. Rather than expanding f(x) as a polynomial series with n terms, expand it up to 2n terms then use polynomial division to write it with the form

\begin{aligned} f(x) = q(x) G_{n}(x) + r(x) \end{aligned}

where q(x) is a quotient function, r(x) is a remainder function and G_{n}(x) is an as-yet unspecified orthogonal polynomial of degree n [1]. Whilst q(x) and r(x) are unspecified their degrees must both be \leq n-1 [1] since f(x) has a degree of 2n -1.

Writing the integral out again

\begin{aligned} \int_{a}^{b} w(x) f(x) \,\textrm{d} x, = \int_{a}^{b} w(x) q(x) G_{n}(x) \,\textrm{d} x, + \int_{a}^{b} w(x) r(x) \,\textrm{d} x \end{aligned}

it is apparent that if q(x) were itself expanded in orthogonal polynomials it would only need to be expanded up to degree n-1, so q(x) = \sum_{i=0}^{n-1} b_i G_{i} and

\begin{aligned} \int_{a}^{b} w(x) f(x) \,\textrm{d} x, = \sum_{i=0}^{n-1} b_i \int_{a}^{b} w(x) G_i(x) G_{n}(x) \,\textrm{d} x, + \int_{a}^{b} w(x) r(x) \,\textrm{d} x. \end{aligned}

This is a very interesting form. From here, the properties of orthogonal polynomials imply that if the limits a and b were the limits of orthogonality for a class of orthogonal polynomials and w(x) happened to be their weight function, then the first integral would be zero. In that case the integral of a function of degree 2n - 1 could be evaluated exactly by performing the integral of r(x) which is only of degree \leq n-1.

Now, when we evaluate this integral numerically we don’t have the form of q(x) or r(x) and we certainly don’t want to have to work them out! What we have is a summation over n sample points

\begin{aligned} \int_{a}^{b} w(x) f(x) \,\textrm{d}x = & w_0 \left( \sum_{i=0}^{n-1} b_i w(x_0) G_i(x_0) G_{n}(x_0) + w(x_0) r(x_0) \right) +\\ & w_1 \left( \sum_{i=0}^{n-1} b_i w(x_1) G_i(x_1) G_{n}(x_1) + w(x_1) r(x_1) \right) +\\ & w_2 \left( \sum_{i=0}^{n-1} b_i w(x_2) G_i(x_2) G_{n}(x_2) + w(x_2) r(x_2) \right) + ... +\\ & w_{n-1} \left( \sum_{i=0}^{n-1} b_i w(x_{n-1}) G_i(x_{n-1}) G_{n}(x_{n-1}) + w(x_{n-1}) r(x_{n-1}) \right) + ...\\ \end{aligned}

Although the integral \sum_{i=0}^{n-1} b_i \int_{a}^{b} w(x) G_i(x) G_{n}(x) \,\textrm{d} x =0, there is nothing to guarantee that this is also true for the summations w_m \sum_{i=0}^{n} w(x_m) G_i(x_m) G_{n}(x_m). To ensure that they are zero, use the trick of making the sample positions x_m where G_n(x_m)=0; these points are available for a large number of polynomials in software packages and tables. The first term on every line then disappears, making the integral

\begin{aligned} \int_{a}^{b} w(x) f(x) \,\textrm{d}x = w_0 w(x_0) r(x_0) + w_1 w(x_1) r(x_1) + w_2 w(x_2) r(x_2) + ... + w_{n-1} w(x_{n-1}) r(x_{n-1}), \end{aligned}

which is n samples of a function with degree n-1, and we know from before that this situation produces an exact answer. To summarise, the function f(n) was expanded up to degree 2n-1. By carefully selecting sample locations, we have found an integral of an (n-1)-degree function which will give the same answer, and we only need n samples to evaluate that integral exactly. So we can evaluate a polynomial of degree 2n - 1 with only n samples of it. Since the zeros and weights only have to be calculated once (and that’s already been done for a lot of polynomials), this method allows fast and accurate evaluation of integrals with very few samples.

Example 1

The Legendre polynomials P_n(x) have the weight function w(x) = 1 and are orthogonal in the range [-1, +1]. The zeros of the third degree Legendre polynomial are (x_0, x_1, x_2) = (-\sqrt{\frac{3}{5}}, 0 , \sqrt{\frac{3}{5}}). These can be substituted into the matrix above giving

\begin{pmatrix} 1 & 1 & 1 \\ -\sqrt{\frac{3}{5}} & 0 & \sqrt{\frac{3}{5}} \\ \frac{3}{5} & 0 & \frac{3}{5} \end{pmatrix} \begin{pmatrix} w_0 \\ w_1 \\ w_2 \end{pmatrix} = \begin{pmatrix} 2 \\ 0 \\ \frac{2}{3} \end{pmatrix},

which can be solved to give the weights (w_0, w_1, w_2) = (\frac{5}{9}, \frac{8}{9}, \frac{5}{9}). The code to achieve this in Python is below; note that the SciPy function roots_legendre automatically returns the weights, so the code calculates them and prints them both out.


import numpy as np
import scipy.special as sp

roots, weights = sp.roots_legendre(3)

a = np.ones((len(roots), len(roots)))
for i in range(1, len(roots)):
    a[i] = roots**i

b = np.array([2, 0, 2/3])

my_weights = np.linalg.solve(a, b)
print(f"Calculated weights: {my_weights}")
print(f"Weights from Scipy: {weights}")

Using these sample points and weights, the function f(x) = 5x^4 + 3x + 2 evaluates to the exact answer of 6. This is expected because with n=3 sample points, results will be exact up to degree 2n-1 = 5. In contrast, if the sample points were (0, 0.5, 1) again then exact results would only be expected up to degree n=3; using them gives the answer 9.833 which is clearly wrong.

Making n samples of the function at the zeros of the Legendre polynomials and using the corresponding weights has increased the degree of the function which can be integrated exactly to 2n-1.

Example 2 – Change of Interval

Requiring limits [-1, +1] may seem restrictive, but the technique can be used for any finite interval by making a change of variables. Define

\begin{aligned} x &= \frac{b-a}{2} \xi + \frac{a+b}{2} \\ \frac{\textrm{d} x}{\textrm{d} \xi} &= \frac{b-a}{2} \end{aligned}

which maps the limits [a, b] to [-1, 1] and makes the integral

\begin{aligned} \int_{a}^{b} f(x) \,\textrm{d}x = \frac{b-a}{2} \int_{-1}^{1} f\left( \frac{b-a}{2}\xi + \frac{a + b}{2}\right) \,\textrm{d} \xi. \end{aligned}

Then by writing f\left( \frac{b-a}{2}\xi + \frac{a + b}{2}\right) = c_0 + c_1 \xi + c_2 \xi^2 + ... + c_{2n-1}\xi^{2n-1} = q(\xi)P_n(\xi) + r(\xi) and repeating the process above, it becomes clear that the sample points can be \xi_i = the ith zero of P_n(\xi) in the range [-1, +1] and that the weights are the same as before. Specifically the function would be sampled as \sum_i w_i f\left( \frac{b-a}{2}\xi_i + \frac{a+b}{2} \right) with (\xi_0, \xi_1, \xi_2) = (-\sqrt{\frac{3}{5}}, 0 , \sqrt{\frac{3}{5}}) and (w_0, w_1, w_2) = (\frac{5}{9}, \frac{8}{9}, \frac{5}{9}).

As a specific example the function f(x) = 5x^4 can be integrated between the limits (-2, 10) with the code below. It gives 100032.0 which is exactly the right answer.


import numpy as np
import scipy.special as sp

def f(x):
    return 5*x**4

roots, weights = sp.roots_legendre(3)

lim_0 = -2
lim_1 = 10
quad_result = (lim_1-lim_0)/2*np.sum(weights * f((lim_1-lim_0)/2*roots + (lim_0+lim_1)/2))
print(f"Quadrature = {quad_result:.8f}")
print(f"Analytical = {100032:.8f}")

In making this step, a method has been demonstrated to evaluate polynomial functions up to degree 2n-1 using only n sample points with any finite limits.

Other limits and Weight Functions

If an integral contains an infinite limit then the formulas above will be difficult to use. Whilst other methods to evaluate these improper integrals do exist [2], there can be complications when using them in practice; for example in [2] the midpoint Riemann sum is used to avoid evaluating \frac{1}{x} and x=0. Furthermore, functions such as e^{-x} do not have finite polynomial expansions and so can’t be integrated exactly with the outlined method. In these cases it might be possible to use the Gaussian quadrature technique using the zeros of other orthogonal polynomials.

Integrals of the form

\begin{aligned} \int_{a}^{b} w(x) f(x) \,\textrm{d} x, \end{aligned}

can be evaluated exactly using this method, if the weights and limits are as given in [3, pp.774; 4]. Once a set of suitable polynomials have been identified, their zeros and the weight coefficients can be found. One example will be given for Gauss-Laguerre quadrature, but other weight functions and limits exist.

Gauss-Laguerre Quadrature

Gauss-Laguerre quadrature can be used for weight functions w(x) = e^{-x} and limits a=0, b=\infty. In this case the sample locations are the zeros of the n^{\textrm{th}} order Laguerre polynomial; for n=3 these are (x_0, x_1, x_2) = (0.41577456, 2.29428036, 6.28994508).

To calculate the weights the integrals \int_{0}^{\infty} x^{i}e^{-x} \,\textrm{d}x have to be computed for i=0,1,2, .... This gives 1, 1, 2 for the first three, which are then plugged into a matrix similar to the one above

\begin{pmatrix} 1 & 1 & 1 \\ x_0 & x_1 & x_2 \\ x_0^2 & x_1^2 & x_2^2 \end{pmatrix} \begin{pmatrix} w_0 \\ w_1 \\ w_2 \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \\ 2 \end{pmatrix}.

Solving this for the weights gives (w_0, w_1, w_2) = (0.71109301, 0.27851773, 0.01038926). The code to achieve this is given below.


import numpy as np
import scipy.special as sp

roots, weights = sp.roots_laguerre(3)

a = np.ones((len(roots), len(roots)))
for i in range(1, len(roots)):
    a[i] = roots**i

b = np.array([1, 1, 2])

my_weights = np.linalg.solve(a, b)
print(f"Calculated weights: {my_weights}")
print(f"Weights from Scipy: {weights}")

As an example the integral

\begin{aligned} \int_{0}^{\infty} e^{-x} L_{10}^{(3)}(x) \,\textrm{d}x \end{aligned}

will be evaluated, where L_{n}^{(\alpha)}(x) is a generalised Laguerre polynomial. Since L_{10}^{(3)}(x) is of degree 10, exact evaluation of this function requires at least 6 samples since 2 \times 6 - 1 = 11 > 10. Computing with five samples gives 65.99603, whilst >=6 samples give 66.0, which is also the analytical answer from L_{10}^{(3)}(0) - L_9^{(3)}(0) = 66 [3, 22.13.12]. The code to demonstrate this is given below.


import numpy as np
import scipy.special as sp

L10 = sp.genlaguerre(10, 3)

def f(x):
    return L10(x)

def integrand(x):
    return np.exp(-x)*f(x)

roots, weights = sp.roots_laguerre(6)
print(f"Quadrature = {np.sum(weights * f(roots)):.5f}")

L9 = sp.genlaguerre(9, 3)
print(f"Analytical = {L10(0) - L9(0):.5f}")

Conclusion

In this post I introduced the general idea of quadrature rules, where n samples of a function can be used to exactly integrate a polynomial of degree n-1. I then explained how this was expanded with Gaussian quadrature, by using the properties of orthogonal polynomials so that n samples can give some exact integrals of polynomials with degree 2n-1, provided the sample locations are chosen as the zeros of certain orthogonal polynomials. I have shown how that is extended to very general integrals with finite limits and how some variations of the technique, such as Gauss-Laguerre quadrature, can evaluate some integrals with infinite limits.

References

[1] – BCcampus Open Publishing, Dividing Polynomials, Accessed 07/2022.
[2] – W. H. Press, “Numerical Recipes: The Art of Scientific Computing”, https://numerical.recipes
[3] – M. Abbramowitz & I. Stegun, “Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables”, https://archive.org/details/AandS-mono600
[4] – Digital Library of Mathematical Functions, https://dlmf.nist.gov/18.3 accessed 07/2022.

Published by ComPhys

Leave a comment