Gaussian Quadratures

Last time I talked about quadratures and orthogonal polynomials, so now it’s time to combine the two.

Theorem: A quadrature scheme with $n$ interpolation points will be exact for polynomials of order up to $2n-1$ if the quadrature points are chosen to be the roots of the associated $n$ th-order orthogonal polynomial.

Note that, for an arbitrary choice of sample points, we would generally expect $n$ points to give a quadrature for polynomials of order $n-1,$ so $2n-1$ is a substantial improvement. For example, a single interpolation point can exactly integrate over polynomials of order one rather than just zero — as we had in our original example — and two interpolation points can integrate over polynomials of order three rather than one, even better than our original example. Three points can integrate up to order five!

Proof: We can split any $(2n-1)$ -order polynomial $P_{2n-1}(x)$ according to the orthogonal polynomial $f_n(x).$ In particular, we can use the decomposition

$P_{2n-1}(x) =P_{n-1}(x)f_n(x)+Q_{n-1}(x),$

where $P_{n-1}$ and $Q_{n-1}$ are polynomials of order $n-1.$ Since $f_n$ is predetermined, we must specify a $(2n-1)$ -order polynomial, with $2n$ coefficients, by choosing the polynomials $P_{n-1}$ and $Q_{n-1},$ with $2n$ coefficients between them. We have as many free coefficients to choose as we need, so this decomposition is always possible.

When we integrate over this expression, we see that

$\int_a^b P_{2n-1}(x)\, \textrm{d}x =\int_a^b (P_{n-1}(x) f_n(x)+Q_{n-1}(x) ) \,\textrm{d}x =\int_a^b Q_{n-1}(x) \,\textrm{d}x ,$

by $f_n(x)$ being orthogonal to any $P_{n-1},$ regardless of our choice of integration points. We therefore have the equality

$\int_a^b P_{2n-1}(x)\, \textrm{d}x =\int_a^b Q_{n-1}(x) \,\textrm{d}x = \sum_{k=1}^n \omega_k Q_{n-1}(x_k) ,$

so we must show that $\sum_{k=1}^n \omega_k P_{2n-1}(x_k)$ is equal to any one of these expressions. But, if we take $x_k$ to be the roots of $f_n(x) ,$ then equality to the third expression follows immediately, since the first term in the decomposition of $P_{2n-1}(x_k)$ disappears for every $x_k,$ and so $P_{2n-1}(x_k) = Q_{n-1}(x_k). \, \square$

We can, therefore, choose our interpolation points to be the roots of an orthogonal polynomial to optimise our quadrature method, with respect to the order of polynomials whose integrals are computed exactly. For our case over the interval $[-1,1],$ this gives the first optimal sets of interpolation points as

$\left\{0\right\}, \, \left\{-\frac{1}{\sqrt{3}} ,\frac{1}{\sqrt{3}}\right\}, \, \left\{-\sqrt{\frac{3}{5}} ,0,\sqrt{\frac{3}{5}} \right\}, \ldots$

Weighted Quadratures

As a final extension, suppose the function we are integrating is not polynomial, but it is a polynomial times some other function — $w(x),$ say — that is explicitly known. In this case, we wish to estimate

$\int_a^b w(x)f(x)\, \textrm{d}x.$

For this case, we can proceed as before, except that now we choose polynomials $f_i(x)$ that are orthogonal, in the sense that

$\langle f_i,f_j\rangle = \int_a^b w(x)f_i(x)f_j(x)\, \textrm{d}x = 0.$

In other words, the known function $w(x)$ becomes a weight function inside the integral, changing our choice of orthogonal polynomials, and of interpolation points. Common examples on the space $[0,1]$ are the Jacobi polynomials for $w(x) = (1-x)^{\alpha}(1+x)^{\beta},$ and the Chebyshev polynomials of first type for $w(x) = 1/\sqrt{1-x^2}.$

That’s about all I have to say. Gaussian quadratures are a reasonably straightforward step-up from basic quadratures, but give a sizeable increase in efficiency if we’re dealing with some function of known type, with an unknown polynomial term.

The Accidental Statistician

Occasional ramblings on statistics

Gaussian Quadratures

Leave a Reply Cancel reply

Share this:

Like this:

Related

Leave a Reply Cancel reply