Gaussian Quadrature

The numerical integration methods described so far are based on a rather simple choice of evaluation points for the function f(x). They are particularly suited for regularly tabulated data, such as one might measure in a laboratory, or obtain from computer software designed to produce tables. If one has the freedom to choose the points at which to evaluate f(x), a careful choice can lead to much more accuracy in evaluating the integral in question. We shall see that this method, called Gaussian or Gauss-Legendre integration, has one significant further advantage in many situations. In the evaluation of an integral on the interval $\alpha$ to $\beta$, it is not necessary to evaluate f(x) at the endpoints, ie. at $\alpha$ or $\beta$, of the interval. This will prove valuable when evaluating various improper integrals, such as those with infinite limits.

Figure 7.4: Comparing trapezoid rule integration and Gaussian integration.

We begin with a simple example illustrated in Figure 7.4.

The simplest form of Gaussian Integration is based on the use of an optimally chosen polynomial to approximate the integrand f(t) over the interval [-1,+1]. The details of the determination of this polynomial, meaning determination of the coefficients of t in this polynomial, are beyond the scope of this presentation. The simplest form uses a uniform weighting over the interval, and the particular points at which to evaluate f(t) are the roots of a particular class of polynomials, the Legendre polynomials, over the interval. It can be shown that the best estimate of the integral is then:

$$\int^1_{-1} f(t) dt = \sum_{i=1}^n w_i f(t_i)$$

where t_i is a designated evaluation point, and w_i is the weight of that point in the sum. If the number of points at which the function f(t) is evaluated is n, the resulting value of the integral is of the same accuracy as a simple polynomial method (such as Simpson's Rule) of about twice the degree (ie. of degree 2n). Thus the carefully designed choice of function evaluation points in the Gauss-Legendre form results in the same accuracy for about half the number of function evaluations, and thus at about half the computing effort.

Gaussian quadrature formulae are evaluating using abscissae and weights from a table like that included here. The choice of value of n is not always clear, and experimentation is useful to see the influence of choosing a different number of points. When choosing to use n points, we call the method an ``n-point Gaussian'' method.

Gauss-Legendre Abscissae and Weights
n	Values of t	Weights	Degree
2	$\pm\sqrt{\frac{1}{3}}$	1.0	3
3	0.0	0.88888889
	$\pm$ 0.77459667	0.55555555	5
4	$\pm$ 0.33998104	0.65214515	7
	$\pm$ 0.86113631	0.34785485
5	0.0	0.56888889	9
	$\pm$ 0.53846931	0.47862867
	$\pm$ 0.90617985	0.23692689
6	$\pm$ 0.23861918	0.46791393	11
	$\pm$ 0.66120939	0.36076157
	$\pm$ 0.93246951	0.17132449
7	0.0	0.41795918	13
	$\pm$ 0.40584515	0.38183005
	$\pm$ 0.74153119	0.27970539
	$\pm$ 0.94910791	0.12948497
8	$\pm$ 0.18343464	0.36268378	15
	$\pm$ 0.52553241	0.31370665
	$\pm$ 0.79666648	0.22238103
	$\pm$ 0.96028986	0.10122854
10	$\pm$ 0.14887434	0.29552422	19
	$\pm$ 0.43339539	0.26926672
	$\pm$ 0.67940957	0.21908636
	$\pm$ 0.86506337	0.14945135
	$\pm$ 0.97390653	0.06667134

The Gauss-Legendre integration formula given here evaluates an estimate of the required integral on the interval for t of [-1,+1]. In most cases we will want to evaluate the integral on a more general interval, say $[\alpha,\beta]$. We will use the variable x on this more general interval, and linearly map the $[\alpha,\beta]$ interval for x onto the [-1,+1] interval for t using the linear tranformation:

$$x = c + mt \ \ \ {\rm where} \ \ \ c = \frac{1}{2}(b + a) \ \ \ {\rm and} \ \ \ m = \frac{1}{2}(b - a)$$

It is easily verified that substituting t = -1 gives x = a and t = 1 gives x = b. We can now write the integral as:

$$I = \int_{a}^{b}\ f(x)\ dx = m \int_{-1}^{+1}\ f(c + mt)\ dt$$

The factor of m in the second integral arises from the change of the variable of integration from x to t, which introduces the factor $\frac{dx}{dt}$. Finally, we can write the Gauss-Legendre estimate of the integral as:

$$I = \int_{a}^{b}\ f(x)\ dx = m \sum_{i=1}^n w_i f(c + m t_i)$$

Consider the evaluation of the integral:

$$I = \int_0^{\pi/2} \sin x \ dx$$

whose value is 1, as we can obtain by explicit integration. Applying the 2-point Gaussian method, and noting that both c and m are $\pi/4$, the table allows us to calculate an approximate value for the integral. The result is 0.998473, which is pretty close to the exact value of one. The calculation is simply:

$$I \approx \frac{\pi}{4} \left ( 1.0 \times \sin \left ( \le... ...+ \frac{1}{\sqrt{3}} \right ) \frac{\pi}{4} \right ) \right )$$

While this example is quite simple, the following table of values obtained for n ranging from 2 to 10 indicates how accurate the estimate of the integral is for only a few function evaluations. The table includes a column of values obtained from Simpson's $\frac{1}{3}$ rule for the same number of function evaluations. The Gauss-Legendre result is correct to almost twice the number of digits as compared to the Simpson's rule result for the same number of function evaluations.

N	Gauss-Legendre	Simpson's 1/3
2	0.9984726135	1.0022798775
4	0.9999999770	1.0001345845
6	0.9999999904	1.0000263122
8	1.0000000001	1.0000082955
10	0.9999999902	1.0000033922