1.15: Gaussian Quadrature - the Algorithm

Gaussian quadrature is an alternative method of numerical integration which is often much faster and more spectacular than Simpson’s rule. Gaussian quadrature allows you to carry out the integration

$$\int_{-1}^1 f(x) dx. \label{1.15.1} \tag{1.15.1}$$

But what happens if your limits of integration are not $\pm 1$? What if you want to integrate

$$\int_a^b F(t) dt? \label{1.15.2} \tag{1.15.2}$$

That is no problem at all – you just make a change of variable. Thus, let

$$x = \frac{2t-a-b}{b-a} , \quad t= \frac{1}{2} [(b-a)x+a+b], \label{1.15.3} \tag{1.15.3}$$

and the new limits are then $x = \pm 1$.

At the risk of being pedagogically unsound I’ll describe first, without any theoretical development, just what you do, with an example – as long as you promise to look at the derivation afterwards, in Section 1.16.

For our example, let’s try to evaluate

$$I = \int_0^{\pi/2} \sin \theta d \theta. \label{1.15.4} \tag{1.15.4}$$

Let us make the change of variable given by Equation $\ref{1.15.3}$ (with $t = \theta$, $a = 0, \ b = \pi/2$), and we now have to evaluate

$$I = \int_{-1}^1 \frac{\pi}{4} \sin \frac{\pi}{4} (x+1) dx . \label{1.15.5} \tag{1.15.5}$$

For a 5-point Gaussian quadrature, you evaluate the integrand at five values of $x$, where these five values of $x$ are the solutions of $P_5(x) = 0$ given in Section 1.14, $P_5$ being the Legendre polynomial. That is, we evaluate the integrand at $x = \pm 0.906 \ 469 \ 514 \ 203, \ \pm 0.538 \ 469 \ 310 \ 106$ and $0$.

I now assert, without derivation (until later), that

$$I = \sum_{i=1}^5 c_{5,i} \ f(x_{5,i}), \label{1.15.6} \tag{1.15.6}$$

where the coefficients $c_{l , i}$ (all positive) are listed with the roots of the Legendre polynomials in Section 1.14.

Let’s try it.

$$\begin{array}{c c c} x_{5,i} & f(x_{5,i}) & c_{5,i} \\ \\ +0.906 \ 179 \ 845 \ 939 & 0.783 \ 266 \ 908 \ 39 & 0.236 \ 926 \ 885 \ 06 \\ +0.538 \ 469 \ 310 \ 106 & 0.734 \ 361 \ 739 \ 69 & 0.478 \ 628 \ 670 \ 50 \\ 0.000 \ 000 \ 000 \ 000 & 0.555 \ 360 \ 367 \ 27 & 0.568 \ 888 \ 888 \ 89 \\ -0.538 \ 469 \ 310 \ 006 & 0.278 \ 501 \ 544 \ 60 & 0.478 \ 628 \ 670 \ 50 \\ -0.906 \ 179 \ 845 \ 939 & 0.057 \ 820 \ 630 \ 35 & 0.236 \ 926 \ 885 \ 06 \\ \end{array}$$

and the expression $\ref{1.15.6}$ comes to $1.000 \ 000 \ 000 \ 04$, and might presumably have come even closer to 1 had we given $x_{l,i}$, and $c_{l,i}$, to more significant figures.

You should now write a computer program for Gaussian quadrature – you will have to store the $x_{l,i}$ and $c_{l,i}$ of course. You have presumably already written a program for Simpson’s rule.

In a text on integration, the author invited the reader to evaluate the following integrals by Gaussian quadrature:

$$\begin{array}{c c c c c} (a) & \int_1^{1.5} x^2 \ln x dx & & (e) & \int_0^{\pi/4} e^{3x} \sin 2x dx \\ (b) & \int_0^1 x^2 x^{-x} dx && (f) & \int_1^{1.6} \frac{2x}{x^2 - 4} dx \\ (c) & \int_0^{0.35} \frac{2}{x^2-4} dx & & (g) & \int_3^{3.5} \frac{x}{\sqrt{x^2 - 4}}dx \\ (d) & \int_0^{\pi/4} x^2 \sin x dx && (h) & \int_0^{\pi/4} \cos^2 x dx \\ \end{array}$$

All of these can be integrated analytically, so I am going to invite the reader to evaluate them first analytically, and then numerically by Simpson’s rule and again by Gaussian quadrature, and to see at how many points the integrand has to be evaluated by each method to achieve nine or ten figure precision. I tried, and the results are as follows. The first column is the answer, the second column is the number of points required by Simpson’s rule, and the third column is the number of points required by Gaussian quadrature.

$$\begin{array}{c c c c c} (a) & 0.192 \ 259 \ 358 && 33 & 4 \\ (b) & 0.160 \ 602 \ 794 && 99 & 5 \\ (c) & −0.176 \ 820 \ 020 && 19 & 4 \\ (d) & 0.088 \ 755 \ 284 \ 4 && 111 & 5 \\ (e) & 2.588 \ 628 \ 633 && 453 & 7 \\ (f) & -0.733 \ 969 \ 175 && 143 & 8 \\ (g) & 0.636 \ 213 \ 346 && 31 & 5 \\ (h) & 0.642 \ 699 \ 082 && 59 & 5 \\ \end{array}$$

Let us now have a look at four of the integrals that we met in Section 1.2.

1. $\int_0^1 \frac{x^4 dx}{\sqrt{2(1+x^2)}}$. This was straightforward. It has an analytic solution of $\frac{\sqrt{18} \ln (1 + \sqrt{2}) -2 }{16} = 0.108 \ 709 \ 465$. I needed to evaluate the integral at 89 points in order to get this answer to nine significant figures using Simpson’s rule. To use Gaussian quadrature, we note that integrand contains only even powers of $x$ and so it is symmetric about $x = 0$, and therefore the integral is equal to $\frac{1}{2} \int_{-1}^1 \frac{x^4 dx}{\sqrt{2(1+x^2)}}$, which makes it immediately convenient for Gaussian quadrature! I give below the answers I obtained for 3- to 7-point Gaussian quadrature.

$$\begin{array}{c c c} & 3 & 0.108 \ 667 \ 036 \\ & 4 & 0.108 \ 711 \ 215 \\ & 5 & 0.108 \ 709 \ 441 \\ & 6 & 0.108 \ 709 \ 463 \\ & 7 & 0.108 \ 709 \ 465 \\ \text{Correct answer} & & 0.108 \ 709 \ 465 \\ \end{array}$$

2. $\int_0^2 \frac{y^2 dy}{\sqrt{2-y}}$. This had the difficulty that the integrand is infinite at the upper limit. We got round this by means of the substitution $y = 2 \sin^2 \theta$, and the integral becomes $\sqrt{128} \int_0^{\pi/2} \sin^5 \theta d \theta$. This has an analytic solution of $\sqrt{8192}/15 = .6 033977866$. I needed 59 points to get this answer to ten significant figures using Simpson’s rule. To use Gaussian quadrature we can let $y = 1 + x$, so that the integral becomes , $\int_{-1}^1 \frac{(1+x)^2 dy}{\sqrt{1-x}},$ which seems to be immediately suitable for Gaussian quadrature. Before we proceed, we recall that the integrand becomes infinite at the upper limit, and it still does so after our change of variable. We note, however, that with Gaussian quadrature, we do not evaluate the integrand at the upper limit, so that this would appear to be a great advantage of the method over Simpson’s method. Alas! – this turns out not to be the case. If, for example, we use a 17-point quadrature, the largest value of $x$ for which we evaluate the integrand is equal to the largest solution of $P_{17} (x) = 0$, which is $0.9906$. We just cannot ignore the fact that the integrand shoots up to infinity beyond this, so we have left behind a large part of the integral. Indeed, with a 17-point Gaussian quadrature, I obtained an answer of $5.75$, which is a long way from the correct answer of $6.03$.

Therefore we have to make a change of variable, as we did for Simpson’s method, so that the upper limit is finite. We chose $y = 2 \sin^2 \theta$ which changed the integral to $\sqrt{128} \int_0^{\pi/2} \sin^5 \theta d \theta$. To make this suitable for Gaussian quadrature, we must now make the further substitution (see Equation $\ref{1.15.3}$) $x = 4\theta/\pi-1$, $\theta = \frac{\pi}{4} (x+1)$. If we wish to impress, we can make the two substitutions in one step, thus: Let $y = 2 \sin^2 \frac{\pi}{4}(1+x)$, $x = \frac{4}{\pi} \sin^{-1} \sqrt{\frac{y}{2}} - 1$. The integral becomes $\sqrt{8} \pi \int_{-1}^1 \sin^5 \frac{\pi}{4}(1+x)dx$, and there are no further difficulties. With a 9-point integration, I obtained the answer, correct to ten significant figures, $6.033 \ 977 \ 866$. Simpson’s rule required 59 points.

3. $\int_0^{\pi/2} \sqrt{\sec \theta} d \theta$. This integral occurs in the theory of a simple pendulum swinging through $90^\circ$. As far as I can tell it has no simple analytical solution unless we have recourse to unfamiliar elliptic integrals, which we would have to evaluate numerically in any case. The integral has the difficulty that the integrand is infinite at the upper limit. We get round this by means of a substitution. Thus let $\sin \phi = \sqrt{2} \sin \frac{1}{2} \theta$. (Did you not think of this?) The integral becomes $\sqrt{2} \int_0^{\pi/2} \frac{d\phi}{\sqrt{1 - \frac{1}{2} \sin^2 \phi}}.$ I needed 13 points by Simpson’s rule to get the answer to ten significant figures, $2.622 \ 057 \ 554$. In order to make the limits $\pm 1$, suitable for Gaussian quadrature, we can make the second substitution (as in example 2), $\phi = \frac{\pi}{4}(x+1)$. If we wish truly to impress our friends, we can make the two substitutions in one step, thus: Let $\sin \frac{\pi}{4} (1+x) = \sqrt{2} \sin \frac{1}{2} \theta$. (No one will ever guess how we thought of that!) The integral becomes $\frac{\pi}{2} \int_{-1}^1 \frac{dx}{\sqrt{2 - \sin^2 \frac{\pi}{4} (x+1)}},$ which is now ready for Gaussian quadrature. I obtained the answer $2.622 \ 057 \ 554$ in a 10-point Gaussian quadrature, which is only a little faster than the 13 points required by Simpson’s rule.

4. $\int_0^\infty \frac{dy}{y^5\left(e^{1/y} - 1 \right)}.$ This integral occurs in the theory of blackbody radiation. It has the difficulty of an infinite upper limit. We get round this by means of a substitution. Thus let $y = \tan \theta$. The integral becomes $\int_0^{\pi/2} \frac{c^3(c^2+1)}{e^c - 1} d\theta$, where $c = \cot \theta$. It has an analytic solution of $\pi^4 /15 = 6.493 \ 939 \ 402$. I needed 261 points by Simpson’s rule to get the answer to ten significant figures. To prepare it for Gaussian quadrature, we can let $\theta = \frac{\pi}{4}(x+1)$, as we did in example 2, so that the integral becomes ,where $\frac{\pi}{4} \int_{-1}^1 \frac{c^3 (c^2+1)}{e^c -1}dx,$ where $c = \cot \frac{\pi}{4} (x+1)$. Using 16- point Gaussian quadrature, I got 6.48. Thus we would need to extend our table of constants for the Gaussian method to much higher order in order to use the method successfully. Doubtless the Gaussian method would then be faster than the Simpson method – but we do not need an extensive (and difficult-to-calculate) set of constants for the latter. A further small point: You may have noticed that it is not immediately obvious that the integrand is zero at the end points, and that some work is needed to prove it. But with the Gaussian method you don’t evaluate the integrand at the end points, so that is one less thing to worry about!

Thus we have found that in most cases the Gaussian method is far faster than the Simpson method. In some cases it is only marginally faster. In yet others it probably would be faster than Simpson’s rule, but higher-order constants are needed to apply it. Whether we use Simpson’s rule or Gaussian quadrature, we have to carry out the integration with successively higher orders until going to higher orders results in no further change to the number of significant figures desired.