Numerical Integration with Legendre Polynomials
A surprising usage of the Legendre polynomials is probably the numerical integration of a function over a bounded interval. Suppose that a function \(f(x)\) defined over \( [-1,1] \) is first approximated using Legendre polynomials at \( n \) points. Particularly, these \( n \) points \( x_i , i = 1,2,\ldots,n \) are chosen to be the \( n \) zeros of the Legendre polynomial of order \( n \): \( P_n(x) \), so that \( P(x_i) = 0 \). Then we may mimic \( f(x) \) by the following polynomial of degree \( n -1 \):
\begin{align}
\tilde{f}(x) = \sum_{i=1}^n \frac{P_n(x)}{(x-x_i)P’_n(x_i)} f(x_i) \tag{1}
\end{align}
This new polynomial will then coincide with the original \(f(x)\) at each assigned point \( x_i \). It is because
\begin{align}
\lim_{x\to x_j} \frac{P_n(x)}{(x-x_i)P’_n(x_i)} =
\begin{cases}
\frac{P_n'(x_i)}{P’_n(x_i)} &= 1 & \text{if \( i = j \)} \\
\frac{P_n(x_j)}{(x_j-x_i)P’_n(x_i)} &= 0 & \text{if \( i \neq j \)}
\end{cases} \tag{2}
\end{align}
where the first case is by L’Hospital’s Rule, so the sum in (1) when \( x = x_i \) only has the \( f(x_i) \) term remaining. And then we will just approximate the integral of \( f(x) \) with
\begin{align}
\int_{-1}^{1} f(x) dx &\approx \int_{-1}^{1} \tilde{f}(x) dx \\
&= \sum_{i=1}^n (\frac{1}{P’_n(x_i)} \int_{-1}^{1} \frac{P_n(x)}{(x-x_i)} dx) f(x_i) = \sum_{i=1}^n w_i f(x_i) \tag{3}
\end{align}
\( w_i \) denotes the weightings in the front. This is then known as the Gauss-Legendre Quadrature rule. There are two key results that we are going to show. (Reference: this StackExchange post) The first one is that such a quadrature can exactly integrate a polynomial \( p(x) \) of degree \( 2n -1 \) or lower. Let’s define the remainder \( R(x) = p(x) -\sum_{i=1}^n \frac{P_n(x)}{(x-x_i)P’_n(x_i)} p(x_i) \). This quantity, by the same logic for (2), vanishes at every \( x_i \), which are the zeros of \( P_n(x) \). Therefore, it must be some (polynomial) multiple of \( P_n(x) \), i.e. \( R(x) = P_n(x) Q(x) \), where \( Q(x) \) is a polynomial of degree \( \leq n -1 \). So
\begin{align}
\int_{-1}^{1} p(x) dx &= \sum_{i=1}^n (\frac{1}{P’_n(x_i)} \int_{-1}^{1} \frac{P_n(x)}{(x-x_i)} dx) p(x_i) + \int_{-1}^{1} R(x) dx \\
&= \sum_{i=1}^n w_i p(x_i) + \int_{-1}^{1} P_n(x) Q(x) dx \tag{4}
\end{align}
Notice that the second integral must be zero due to orthogonality: the polynomial \( Q(x) \) can be written exactly as a Legendre series of order \( n -1 \) by the previous tutorial, and they are all orthogonal to \( P_n(x) \). So the exact value of \( \int_{-1}^{1} p(x) dx \) will be equal to the first integral, which is just the Gauss-Legendre Quadrature (3) for \( p(x) \).
Next, we derive an alternative expression of the weightings in (3):
\begin{align}
w_i = \frac{2}{(1-x_i^2)[P_n'(x_i)]^2} \tag{5}
\end{align}
To this end, we consider the integral of
\begin{align}
\int_{-1}^{1} P_n'(x)\frac{P_n(x)}{x-x_i} dx \tag{6}
\end{align}
The integrand is a polynomial of degree \( 2n-2\) (both \( P_n'(x) \) and \( P_n(x)/(x-x_i) \) are of degree \( n-1 \)), so Gauss-Legendre quadrature will yield the exact value as shown above:
\begin{align}
\int_{-1}^{1} P_n'(x)\frac{P_n(x)}{x-x_i}dx = \sum_{j=1}^n w_j P_n'(x_j)\frac{P_n(x_j)}{x_j-x_i} \tag{7}
\end{align}
By the same essence of (2), it is equal to
\begin{align}
\sum_{j=1}^n w_j P_n'(x_j)\frac{P_n(x_j)}{x_j-x_i} = w_i[P_n'(x_i)]^2 \tag{8}
\end{align}
The second way to compute (6) is via the usual integration by parts:
\begin{align}
\int_{-1}^{1} P_n'(x)\frac{P_n(x)}{x-x_i}dx &= [P_n(x)\frac{P_n(x)}{x-x_i}]_{-1}^{1} -\int_{-1}^{1} P_n(x) (\frac{P_n(x)}{x-x_i})’ dx \\
&= [\frac{P_n^2(x)}{x-x_i}]_{-1}^{1} -0 \\
&= \frac{1}{1-x_i} -\frac{1}{(-1)+x_i} = \frac{2}{1-x_i^2} \tag{9}
\end{align}
where the second integral is zero again due to orthogonality, and the Legendre polynomials follow the normalization \( P_n^2(\pm 1) = 1 \). Comparing (8) and (9) immediately leads to (5).
Example
Numerically integrate \( \cos(\pi x/2) \) over \([-1,1]\) with \( 7 \)-points Gauss-Legendre quadrature and compare with the analytical value. (A table of reference weighting and zero values can be found here)
According to (3) and the table provided (\(n=7\), we have the numerical answer of
\begin{align}
\sum_{i=1}^7 w_i f(x_i) &\approx 0.41796 \cos(\pi\frac{0}{2}) \\
&\quad + 0.38183 \cos(\pi\frac{0.40584}{2}) + 0.38183 \cos(\pi\frac{-0.40584}{2}) \\
&\quad + 0.27971 \cos(\pi\frac{0.74153}{2}) + 0.27971 \cos(\pi\frac{-0.74153}{2}) \\
&\quad + 0.12948 \cos(\pi\frac{0.94911}{2}) + 0.12948 \cos(\pi\frac{-0.94911}{2}) \\
&\approx 1.27325 \tag{10}
\end{align}
Notice that by the symmetry of the cosine function and the Legendre polynomial, the paired terms are identical so we can save some calculations. The exact value of the integral \( \int_{-1}^1 \cos(\pi x/2) dx \) is \( 4/\pi \approx 1.27324 \), so our approximation works quite well in this case.
Exercise
Use the Gauss-Legendre quadrature of \( n=6 \) to integrate \( e^x \) over \( [-1/2,1/2] \) numerically. Since the integration interval is not \( [-1,1] \) as before, we will need to do a simple change of variable.
Answer
To broadcast the given interval into \( [-1,1] \), we simply let \( t = 2x \) and then the integral becomes
\begin{align}
\int_{-\frac{1}{2}}^{\frac{1}{2}} e^x dx = \int_{-1}^{1} e^\frac{t}{2} d(\frac{t}{2}) = \frac{1}{2} \int_{-1}^{1} e^\frac{t}{2} dt
\end{align}
and we numerically integrate the integral in \(t\) to the right. By (3) and the table in the link provided above, we have
\begin{align}
\sum_{i=1}^6 w_i f(x_i) &\approx 0.36076 e^{0.66121/2} + 0.36076 e^{-0.66121/2} \\
&\quad + 0.46791 e^{0.23862/2} + 0.46791 e^{-0.23862/2} \\
&\quad + 0.17132 e^{0.93247/2} + 0.17132 e^{-0.93247/2} \\
&\approx 2.08436
\end{align}
And hence the required value will be \( (1/2)(2.08436) = 1.04218 \), whereas the exact result should be \( e^{1/2} -e^{-1/2} \approx 1.04219 \).
Leave a Reply