Gaussian quadrature

GUIDING QUESTION: How can I compute an integral numerically?

Motivation

How can we improve the Newton-Cotes rules?
We can increase the degree of precision of our quadrature rules by choosing the nodes optimally.

Déjà vu?

Sound familiar?
Recall that finding better nodes for polynomial interpolation led us to the Chebyshev points.

Guiding philosophy: Find $n$ nodes $x_\color{var(--emphColor)}{1}, \ldots, x_\color{var(--emphColor)}{n}$ so that the quadrature rule $$\int_a^b f(x) dx \approx \sum_i w_i f(x_i)$$ has DOP $\bbox[3pt, border: 3pt solid var(--emphColor)]{2n - 1}$.

We are almost doubling the DOP of an $n$-point NC rule.
We follow Bradie's convention here, so we use $1$-based indexing for our nodes when discussing Gaussian quadrature rules.
When discussing NC rules, we used $0$-based indexing.

Fixing the DOP

We require the DOP of our $n$-point Gaussian quadrature rule to be $\bbox[3pt, border: 3pt solid var(--emphColor)]{2n - 1}$,
so our nodes and weights must satisfy \begin{align*} \int_a^b p(x) dx = \sum_{i = 1}^n w_i p(x_i), \end{align*} for all polynomials $p(x)$ with degree at most $\bbox[3pt, border: 3pt solid var(--emphColor)]{2n - 1}$.

Fixing the DOP

Every such polynomial can be written as $$p(x) = \sum_{j = 0}^{2n-1} a_j x^j,$$ for some coefficients $a_j$.
Therefore, our quadrature rule will integrate $p$ exactly if and only if it can integrate each $\phi_i(x) = x^k$ exactly for $k = 0, 1, \ldots, \bbox[3pt, border: 3pt solid var(--emphColor)]{2n - 1}$.

Fixing the DOP

We require the DOP of our $n$-point Gaussian quadrature rule to be $\bbox[3pt, border: 3pt solid var(--emphColor)]{2n - 1}$,
so our nodes and weights must satisfy \begin{align*} \int_a^b x^k dx = \sum_{i = 1}^n w_i p(x_i), \end{align*} for $k = 0, 1, \ldots, \bbox[3pt, border: 3pt solid var(--emphColor)]{2n - 1}.$
Our nodes and weights satisfy a nonlinear system. We forced the DOP to be $\bbox[3pt, border: 3pt solid var(--emphColor)]{2n - 1}$ because this choice gives $2n$ equations to solve for the $2n$ variables $x_\color{var(--emphColor)}{1}, \ldots, x_\color{var(--emphColor)}{n}$ and $w_\color{var(--emphColor)}{1}, \ldots, w_\color{var(--emphColor)}{n}$.

Let's simplify things

From now on, we'll focus on integrals over the interval $[\color{var(--emphColor)}{-1}, \color{var(--emphColor)}{1}]$.
The change of variables $$x = \frac{b-a}{2}t + \frac{a+b}{2}$$ transforms an integral over any interval $[a, b]$ into an integral over $[\color{var(--emphColor)}{-1}, \color{var(--emphColor)}{1}]$:
$$\int_a^b f(x) d\color{var(--emphColor)}{x} = \int_{\color{var(--emphColor)}{-1}}^\color{var(--emphColor)}{1} f\bigg( \frac{b-a}{2}t + \frac{a+b}{2} \bigg) d\color{var(--emphColor)}{t} $$

Finding the nodes

Now notice that \begin{equation*} \int_{-1}^1 x^k dx = \begin{cases} {2 \over k+1}, & k \text{ is even,} \\ 0, & k \text{ is odd.} \end{cases} \end{equation*}

Finding the nodes

Therefore, the system nodes and weights we seek satisfy \begin{align} w_1 + \cdots + w_n &= {2 \over \color{var(--emphColor)}{0} + 1} \\ w_1 x_1^\color{var(--emphColor)}{1} + \cdots + w_n x_n^\color{var(--emphColor)}{1} &= \,\,\,\,\, \color{var(--emphColor)}{0}\\ w_1 x_1^\color{var(--emphColor)}{2} + \cdots + w_n x_n^\color{var(--emphColor)}{2} &= {2 \over \color{var(--emphColor)}{2} + 1} \\ w_1 x_{1}^\color{var(--emphColor)}{3} + \cdots + w_n x_n^\color{var(--emphColor)}{3} &= \,\,\,\,\, \color{var(--emphColor)}{0} \\ &\,\,\vdots \end{align}
There's one row for each power $\color{var(--emphColor)}{k} = 0, 1, \ldots, 2n-1$.
This is a nonlinear system in the nodes, but it is linear in the weights $w_j$, with Vandermonde coefficient matrix!

Legendre polynomials

It turns out that the nodes we seek are roots of the Legendre polynomials,
defined by the recursive relation \begin{equation*} p_{n+1}(x) = \frac{2n+1}{n+1}xp_n(x) - \frac{n}{n+1}p_{n-1}(x), \end{equation*}
with $p_0(x) = 1$, $p_1(x) = x$.

Legendre polynomials

The Legendre polynomial $p_n$ is even when $n$ is even, and it is odd when $n$ is odd.
Just like the Chebyshev polynomials!
All the roots of $p_n$ lie in the interval $[-1, 1]$, so we may always assume the placement of the nodes is symmetric about $x = 0$.
We may also assume symmetry in the weights!
($2$-point Gauss quadrature) In this case, symmetry implies $$x_2 = -x_1, \text{ and } w_1 = w_2.$$
The first equation in the DOP system gives $$w_1 + w_2 = 2,$$ so $w_1 = w_2 = 1$.
The third equation then gives \begin{equation*} 2x_1^2 = w_1 x_1^2 + w_2x_2^2 = \frac{2}{3} \end{equation*}
($2$-point Gauss quadrature) So the $2$-point Gaussian quadrature approximation suggests $$\int_\color{var(--emphColor)}{-1}^\color{var(--emphColor)}{1} f(x) dx \approx f\bigg({-1 \over \sqrt{3}}\bigg) + f\bigg({1 \over \sqrt{3}}\bigg).$$
In general, it is best to find the nodes first and then compute the weights by solving a linear system.
Like the weights defining an NC rule, the weights defining Gaussian quadrature rules are also given by definite integrals of appropriate Lagrange polynomials!
However, it's easier to obtain the weights by solving the linear system than to compute the definite integrals.
Congratulations! You reached the end of this lecture.
Take me home