1.16: Gaussian Quadrature - Derivation

In order to understand why Gaussian quadrature works so well, we first need to understand some properties of polynomials in general, and of Legendre polynomials in particular. We also need to remind ourselves of the use of Lagrange polynomials for approximating an arbitrary function.

First, a statement concerning polynomials in general: Let \(P\) be a polynomial of degree \(n\), and let \(S\) be a polynomial of degree less than \(2n\). Then, if we divide \(S\) by \(P\), we obtain a quotient \(Q\) and a remainder \(R\), each of which is a polynomial of degree less than \(n\).

That is to say: \[\frac{S}{P} = Q + \frac{R}{P}. \label{1.16.1} \tag{1.16.1}\]

What this means is best understood by looking at an example, with \(n = 3\). For example,

let \[P = 5x^3 - 2x^2 + 3x + 7 \label{1.16.2} \tag{1.16.2}\]

and \[S = 9x^5 + 4x^4 - 5x^3 + 6x^2 + 2x - 3. \label{1.16.3} \tag{1.16.3}\]

If we carry out the division \(S ÷ P\) by the ordinary process of long division, we obtain

\[\frac{9x^5 + 4x^4 - 5x^3 + 6x^2 + 2x - 3}{5x^3 - 2x^2 + 3x + 7} = 1.8x^2 + 1.52x - 1.472 - \frac{14.104x^2 + 4.224x - 7.304}{5x^3 - 2x^2 + 3x + 7}. \label{1.16.4} \tag{1.16.4}\]

For example, if \(x = 3\), this becomes

\[\frac{2433}{133} = 19.288 - \frac{132.304}{133}. \]

The theorem given by Equation \(\ref{1.16.1}\) is true for any polynomial \(P\) of degree \(l\). In particular, it is true if \(P\) is the Legendre polynomial of degree \(l\).

__________________________________

Next an important property of the Legendre polynomials, namely, if \(P_n\) and \(P_m\) are Legendre polynomials of degree \(n\) and \(m\) respectively, then

\[\int_{-1}^1 P_n P_m dx = 0 \quad \text{unless } m = n. \label{1.16.5} \tag{1.16.5}\]

This property is called the orthogonal property of the Legendre polynomials.

I give here a proof. Although it is straightforward, it may look formidable at first, so, on first reading, you might want to skip the proof and go on the next part (after the next short horizontal dividing line).

From the symmetry of the Legendre polynomials (see figure \(\text{I.7}\)), the following are obvious:

\[\int_{-1}^1 P_n P_m dx \neq 0 \quad \text{if } m=n\]

and \[\int_{-1}^1 P_n P_m = 0 \quad \text{if one (but not both) of } m \text{ or } n \text{ is odd}.\]

In fact we can go further, and, as we shall show,

\[\int_{-1}^1 P_n P_m dx = 0 \quad \text{unless } m = n , \text{ whether } m \text{ and } n \text{ are even or odd}.\]

Thus \(P_m\) satisfies the differential Equation (see Equation 1.14.7)

\[(1-x^2) \frac{d^2P_m}{dx^2} - 2x \frac{dP_m}{dx} + m(m+1) P_m = 0, \label{1.16.6} \tag{1.16.6}\]

which can also be written

\[\frac{d}{dx} \left[ (1-x^2) \frac{dP_m}{dx} \right] + m(m+1) P_m = 0. \label{1.16.7} \tag{1.16.7}\]

Multiply by \(P_n\):

\[P_n \frac{d}{dx} \left[ (1-x^2) \frac{dP_m}{dx} \right] + m(m+1)P_m P_n = 0, \label{1.16.8} \tag{1.16.8}\]

which can also be written

\[\frac{d}{dx} \left[ (1-x^2) P_n \frac{dP_m}{dx} \right] - (1-x^2) \frac{dP_n}{dx} \frac{dP_m}{dx} + m(m+1) P_m P_n = 0. \label{1.16.9} \tag{1.16.9}\]

In a similar manner, we have

\[\frac{d}{dx} \left[ (1 - x^2) P_m \frac{dP_n}{dx} \right] - (1-x^2) \frac{dP_n}{dx} \frac{dP_m}{dx} + n(n+1) P_m P_n = 0. \label{1.16.10} \tag{1.16.10}\]

Subtract one from the other:

\[\frac{d}{dx} \left[ (1-x^2) \left( P_n \frac{dP_m}{dx} - P_m \frac{dP_n}{dx} \right) \right] + [m(m+1) - n(n+1)]P_m P_n = 0. \label{1.16.11} \tag{1.16.11}\]

Integrate from \(−1\) to \(+1\):

\[ \left[ (1-x^2) \left( P_n \frac{dP_m}{dx} - P_m \frac{dP_n}{dx} \right) \right]_{-1}^1 = [n(n+1) - m(m+1)] \int_{-1}^1 P_m P_n dx. \label{1.16.12} \tag{1.16.12}\]

The left hand side is zero because \(1 − x^2\) is zero at both limits.

Therefore, unless \(m = n\),

\[\int_{-1}^1 P_m P_n dx = 0 . \quad \quad \text{Q.E.D.} \label{1.16.13} \tag{1.16.13}\]

___________________________________

I now assert that, if \(P_l\) is the Legendre polynomial of degree \(l\), and if \(Q\) is any polynomial of degree less than \(l\), then

\[\int_{-1}^1 P_l Q dx = 0. \label{1.16.14} \tag{1.16.14}\]

I shall first prove this, and then give an example, to see what it means.

To start the proof, we recall the recursion relation (see Equation 1.14.4 – though here I am substituting \(l − 1\) for \(l\)) for the Legendre polynomials:

\[lP_l = (2l-1) xP_{l-1} - (l-1) P_{l-2} . \label{1.16.15} \tag{1.16.15}\]

The proof will be by induction.

Let \(Q\) be any polynomial of degree less than l. Multiply the above relation by \(Qdx\) and integrate from \(−1\) to \(+1\):

\[l \int_{-1}^1 P_l Q dx = (2l-1) \int_{-1}^1 x P_{l-1} Q dx - (l-1) \int_{-1}^1 P_{l-2} Q dx. \label{1.16.16} \tag{1.16.16}\]

If the right hand side is zero, then the left hand side is also zero.

A correspondent has suggested to me a much simpler proof. He points out that you could in principle expand \(Q\) in Equation \(\ref{1.16.14}\) as a sum of Legendre polynomials for which the highest degree is \(l-1\). Then, by virtue of Equation \(\ref{1.16.13}\), every term is zero.

For example, let \(l = 4\), so that

\[P_{l-2} = P_2 = \frac{1}{2} (3x^2 - 1) \label{1.16.17} \tag{1.16.17}\]

and \[xP_{l-1} = xP_3 = \frac{1}{2} (5x^4 - 3x^2), \label{1.16.18} \tag{1.16.18}\]

and let \[Q = 2(a_3 x^3 + a_2 x^2 + a_1 x + a_0 ) . \label{1.16.19} \tag{1.16.19}\]

It is then straightforward (and only slightly tedious) to show that

\[\int_{-1}^1 P_{l-2} Q dx = \left( \frac{6}{5} - \frac{2}{3} \right) a_2 \label{1.16.20} \tag{1.16.20}\]

and that \[\int_{-1}^1 xP_{l-1} Q dx = \left( \frac{10}{7} - \frac{6}{5} \right) a_2 . \label{1.16.21} \tag{1.16.21}\]

But \[7 \left( \frac{10}{7} - \frac{6}{5} \right) a_2 - 3 \left( \frac{6}{5} - \frac{2}{3} \right) a_2 = 0, \label{1.16.22} \tag{1.16.22}\]

and therefore \[\int_{-1}^1 P_4 Q dx = 0 . \label{1.16.23} \tag{1.16.23}\]

We have shown that \[l \int_{-1}^1 P_l Q dx = (2l - 1) \int_{-1}^1 x P_{l-1} Q dx - (l - 1) \int_{-1}^1 P_{l-2} Q dx = 0 \label{1.16.24} \tag{1.16.24}\]

for \(l = 4\), and therefore it is true for all positive integral \(l\).

You can use this property for a parlour trick. For example, you can say: “Think of any polynomial. Don’t tell me what it is – just tell me its degree. Then multiply it by (here give a Legendre polynomial of degree more than this). Now integrate it from \(−1\) to \(+1\). The answer is zero, right?” (Applause.)

Thus: Think of any polynomial. \(3x^2 - 5x + 7\). Now multiply it by \(5x^3 - 3x\). OK, that’s \(15x^5 - 25x^4 - 2x^3 + 15x^2 - 21x\). Now integrate it from \(−1\) to \(+1\). The answer is zero.

__________________________________________

Now, let \(S\) be any polynomial of degree less than \(2l\). Let us divide it by the Legendre polynomial of degree \(l\), \(P_l\), to obtain the quotient \(Q\) and a remainder \(R\), both of degree less than \(l\). Then I assert that

\[\int_{-1}^1 Sdx = \int_{-1}^1 Rdx. \label{1.16.25} \tag{1.16.25}\]

This follows trivially from Equations \(\ref{1.16.1}\) and \(\ref{1.16.14}\). Thus

\[\int_{-1}^1 S dx = \int_{-1}^1 (QP_l + R) dx = \int_{-1}^1 Rdx . \label{1.16.26} \tag{1.16.26}\]

Example: Let \(S = 6x^5 - 12x^4 + 4x^3 + 7x^2 - 5x + 7\). The integral of this from \(−1\) to \(+1\) is \(13.86\). If we divide \(S\) by \(\frac{1}{2} (5x^3 - 3x)\), we obtain a quotient of \(2.4x^2 - 4.8x + 3.04\) and a remainder of \(-0.2x^2 - 0.44x + 7\). The integral of the latter from \(−1\) to \(+1\) is also \(13.86\).

______________________________________

I have just described some properties of Legendre polynomials. Before getting on to the rationale behind Gaussian quadrature, let us remind ourselves from Section 1.11 about Lagrange polynomials. We recall from that section that, if we have a set of n points, the following function:

\[y = \sum_{i=1}^n y_i L_i (x) \label{1.16.27} \tag{1.16.27}\]

(in which the \(n\) functions \(L_i(x )\), \(i = 1,n\), are Lagrange polynomials of degree \(n-1)\) is the polynomial of degree \(n-1\) that passes exactly through the \(n\) points. Also, if we have some function \(f(x)\) which we evaluate at \(n\) points, then the polynomial

\[ y = \sum_{i=1}^n f(x_i) L_i (x) \label{1.16.28} \tag{1.16.28}\]

is a jolly good approximation to \(f(x)\) and indeed may be used to interpolate between nontabulated points, even if the function is tabulated at irregular intervals. In particular, if \(f(x)\) is a polynomial of degree \(n − 1\), then the expression \(\ref{1.16.28}\) is an exact representation of \(f(x)\).

________________________________

We are now ready to start talking about quadrature. We wish to approximate \(\int_{-1}^1 f(x) dx\) by an \(n\)-term finite series

\[\int_{-1}^1 f(x) dx \approx \sum_{i=1}^n c_i f(x_i), \label{1.16.29} \tag{1.16.29}\]

where \(-1 < x_i < 1\). To this end, we can approximate \(f(x)\) by the right hand side of Equation \(\ref{1.16.28}\), so that

\[\int_{-1}^1 f(x) dx \approx \int_{-1}^1 \sum_{i=1}^n f(x_i) L_i (x) dx =\sum_{i=1}^nf(x_i)\int_{-1}^1L_i(x)dx. \label{1.16.30} \tag{1.16.30}\]

Recall that the Lagrange polynomials in this expression are of degree \(n − 1\).

The required coefficients for Equation \(\ref{1.16.29}\) are therefore

\[c_i = \int_{-1}^1 L_i (x) dx. \label{1.16.31} \tag{1.16.31}\]

Note that at this stage the values of the \(x_i\) have not yet been chosen; they are merely restricted to the interval [−1 , 1].

__________________________________

Now let’s consider \(\int_{-1}^1 S(x) dx\), where \(S\) is a polynomial of degree less than \(2n\), such as, for example, the polynomial of Equation \(\ref{1.16.3}\). We can write

\[\int_{-1}^1 S(x) dx = \int_{-1}^1 \sum_{i=1}^n S(x_i) L_i (x) dx = \int_{-1}^1 \sum_{i=1}^n L_i (x) [Q(x_i) P(x_i) + R(x_i)] dx. \label{1.16.32} \tag{1.16.32}\]

Here, as before, \(P\) is a polynomial of degree \(n\), and \(Q\) and \(R\) are of degree less than \(n\).

If we now choose the \(x_i\) to be the roots of the Legendre polynomials, then

\[\int_{-1}^1 S(x) dx = \int_{-1}^1 \sum_{i=1}^n L_i (x) R(x_i) dx. \label{1.16.33} \tag{1.16.33}\]

Note that the integrand on the right hand side of Equation \(\ref{1.16.33}\) is an exact representation of \(R(x)\). But we have already shown (Equation \(\ref{1.16.26}\)) that \(\int_{-1}^1 S(x) dx = \int_{-1}^1 R(x) dx\), and therefore

\[\int_{-1}^1 S(x) dx = \int_{-1}^1 R(x) dx = \sum_{i=1}^n c_i R(x_i) = \sum_{i=1}^n c_i S(x_i) . \label{1.16.34} \tag{1.16.34}\]

It follows that the Gaussian quadrature method, if we choose the roots of the Legendre polynomials for the \(n\) abscissas, will yield exact results for any polynomial of degree less than \(2n\), and will yield a good approximation to the integral if \(S(x)\) is a polynomial representation of a general function \(f(x)\) obtained by fitting a polynomial to several points on the function.