Why is the Gaussian-Legendre Quadrature so effective?

Question

I understand how it works, how its derived, etc. The proof of it has been shown to me. That is to say, I know how Legendre polynomials are derived, I know they are orthogonal, I know we sample a function at the roots of the polynomial, I know how the weights are calculated, and I see that it is basically a Riemann sum in that we are approximating the area with rectangles.

But the Riemann left sum, midpoint, right sum, simpsons rule, etc., have all been shown to me too. But they are not as impressive. What I dont understand is what makes the G-L Q so effective. Where does it derive its stellar accuracy? And why should it be so accurate anyway? The derivation is an interesting proof of concept that a numeric approximation can be done this way and is interesting in its own right... but to see that it is drastically so much more accurate than a Riemann sum is a completely different matter.

I dont see why the orthogonality of the Legendre polynomials is so important and how it plays such a key role in the accuracy. I dont see why we choose the roots of the polynomials as the sampling points of a function, as opposed to any other set of points with appropriate weights derived. Is the norm used in deriving the Legendre polynomials unique or can other norms and polynomials prove equally (or more) effective? Does the G-L Q work by approximating an arbitrary function by a polynomial? And if so, wouldnt we get more accuracy making those polynomials generic rather than by restricting them to linear combinations of Legendre?

Answer 1

I find it natural. We approximate an integral of function $f(x)$ over interval $[-1,1]$: \begin{align} \int_{-1}^1f(x)\,\text{d}x\approx\sum_{i=1}^nw_if(x_i),~~~~~~~(1) \end{align} where $x_i,~i=1,\dots,n$ are points and $w_i,~i=1,\dots,n$ are weights. Other quadrature rules (Newton-Cottes quadrature, Clenshaw-Curtis quadrature etc.) see Eq.(1) as an optimization problem \begin{align} \min_{w_1,\dots,w_n}\left[\max_{f(x) \in \left\{1,x,\dots,x^{2n-1}\right\}}\Bigg(\int_{-1}^1f(x)\,\text{d}x-\sum_{i=1}^nw_if(x_i)\Bigg)^2\right], \end{align} where $x_i$ are given and only $w_i$ are free. The Gauss quadrature sees Eq.(1) as an optimization problem, where both $x_i$ and $w_i$ are free: \begin{align} \min_{x_1,\dots,x_n,~w_1,\dots,w_n}\left[\max_{f(x) \in \left\{1,x,\dots,x^{2n-1}\right\}}\Bigg(\int_{-1}^1f(x)\,\text{d}x-\sum_{i=1}^nw_if(x_i)\Bigg)^2\right], \end{align} so you can find in general a better optimum (=satisfy Eq.(1) for more polynomials exactly). The fact, that $x_i$ are root's of the Legendre polynomial, can be seen as an interesting coincidence.

Answer 2

There is a key estimate which is commonly used in the analysis of quadrature rules for relatively nice functions:

$$|I(f)-Q(f)| \leq |I(f)-I(p)| + |I(p)-Q(p)| + |Q(p)-Q(f)|$$

where $I$ is the exact integral, $Q$ is the quadrature rule, $f$ is your integrand of interest, and $p$ can be any reasonable function, but for this bound to be useful it should be close to $f$ in some suitable sense. Generally we consider the case where $p$ is some polynomial.

Now on a bounded domain with an integrable weight function, Weierstrass's theorem lets you make the first term small as long as you're allowed to choose $p$ to be a polynomial with a large degree. If $Q$ is exact for polynomials of degree all the way up to the degree of $p$, then the second term is just zero. Thus many different rules behave nicely in the first two terms.

If you proceed this way, then the issue is the third term. In general, you can do this computation:

$$|Q(p)-Q(f)|=|Q(p-f)|=\left | \sum_{i=1}^n w_i (p(x_i)-f(x_i)) \right | \leq \left ( \max_i |p(x_i)-f(x_i)| \right ) \sum_{i=1}^n |w_i|.$$

Again Weierstrass's theorem lets you make that maximum as small as you want. But what about the sum? In general there is no reason to expect that sum to be bounded. But it is bounded if the weights have all the same sign, because $I(1)=Q(1)=\sum_{i=1}^n w_i$. For Gauss quadrature these weights are in fact positive, as you can show by noting that

$$0<I \left ( \prod_{j \neq i} (x-x_j)^2 \right ) = Q \left ( \prod_{j \neq i} (x-x_j)^2 \right )=w_i \prod_{j \neq i} (x_i-x_j)^2$$

which is because these polynomials all have degree $2n-2$, are everywhere nonnegative, and vanish at all but one of the nodes.

Note that from this point of view, the theory of orthogonal polynomials doesn't really come into play; the fact that Gauss quadrature converges follows from merely knowing the exactness behavior. Orthogonal polynomials only really come into play when you ask for the nodes, weights, and convergence rate.

There is a subtlety that I glossed over. What if you have some other method and instead of using Weierstrass's theorem to select $p$ so that the first term is small, you select $p$ to just interpolate $f$, so that the third term is zero? In this case the second term can again be made to be zero. But does the first term necessarily go to zero? The answer is no, which is closely related to Runge's phenomenon, which you can read more about on Wikipedia.