Gaussian quadrature and orthogonal polynomials

The proofs here are based on the lecture notes from a course in approximation theory [1]. I think most of the results and proofs are standard, so I expect most of the texts on this subject to contain similar results and proofs.

Suppose is a nonnegative weight function and we want to compute the integral

for many different functions .

Depending on what and are, we might be able to evaluate the integral analytically (i.e. by integrating and obtaining an exact answer). However, in computational contexts it is often both more efficient and simpler to use a quadrature rule.

Definition 1: A quadrature rule is a collection of weights and positions such that

If, for some , the approximation is exact (i.e. the approximation symbol in the definition can be replaced by an equality symbol), we say that the quadrature rule is exact for .

Since the approximation given by the quadrature rule only depends on the value of on , the approximation will be better for smoother functions, for which the value of doesn’t fluctuate too wildly, so the behaviour is captured well by it’s value on some positions.

In practice, we often require the quadrature rule to be exact for some set of functions, often polynomials up to some maximum degree. We call this the exactness conditions. We can write the exactness condition as a finite number of equations by taking a basis of . By the linearity of the integral the quadrature rule is also exact for any linear combination of , which means that the quadrature rule is exact for any .

In the following, we will denote the space of polynomials with degree up to as .

Quadrature rules, the easy way

It is not too complicated to derive quadrature rules that are exact for any polynomial for a given . Given a function and distinct points , we can use Lagrange interpolation to construct a polynomial of degree that interpolates at . That is,

The polynomial is defined as , where are the Lagrange polynomials

We can now integrate the interpolation polynomial of as

Bringing the summation and the terms out of the integral, we see that

So, if we set , we have

Now, since interpolates , we expect that , so that

So we have found a quadrature rule that is exact for any .

It is also possible to find the weights in a more computational way. If we want the quadrature rule to be exact for function space of dimension , we can pick a basis of functions, pick distinct positions , write out the exactness requirements, and interpret the resulting system of equations with the weights as the unknowns.

For example, suppose that we want the quadrature rule to be exact for any . We pick basis functions (for example, would work) and distinct positions . Then, the exactness requirements become

which is a linear system of equations and unknowns . In particular, if we use the resulting system is a Vandermonde system.

For now, we have regarded the positions as given. It is usually a good choice to pick as Chebyshev nodes, since this ensures that the polynomial that interpolates is close to (it minimizes oscillatory behaviour of ).

If we change our point of view and interpret the positions as unknowns we can interpret the exactness requirements as a nonlinear system in unknowns. If we expect this system to behave like a non-signular linear one (there is no good reason to expect this, but let’s be optimistic), we would expect to be able to solve a system of equations with it. This is by analogy with a linear system, where a non-singular system has a unique solution if there are as many unknowns as equations.

So, this would mean that quadrature rules with points exists that are exact for any . These quadrature uses roughly half the number of points of the quadrature rules derived via the method presented earlier, at the cost of not having the freedom of choosing the positions freely.

In fact, it is possible to show that such quadrature rules exist, using orthogonal polynomials.

Orthogonal polynomials

Consider the following inner product on a space of functions:

Like before, is a nonnegative weight function. This forms a vector space, and the usual notion of orthogonality applies.

Definition 2: A sequence of polynomials is orthogonal if

  1. whenever

Like with vectors, we call the sequence orthonormal if additionally .

In addition, we will need the following two lemma’s.

Lemma 3: If is a sequence of orthogonal polynomials and is a polynomial with , then

Proof: The polynomials form a basis of . Since , we can write . It follows that . Since is an orthogonal sequence of polynomials and , all terms in the summation are zero, so it follows that .

The following lemma states that polynomials can be easily expressed as a linear combination of polynomials from a sequence of orthonormal polynomials.

Lemma 4: If is a sequence of orthonormal polynomials and is a polynomial with , then

Proof: Since is a basis of and , we can write as for some . Taking the inner product of and , we see

So . Substituting this in gives the required result.

A sequence of orthonormal polynomials can be computed by using Gram-Schmidt method, but it is more efficient to use the following result.

Theorem 5: If is a sequence of orthonormal polynomials then

Proof: From lemma 2 we have . Now

So we have

By lemma 1, is zero when , so are the only for which is nonzero. So we have

Using theorem 5, we can compute the orthonormal polynomial from the previous polynomials as follows. First we compute a polynomial which is orthogonal to but not normalized as

Then we compute by normalizing :

There is the following nice result about the roots of orthogonal polynomials.

Lemma 6: Let be a sequence of orthogonal polynomials with respect to the inner product , and let the support of be contained in . Then has distinct real roots .

Proof: Let be the zeroes of with odd multiplicity and lying in . Define . Now does not change sign on , so . By lemma 1, it follows that . So , which implies that has distinct roots in .

Now we are ready to prove the following theorem.

Theorem 7: Let be a sequence of orthogonal polynomials with respect to the inner product , and let be the zeroes of . Define for with the Lagrange polynomials for . Then


Proof: Let be the polynomial of degree that interpolates at . If we can show that , it will follow that .

First, we note that assumes the same values as for . So has zeroes at . Since are the zeroes from , the function has the same zeroes as and must be some multiple of , say, . Then we can write for some .

Now, we can write the integral as . Since we have by lemma 1. So we see that . It was shown earlier that , so it follows that

Sometimes, we want to compute the integral

For many different functions and . Using theorems 5 and 7 (and some root-finding method), we can compute a different quadrature rule for every . The inconvenient thing about this method is that the positions will be different for every . In this case, the benefit of being able to re-use the quadrature points probably outweighs the cost of having more quadrature points.

As described in [2], this can be applied to compute quadrature rules that can be used for the finite element method (FEM). The resulting quadrature rules are called ‘weighted quadrature rules’.


[1] Approximation Theory. Walter Groenevelt. Lecture notes WI4415, TU Delft, 2016.

[2] Matrix-free weighted quadrature for a computationally efficient isogeometric -method. Giancarlo Sangalli and Mattia Tani, 2017.