Gaussian quadrature

Introduction

Numerical integration of function is pervasive in scientific computations. It is also termed as quadrature and is generally evaluated on discretized grids. The simplest example would be

$$I={\int }_{-1}^1\text{dd}xf(x)\approx \sum _{n=1}^N{w}_{n}f(x_n)$$

where $x_n \in [-1,1]$ with integration weights $w_n$. An $N$-point grid with fixed $\{x_n\}$ points can be constructed to ensure exact integration for functions with Taylor polynomials up to order $N-1$. This can be shown by comparing the expansion of the exact and approximated forms:

$$\begin{array}{rl}I\equiv & {\int }_{-1}^1\text{dd}xf(x)={\int }_{-1}^1\text{dd}x\sum _{p=0}^{\mathrm{\infty }}\frac{f^{(p)}(0)}{p!}{x}^p=\sum _{p=0}^{\mathrm{\infty }}\frac{f^{(p)}(0)}{p!}{\int }_{-1}^1\text{dd}x{x}^p\\ \approx & \sum _{n=1}^N{w}_{n}f(x_n)=\sum _{n=1}^N{w}_n\sum _{p=0}^{\mathrm{\infty }}\frac{f^{(p)}(0)}{p!}{x}_n^p=\sum _{p=0}^{\mathrm{\infty }}\frac{f^{(p)}(0)}{p!}\sum _{n=1}^N{w}_n{x}_n^p\end{array}$$

It is exact if for any $p\in\mathbb{N}$, there exists

$$\sum _{n=1}^N{w}_n{x}_n^p={\int }_{-1}^1\text{dd}x{x}^p$$

For a fixed $\{x_n\}$, this is equivalent to solving a set of linear equations for $\{\omega_n\}$. A unique solution can be found when $p$ is truncated at $N-1$, i. e. $N$ equations for $N$ variables.

If we relax the constraint of fixed nodes $\{x_n\}$ and treat them as unknown variables, the degrees of freedom becomes $2N$ and $p$ needs to go to $2N-1$. This leads to a better estimate of an integral, as it now can integrate exactly polynomials up to order of $2N-1$. It is proposed by Gauss and hence called Gaussian quadrature.

Gauss-Legendre quadrature

To determine the nodes $\{x_n\}$ and weights $\{w_n\}$ for Gaussian quadrature, we consider the integration of polynomials up to $2N-1$ order. Such function can be written in the following form

$$f(x)=Q(x)P_N(x)+R(x)$$

where $P_N(x)$ is the $N$-th-order Legendre polynomial, and $Q(x)$ and $R(x)$ are general polynomials of $N-1$ or less order. According to the orthogonality of Legendre polynomials

$${\int }_{-1}^1\text{dd}x{P}_i(x)P_j(x)=\frac{2}{2i+1}{\delta }_{ij},$$

$P_N$ is orthogonal to $Q$ and $R$, therefore the formal integral

$$I\equiv {\int }_{-1}^1\text{dd}xf(x)={\int }_{-1}^1\text{dd}x[Q(x)P_N(x)+R(x)]={\int }_{-1}^1\text{dd}xR(x)$$

For the $n$-series, we assign the roots of $P_{N}(x)$ to $\{x_n\}$, so that $P_N(x_n) = 0, \forall n = 1, \cdots N$

$$\sum _{n=1}^N{w}_{n}f(x_n)=\sum _{n=1}^N{w}_n[Q(x_n)P_N(x_n)+R(x_n)]=\sum _{n=1}^N{w}_{n}R(x_n)$$

These two equations are the same as Eq. $\text{(???)}$ except with $R(x)$ rather than $f(x)$. Since now the nodes are settled down, we only have to solve the weights. They can be solved by letting

$$f(x)=\frac{P_N(x)P_N^{\prime }(x)}{x-x_n}$$

where $x_n$ is the $n$-th root of $P_N(x)$, and $P’_N$ being its derivative. Insert it to the $n$-series

$$\begin{array}{rl}\sum _{n^{\prime }=1}^N{w}_{n^{\prime }}f(x_{n^{\prime }})=& \sum _{n^{\prime }=1,n^{\prime }\ne n}^N{w}_{n^{\prime }}\frac{P_N(x_{n^{\prime }})P_N^{\prime }(x_{n^{\prime }})}{x_{n^{\prime }}-x_n}+w_n\underset{x\to x_n}{lim}\frac{P_N(x)P_N^{\prime }(x)}{x-x_n}\end{array}$$

Each term in the summation over $n\ne n’$ is zero as we have chosen $\{x_n\}$ be the roots of $P_{N}$ and the denominator is nonzero. The limit, using L’Hospital’s rule, is

$$\begin{array}{r}\underset{x\to x_n}{lim}\frac{P_N(x)P_N^{\prime }(x)}{x-x_n}=\underset{x\to x_n}{lim}{[P_N(x)P_N^{\prime }(x)]}^{\prime }=P_N^{\prime }(x_n)P_N^{\prime }(x_n)+P_N(x_n)P_N^{″}(x_n)={[P_N^{\prime }(x_n)]}^2\end{array}$$

Therefore

$$\sum _{n^{\prime }=1}^N{w}_{n^{\prime }}f(x_{n^{\prime }})=w_n{[P_N^{\prime }(x_n)]}^2$$

On the other hand, the integral

$$\begin{array}{rl}{\int }_{-1}^1\text{dd}xf(x)=& {\int }_{-1}^1\text{dd}x\frac{P_N(x)P_N^{\prime }(x)}{x-x_n}={\int }_{-1}^1\text{dd}[P_N(x)]\frac{P_N(x)}{x-x_n}\\ =& {\frac{P_N^2(x)}{x-x_n}|}_{-1}^1-{\int }_{-1}^1{P}_N(x)\text{dd}\frac{P_N(x)}{x-x_n}\end{array}$$

For the second term, as $x_n$ is a root of $P_N$, $P_N(x)/(x-x_n)$ must be a polynomial of order $N-1$. Taking derivative will lead to a polynomial of order $N-2$ in the integrand, which can always be expressed as a linear combination using Legendre polynomials of order no more than $N-2$. Thus the integral vanishes due to their orthogonality with $P_N$, and

$$\begin{array}{rl}{\int }_{-1}^1\text{dd}xf(x)=& {\frac{P_N^2(x)}{x-x_n}|}_{-1}^1=\frac{P_N^2(1)}{1-x_n}+\frac{P_N^2(-1)}{1+x_n}=\frac{2}{1-x_n^2}\end{array}$$

where we have used the property of Legendre polynomials: $P_N(1)=1, P_{N}(-1) = (-1)^N$. Equating Eqs. $\text{(???)}$ and $\text{(???)}$ results in

$$w_n=\frac{2}{(1-x_n^2){[P_N^{\prime }(x_n)]}^2}$$

which is the desired weight.

Gauss-Legendre quadrature is available in SciPy. The roots and weights can be obtained by calling roots_legendre in the scipy.special module. Below is an example to integrate $6x^8 + 4x^2 - 1$. The exact integral between $[-1,1]$ is 2.

import scipy
import numpy as np

def f(x):
    return 6. * x ** 8 + 4. * x ** 2 - 1.

for n in [2, 3, 4, 5, 6]:
    xs, ws = scipy.special.roots_legendre(n)
    integral = np.sum(ws * f(xs))
    print("N = {:d}, 2N-1 = {:2d}, integral = {:.10f}"
          .format(n, 2 * n - 1, integral))

N = 2, 2N-1 =  3, integral = 0.8148148148
N = 3, 2N-1 =  5, integral = 1.5306666667
N = 4, 2N-1 =  7, integral = 1.9303401361
N = 5, 2N-1 =  9, integral = 2.0000000000
N = 6, 2N-1 = 11, integral = 2.0000000000

Extension to any interval

For any function $f$ defined on closed interval $[a,b]$, we can do the following change of variables

$$x=\frac{b-a}{2}t+\frac{a+b}{2},$$

therefore

$$I={\int }_a^b\text{dd}xf(x)=\frac{b-a}{2}{\int }_{-1}^1\text{dd}tf(\frac{b-a}{2}t+\frac{a+b}{2})=\frac{b-a}{2}{\int }_{-1}^1\text{dd}tg(t).$$

The $t$ integral can then be approximated by the Gauss-Legendre quadrature

$$\begin{array}{rl}I\approx & \frac{b-a}{2}\sum _{n=1}^N{w}_{n}g(x_n)=\sum _{n=1}^N\left(\frac{b-a}{2}{w}_n\right)f(\frac{b-a}{2}{x}_n+\frac{a+b}{2})\end{array}$$

Therefore for interval $[a,b]$, the nodes and weights are mapped to

$$\begin{array}{rl}{w}_n\to & \frac{b-a}{2}{w}_n\\ x_n\to & \frac{b-a}{2}{x}_n+\frac{a+b}{2}\end{array}$$

It can be also extended to some open intervals. For example,

$$x=\frac{1+t}{1-t}+c$$

can be used for semi-infinite integral $[c, \infty)$. In this case,

$$\begin{array}{rl}I=& {\int }_0^{\mathrm{\infty }}\text{dd}xf(x)={\int }_0^{\mathrm{\infty }}\text{dd}\left(\frac{1+t}{1-t}\right)f(\frac{1+t}{1-t}+c)\\ =& {\int }_{-1}^1\text{dd}t\frac{2}{{(1-t)}^2}f(\frac{1+t}{1-t}+c)\approx \sum _{n=1}^N\frac{2w_n}{{(1-x_n)}^2}f(\frac{1+x_n}{1-x_n}+c)\end{array}$$

therefore

$$\begin{array}{rl}{w}_n\to & \frac{2}{{(1-x_n)}^2}{w}_n\\ x_n\to & \frac{1+x_n}{1-x_n}+c\end{array}$$

In practice, the semi-infinite integral can be divided into separate regions, and each region is integrated using its own Gauss-Legendre grids. [1]

See also

• Gaussian quadrature - Wikipedia
• Legendre polynomials - Wikipedia
• Gauss-Legendre quadrature - Wikipedia
• roots_legendre - SciPy Manual

References

[1] R. W. Godby, M. Schlüter, and L. J. Sham, Phys. Rev. B 37, 10159 (1988) [DOI].