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 $μ (x)$ is a nonnegative weight function and we want to compute the integral

\int f (x) μ (x) d x

for many different functions $f$ .

Depending on what $f$ 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 $w_{1}, w_{2}, . . ., w_{n}$ and positions $x_{1}, x_{2}, . . ., x_{n}$ such that

k = q \sum n w_{k} f (x_{k}) \approx \int f (x) μ (x) d x

Since the approximation given by the quadrature rule only depends on the value of $f$ on $x_{1}, x_{2}, . . ., x_{n}$ , the approximation will be better for smoother functions, for which the value of $f$ doesn’t fluctuate too wildly, so the behavior is captured well by its value on some positions.

In practice, we often require the quadrature rule to be exact for some set $S$ of functions. Often, we take $S$ as the set of polynomials up to some maximum degree. We call this the exactness conditions. We can write the exactness conditions as a finite number of equations by taking a basis $f_{1}, f_{2}, . . .$ of $S$ . By the linearity of the integral the quadrature rule is also exact for any linear combination of $f_{1}, f_{2}, . . .$ which means that the quadrature rule is exact for any $f \in S$ .

In the following, we will denote the space of polynomials with degree up to $n$ as $P_{n}$ .

Quadrature rules, the easy way

It is not too complicated to derive quadrature rules that are exact for any polynomial $p \in P_{n}$ for a given $n$ . Give a function $f$ and $d + 1$ distinct points $x_{0}, x_{1}, . . ., x_{d}$ we can use Lagrange interpolation to construct a polynomial $f_{d}$ of degree $d$ that interpolates $f$ at $x_{0}, x_{1}, . . ., x_{d}$ . That is,

f_{d} (x_{k}) = f (x_{k}) for k = 0, 1, ..., d

The polynomial $f_{d}$ is defined as $f_{d} (x) = \sum_{k = 0}^{n} f (x_{k}) l_{k} (x)$ , where $l_{k}$ are the Lagrange polynomials defined by

l_{k} (x) = \frac{( x - x _{0} ) ( x - x _{1} ) \dots ( x - x _{k - 1} ) ( x - x _{k + 1} ) \dots ( x - x _{d} )}{( x _{k} - x _{0} ) ( x _{k} - x _{1} ) \dots ( x _{k} - x _{k - 1} ) ( x _{k} - x _{k + 1} ) \dots ( x _{k} - x _{d} )}

We can now integrate the interpolation polynomial $f_{d}$ of $f$ as

\int f_{d} (x) μ (x) d x = \int k = 0 \sum n f (x_{k}) l_{k} (x) μ (x) d x

Bringing the summation and the terms $f (x_{k})$ outside of the integral, we see

\int f_{d} (x) μ (x) d x = k = 0 \sum n f (x_{k}) \int l_{k} (x) μ (x) d x

So if we set $w_{k} = \int l_{k} (x) μ (x) d x$ , we have

\int f_{d} (x) μ (x) d x = k = 0 \sum n w_{k} f (x_{k})

Now, since $f_{d}$ interpolates $f$ , we expect that $f_{d} \approx f$ , so that

\int f (x) μ (x) d x \approx \int f_{d} (x) μ (x) d x = k = 0 \sum n w_{k} f (x_{k})

So we have found a quadrature rule that is exact for any $f \in P_{d}$ .

It is also possible to find the weights in a more computational way. If we want the quadrature rule to be exact for a function space of dimension $d$ , we can pick a basis of $d$ functions, pick $d$ distinct positions $x_{1}, x_{2}, . . ., x_{d}$ , write out the exactness conditions, and solve the resulting linear system of equations for the weights.

For example, suppose that we want the quadrature rule to be exact for any $p \in P_{d}$ . We pick $d + 1$ basis functions $b_{0}, b_{1}, . . ., b_{d}$ (for example $b_{k} (x) = x^{k}$ would work) and distinct positions $x_{0}, x_{1}, . . ., x_{d}$ . Then the exactness requirements become

j = 0 \sum d w_{j} b_{k} (x_{j}) = \int b_{k} (x) μ (x) d x for k = 0, 1, ..., d

which is a linear system of $d + 1$ equations in the $d + 1$ unknowns $w_{0}, w_{1}, . . . w_{d}$ . In particular, if we use $b_{k} (x) = x^{k}$ the resulting system is a Vandermonde system.

For now, we have regarded the positions $x_{1}, x_{2}, . . ., x_{d}$ as given. It is usually a good choice to pick them as Chebyshev nodes, since this ensures that the polynomial $f_{d}$ that interpolates $f$ is actually close to $f$ – it minimizes oscillatory behavior of $f_{d}$ .

If we change our point of view and interpret the positions $x_{1}, x_{2}, . . ., x_{d}$ as unknowns we can interpret the exactness requirements as a nonlinear system in $2 d$ unknowns. If we expect this system to behave like a non-singular linear one (which is not a given, but let’s be optimistic), we would expect to be able to solve a system of $2 d$ 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 for each $n$ there exists a quadrature rule with $n$ points that is exact for $P_{2 n - 1}$ . 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 $x_{0}, x_{2}, . . ., x_{n - 1}$ freely.

It is possible to derive quadrature rules like this using orthogonal polynomials.

Orthogonal polynomials

Consider the following inner product on a space of functions

⟨ f, g ⟩ = \int f (x) g (x) μ (x) d x

Like before, $μ (x)$ is a nonnegative weight function. This inner product makes the function space into a vector space.

The usual notions of orthogonality and orthonormality apply.

Definition 2: A sequence of polynomials $p_{0}, p_{1}, . . .$ is orthogonal if

$deg (p_{k}) = k$
$⟨ p_{i}, p_{j} ⟩ = 0$ whenever $i \neq = j$

Definition 3: A sequence of polynomials $p_{0}, p_{1}, . . .$ is orthonormal if it is orthogonal and $⟨ p_{k}, p_{k} ⟩ = 1$ for every $k = 0, 1, . . .$

We will use the following lemmas, which are easy to understand and prove.

Lemma 4: If $p_{0}, p_{1}, . . .$ is a sequence of orthogonal polynomials and $p$ is a polynomial with $deg (p) < n$ , then

⟨ p, p_{n} ⟩ = 0

Proof: The polynomials $p_{0}, p_{1}, . . ., p_{n - 1}$ form a basis of $P_{n}$ . Since $p \in P_{n}$ , we can write $p (x) = \sum_{k = 0}^{n - 1} c_{k} p_{k} (x)$ . It follows that $⟨ p, p_{n} ⟩ = ⟨ \sum_{k = 0}^{n - 1} c_{k} p_{k}, p ⟩ = \sum_{k = 0}^{n - 1} c_{k} ⟨ p_{k}, p_{n} ⟩$ . Since $p_{0}, p_{1}, . . .$ is an orthogonal sequence of polynomials and $k < n$ , all terms in the summation are zero, so it follows that $⟨ p, p_{n} ⟩ = 0$ .

$□$

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

Lemma 5: If $p_{0}, p_{1}, . . .$ is a sequence of orthonormal polynomials and $p$ is a polynomial with $deg (p) = n$ , then

p (x) = k = 0 \sum n ⟨ p, p_{k} ⟩ p_{k} (x)

Proof: Since $p_{0}, p_{1}, . . ., p_{n}$ is a basis of $P_{n}$ and $p \in P_{n}$ , we can write $p$ as $p (x) = \sum_{k = 0}^{n} c_{k} p_{k} (x)$ for some $c_{0}, c_{1}, . . ., c_{n} \in R$ . Taking the inner product of $p$ and $p_{k}$ , we see

⟨ p, p_{k} ⟩ = ⟨ j = 0 \sum n c_{j} p_{j}, p_{k} ⟩ = j = 0 \sum n c_{j} ⟨ p_{j}, p_{k} ⟩ = c_{k}

So $c_{k} = ⟨ p, p_{k} ⟩$ . Substituting this in $p(x) = \sum_{k = 0}^n c_k p_k(x) gives the required result.

$□$

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

Theorem 6: If $p_{0}, p_{1}, . . .$ is a sequence of polynomials then

⟨ p_{n}, x p_{n + 1} ⟩ p_{n + 1} (x) = x p_{n} (x) - ⟨ p_{n}, x p_{n - 1} ⟩ p_{n - 1} (x) - ⟨ p_{n} x p_{n} ⟩ p_{n} (x)

Proof: From lemma 5 we have $x p_{n} (x) = \sum_{k = 0}^{n + 1} ⟨ x p_{n}, p_{k} ⟩ p_{k}$ . Now

⟨ x p_{n}, p_{k} ⟩ = \int x p_{n} (x) p_{k} (x) μ (x) d x = ⟨ p_{n}, x p_{k} ⟩

So we have

x p_{n} = k = 0 \sum n + 1 ⟨ p_{n}, x p_{k} ⟩ p_{k}

By lemma 4, $⟨ p_{n}, x p_{k} ⟩$ is zero when $deg (x p_{k}) = k + 1 < n$ , so $k = n - 1, n, n + 1$ are the only $k \leq n + 1$ for which $⟨ p_{n}, x p_{k} ⟩$ is nonzero. So we have

x p_{n} = ⟨ p_{n}, x p_{n - 1} ⟩ p_{n - 1} (x) + ⟨ p_{n}, x p_{n} ⟩ p_{n} (x) + ⟨ p_{n}, x p_{n + 1} ⟩ p_{n + 1} (x)

$□$

Using theorem 6, we can compute the orthonormal polynomial $p_{n + 1}$ from the previous polynomials as follows. First, we compute a polynomial $q_{n + 1} (x)$ which is orthogonal to $p_{0}, p_{1}, . . ., p_{n}$ as

q (x) = x p_{n} (x) - ⟨ p_{n}, x p_{n - 1} ⟩ p_{n - 1} (x) - ⟨ p_{n}, x p_{n} ⟩ p_{n} (x)

Then, we compute $p_{n + 1}$ by normalizing $q_{n + 1}$ :

p_{n + 1} (x) = \frac{q _{n + 1} ( x )}{⟨ q _{n + 1} , q _{n + 1} ⟩}

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

Lemma 7: Let $p_{0}, p_{1}, . . .$ be a sequence of orthonal polynomials with respect to the inner product $⟨ f, g ⟩ = \int f (x) f (x) μ (x) d x$ , and let the support of $μ$ be contained in $[a, b]$ . Then $p_{n}$ has $n$ distinct real roots $x_{1}, x_{2}, . . ., x_{n} \in (a, b)$ .

Proof: Let $x_{1}, x_{2}, . . ., x_{m}$ be the zeroes of $p_{n}$ that have odd multiplicity and lie in $(a, b)$ . Define $q (x) = (x - x_{1}) (x - x_{2}) \dots (x - x_{m})$ . Now $q (x) p_{n} (x)$ does not change sign on $[a, b]$ , so $⟨ q, p_{n} ⟩ \neq = 0$ . It follows from $m \leq n$ and lemma 4 that $deg (q) = n$ . So $m = n$ , which implies that $p_{n}$ has $n$ distinct roots in $(a, b)$ .

Now we are ready to prove the following theorem.

Theorem 8: Let $p_{0}, p_{1}, . . .$ be a sequence of orthonal polynomials with respect to the inner product $⟨ f, g ⟩ = \int f (x) f (x) μ (x) d x$ and let $x_{1}, x_{2}, . . ., x_{n}$ be the zeroes of $p_{n}$ . Define $w_{k} = \int l_{k} (x) μ (x) d x$ for $k = 1, 2, . . ., n$ , where $l_{1}, l_{2}, . . ., l_{n}$ are the Lagrange polynomials for $x_{1}, x_{2}, . . ., x_{n}$ . Then

k = 1 \sum n w_{k} f (x_{k}) = \int f (x) μ (x) d x

for $f \in P_{2 n - 1}$ .

Proof: Let $f_{n} \in P_{n}$ be the polynomial of degree $n$ that interpolates $f$ at $x_{1}, x_{2}, . . ., x_{n}$ . If we can show that $\int (f (x) - f_{n} (x)) μ (x) d x = 0$ , it follows that $\int f (x) μ (x) d x = \int f_{n} (x) μ (x) d x$ .

First, we note that $f_{n} (x)$ assumes the same values as $f (x)$ for $x = x_{1}, x_{2}, . . ., x_{n}$ . So $f (x) - f_{n} (x)$ has zeroes at $x_{1}, x_{2}, . . ., x_{n}$ . Since $x_{1}, x_{2}, . . ., x_{n}$ are the zeroes from $p_{n}$ , the function $f (x) - f_{n} (x)$ has the same zeroes as $p_{n}$ and must be some multiple of $p_{n}$ , say $f (x) - f_{n} (x) = c p_{n} (x)$ . Then we can write $f (x) - f_{n} (x) = c p_{n} (x) q (x)$ for some $q \in P_{n - 1}$ .

Now, we can write the integral as $\int (f (x) - f_{n} (x)) μ (x) d x = c \int p_{n} (x) q (x) μ (x) d x = c ⟨ p_{n}, q ⟩$ . Since $deg (q) < n$ we have $⟨ p_{n}, q ⟩ = 0$ by lemma 4. So we see that $\int f (x) μ (x) d x = \int f_{n} (x) μ (x) d x$ . It was shown ealier that $\int f_{n} (x) μ (x) d x = \sum_{k = 1}^{n} w_{k} f (x_{k})$ , so it follows that

\int f (x) μ (x) d x = \int f_{n} (x) μ (x) d x = k = 1 \sum n w_{k} f (x_{k})

$□$

Sometimes, we want to compute the integral

\int p (x) q (x) μ (x) d x

for many different $p$ and $q$ . Using theorem 6 and 8 and a root-finding method, we can compute a different quadrature rule for every $q$ . The inconvenient thing about this method is that the quadrature points $x_1, x_2, … will be different for every $q$ .

In [2], it is proposed to instead pick one set of quadrature points and re-use these for every $q$ . This means we have to use more quadrature points per $q$ , but since we can use the same set of quadrature points we only have to evaluate $p$ on this one set, which is more efficient in the end. It also means that we can use the ‘simple’ way of computing quadrature rules outlined in the section ‘quadrature rules, the easy way’.

References

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

[2] Matrix-free weighted quadrature for a computationally efficient isogeometric k-method. Giancarlo Sangalli and Mattia Tani, 2017. https://arxiv.org/abs/1712.08565