Introduction to Quadrature

In this section, we aim to approximate the definite integral

\[ \mathcal{I}(f) := \int_a^b f(x) dx \]

on the interval \([a,b] \subset \mathbb{R}\), using a method known as quadrature. A quadrature formula is an approximation of the form

\[ \mathcal{Q}_n(f) := \sum_{k=0}^n \sigma_k f(x_k), \]

where \(\{x_k\}\) are called nodes, and \(\sigma_k\) are corresponding weights or coefficients. The choice of nodes and coefficients is crucial for the accuracy of the approximation. Importantly, quadrature formulas share a key property with integrals: linearity. This means

\[ \mathcal{Q}_n(\alpha f + \beta g) = \alpha \mathcal{Q}_n(f) + \beta \mathcal{Q}_n(g) \quad \forall \, \alpha, \beta \in \mathbb{R}. \]

The degree of exactness of a quadrature formula \(\mathcal{Q}_n\) is defined as the highest degree \(r\) for which it exactly integrates all polynomials up to degree \(r\), but not for degree \(r+1\). Formally,

\[\begin{split} \begin{aligned} \mathcal{Q}_n(x^j) & = \mathcal{I}(x^j) \quad \text{for } j = 0, \ldots, r,\\ \mathcal{Q}_n(x^{r+1}) & \neq \mathcal{I}(x^{r+1}). \end{aligned} \end{split}\]

Therefore, a quadrature formula \(\mathcal{Q}_n\) with a degree of exactness \(r\) perfectly approximates the integral of any polynomial in \(\mathcal{P}_r\), the set of all polynomials of degree up to \(r\).

Interpolatory Quadratures

Interpolatory quadrature involves approximating a function \(f\) by a polynomial \(p\) that interpolates \(f\) at given nodes \(\{x_k\}\), and then integrating this polynomial. Utilising Lagrange interpolation, \(p(x)\) is expressed as

\[ p(x) = \sum_{k=0}^n f(x_k) L_k(x), \quad \text{where } L_k(x) = \prod_{l\ne k} \frac{x-x_l}{x_k-x_l}. \]

This leads to the quadrature formula

\[ \mathcal{Q}_n(f) := \int_a^b p(x) dx = \sum_{k=0}^n f(x_k) \int_a^b L_k(x) dx, \]

with the weights

\[ \sigma_k = \int_a^b L_k(x) dx. \]

Hence, the interpolatory quadrature formula becomes \(\mathcal{Q}_n(f) = \sum_{k=0}^n \sigma_k f(x_k)\). The degree of exactness of \(\mathcal{Q}_n\) is at least \(n\), and it for \(n \geq 0\) holds that

\[ \sum_{k=0}^n \sigma_k = b-a, \qquad\mathcal{Q}_n(1) = \mathcal{I}(1). \]

In the case of equidistant nodes, we obtain the Newton-Cotes formulas, and if these nodes include the endpoints \(a\) and \(b\), the formulas are termed closed Newton-Cotes formulas.

Example 6 (Interpolatory quadrature formulas)

• Rectangle Method:

  \[ \mathcal{Q}_0(f) = (b-a)f\Bigl(\frac{a+b}{2}\Bigr). \]

• Trapezoid Rule (a closed Newton-Cotes formula for \(n=1\)):

  \[ \mathcal{Q}_1(f) = (b-a)\frac{f(a)+f(b)}{2}. \]

• Simpson’s Rule (a closed Newton-Cotes formula for \(n=2\)):

  \[\mathcal{Q}_2(f) = (b-a)\,\Bigl[ \frac{1}{6} f(a) + \frac{2}{3} f\Bigl(\frac{a+b}{2}\Bigr) + \frac{1}{6} f(b) \Bigr].\]

Error bound

Using the interpolation error bound, we obtain

\[\begin{split} \begin{aligned} |\mathcal{I}(f)-\mathcal{Q}_n(f)| &= |\mathcal{I}(f)-\mathcal{Q}_n(p)| = |\mathcal{I}(f)-\mathcal{I}(p)| = |\mathcal{I}(f-p)|\\ &= \Bigl| \int_a^b \frac{\omega_{n+1}(x)}{(n+1)!}\,f^{(n+1)}(\xi(x)) dx \Bigr|\\ &\le \frac1{(n+1)!}\,\max_{\xi\in[a,b]} |f^{(n+1)}(\xi)| \,\int_a^b |\omega_{n+1}(x)| dx. \end{aligned} \end{split}\]

This error bound suggests that the accuracy of \(\mathcal{Q}_n\) is limited both by the differentiability of \(f\), which is usually out of our control as numerical analysts, and by the size of \(|\omega_{n+1}|\), which depends on the choice of the nodes \(\{x_k\}\) and thus may be improved. We will address these two issues in the next section.

Runge phenomenon

Even when \(|\omega_{n+1}|\) is bounded, and when \(f\) is smooth, interpolatory quadrature may still suffer from Runge’s phenomenon. We find negative values for some coefficients \(\sigma_k^{(n)}\) even for moderate values of \(n\), 8 in the case of closed Newton-Cotes formulas. In fact, one can prove that, for the closed Newton-Cotes formulas,

\[ \sum_{k=0}^n \bigl|\sigma_k^{(n)}\bigr| \to \infty \qquad n \to \infty. \]

This makes Newton-Cotes formulas vulnerable to instability for large \(n\).

Python skills

Here is an implementation of the trapezoidal and Simpson’s rule in Python.

def trapezoidal_rule(f, a, b):
    return (b - a) / 2 * (f(a) + f(b))

def simpsons_rule(f, a, b):
    return (b - a) / 6 * (f(a) + 4 * f((a + b) / 2) + f(b))

# Example usage
# Define a function to integrate
def my_function(x):
    return x ** 2

# Apply the quadrature rules
a = 0  # Lower limit
b = 1  # Upper limit
print("Trapezoidal Rule:", trapezoidal_rule(my_function, a, b))
print("Simpson's Rule:", simpsons_rule(my_function, a, b))

Self-check questions

Question

Calculate the weights and nodes for the closed Newton-Cotes formula with \(n = 3\).

Answer

The nodes are \(x_0 = a\), \(x_1 = a + h\), \(x_2 = a + 2h\) and \(x_3 = b\).

The weights are determined with the integration of the Lagrange interpolating polynomial:

• \(\sigma_0 = \sigma_3 = 3h/8\)

• \(\sigma_1 = \sigma_2 = 9h/8\)

Thus, the integral approximation over the interval \([a, b]\) is

\[\begin{split} \begin{align} \int_a^b f(x) \, dx & \approx \frac{3h}{8} [f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3)]\\ & = \frac{b - a}{8} [f(a) + 3f(a + h) + 3f(a + 2h) + f(b)]. \end{align} \end{split}\]

The rule is known as Simpson’s 3/8 Rule.

Question

A quadrature rule \(\mathcal{Q}\) has a negative weight \(\sigma_k\). Show that there exists a non-negative function \(g\) such \(\mathcal{Q}(g) < 0\). This means that the quadrature rule assigns a negative estimate of the integral to the non-negative function.

Answer

Let \(x_i\) be the \(i\)th node of \(\mathcal{Q}\). We define \(g : [a,b] \to \mathbb{R}\) to be \(\delta_{ik}\) at nodes \(x_i\), and on the restrictions \([x_i, x_{i+1}]\) between nodes be the linear interpolant of those values.

Then \(g\) is non-negative; however,

\[ \mathcal{Q}(g) = \sigma_k < 0. \]

Question

For a family \(\mathcal{Q}_n\) of quadrature methods (interpolatory or otherwise) one observes

\[ \sum_{k=0}^n \bigl|\sigma_k^{(n)}\bigr| \to \infty \qquad n \to \infty. \]

Why does the relative forward error, observed over the set of continuous functions, grow unboundedly as \(n \to \infty\)?

Answer

Let \(x_i^{(n)}\) be the \(i\)th node of \(\mathcal{Q}_n\). We define \(g^{(n)} : [a,b] \to \mathbb{R}\) to be \(\text{sign}(\sigma_k^{(n)})\) at nodes \(x_i^{(n)}\), and on the restrictions \([x_i^{(n)}, x_{i+1}^{(n)}]\) between nodes be the linear interpolant of those values.

Then

\[ \mathcal{Q}_n(g^{(n)}) = \sum_{k=0}^n \sigma_k^{(n)} g^{(n)}(x_k) = \sum_{k=0}^n | \sigma_k | \to \infty \qquad n \to \infty. \]

However, \(|\mathcal{I}(g^{(n)})| \leq (b-a)\). In the notation of the Foundation chapter on the forward error, the exact operator is \(f: C[a,b] \to \mathbb{R}, g \mapsto \mathcal{I}(g)\), while the \(n\)th quadrature approximation is \(\tilde{f}_n: C[a,b] \to \mathbb{R}, g \mapsto \mathcal{Q}_n(g)\).

Thus, we established that the largest relative forward error, observed over the set of continuous functions, diverges as \(n \to \infty\):

\[ \sup_{g \in C[a,b]} \frac{| f(g) - \tilde{f}_n(g) |}{|f(g)|} \geq \frac{| f(g^{(n)}) - \tilde{f}_n(g^{(n)}) |}{|f(g^{(n)})|} \to \infty, \]

using \(0 < f(g^{(n)}) \leq b - a\).

Question

Given a set of nodes \(a \leq x_0 < x_1 < \cdots < x_n \leq b\), demonstrate that there exist unique real numbers \(a_0, a_1, \ldots, a_n\) such that for every polynomial \(p\) in \(\mathcal{P}_n\) the following equality holds:

\[ \sum_{k = 0}^n a_k \, p(x_k) = \int_a^b p(x) \, dx. \]

Answer

The existence of the coefficients follows directly from the results about interpolatory quadrature formulas. For completeness, we summarise that the Lagrange basis polynomials

\[ L_k(x) = \prod_{l\ne k} \frac{x-x_l}{x_k-x_l} \]

constitute a basis for \(\mathcal{P}_n\). Any polynomial \(p\) in \(\mathcal{P}_n\) can be represented as a linear combination of these basis polynomials:

\[ p = \sum_{k = 0}^n p(x_k) \, L_k(x). \]

By integrating both sides over the interval \([a, b]\), we get:

\[ \int_a^b p(x) \, dx = \sum_{k = 0}^n p(x_k) \int_a^b L_k(x) \, dx. \]

The coefficients \(a_k\) can then be defined as:

\[ a_k := \int_a^b L_k(x) \, dx. \]

This definition ensures the existence of coefficients that satisfy the given relationship for all \(p\) in \(\mathcal{P}_n\).

To prove the uniqueness of the coefficients \(a_0, a_1, \ldots, a_n\), assume that another set of coefficients satisfies the question’s condition. Now the property of the Lagrange basis polynomials, \(L_k(x_j) = \delta_{kj}\) leads to:

\[ \int_a^b L_k(x) \, dx = \sum_{j = 0}^n a_j \, L_k(x_j) = a_k. \]

This implies that any other set of coefficients satisfying the condition must be identical to the \(a_k\)’s defined above, proving their uniqueness.

Question

a) Suppose the interpolatory quadrature formula \(\mathcal{Q}_1\) has the nodes \(1/4\) and \(2/3\) and the integration domain \([0,1]\). Compute the weights of \(\mathcal{Q}_1\).

b) Calculate \(\mathcal{Q}_1(- \ln)\) to approximate the integral

\[ - \int_0^1 \ln x \, dx. \]

Answer

a) The weights of \(\mathcal{Q}_1\) are given by

\[ \begin{eqnarray*} \sigma_k = \int_0^1 \prod_{m=0, \atop m \neq k}^n \frac{x-x_m}{x_k - x_m} \, dx. \end{eqnarray*} \]

Hence

\[\begin{split} \begin{eqnarray*} \sigma_0 &=& \int_0^1 \frac{x-2/3}{1/4 - 2/3} \, dx = 2/5,\\ \sigma_1 &=& \int_0^1 \frac{x-1/4}{2/3 - 1/4} \, dx = 3/5. \end{eqnarray*} \end{split}\]

b) In general, \(\mathcal{Q}_1(f) = \sigma_0 f(x_0) + \sigma_1 f(x_1)\). Here

\[ \begin{eqnarray*} \mathcal{Q}_1(- \ln) &=& - \frac{2}{5} \ln (1/4) - \frac{3}{5} \ln (2/3) \approx 0.79779680931285. \end{eqnarray*} \]

Note: The exact integral is

\[ \int_0^1 - \ln(x) dx = x - x \ln x\big|_{x \to 0}^1 = 1. \]

Question

Using \(\int_0^1 \frac{1}{\sqrt{x}} \, dx = 2\) and Simpson’s rule, find the approximate value of

\[ \begin{eqnarray*} \int_0^1 \frac{1}{\sqrt{\sin x}} \, dx. \qquad (*) \end{eqnarray*} \]

Answer

We cannot use Simpson’s rule directly to integrate \((*)\) numerically. The reason is that \(1/\sqrt{\sin x}\) diverges at the node \(x_0 = 0\). However, we can reformulate the problem as follows:

\[\begin{split} \begin{eqnarray*} \int_0^1 \frac{1}{\sqrt{\sin x}} \, dx &=& \int_0^1 \frac{1}{\sqrt{x}} \, dx + \int_0^1 \frac{1}{\sqrt{\sin x}} - \frac{1}{\sqrt{x}} \, dx\\ &=& 2 + \int_0^1 \frac{1}{\sqrt{\sin x}} - \frac{1}{\sqrt{x}} \, dx. \end{eqnarray*} \end{split}\]

The function \(\frac{1}{\sqrt{\sin x}} - \frac{1}{\sqrt{x}}\) is continuous and bounded on \([0,1]\) so that Simpson’s rule can be applied:

\[\begin{split} \begin{eqnarray*} \mathcal{Q}_2(\frac{1}{\sqrt{\sin x}} - \frac{1}{\sqrt{x}}) &=& (1-0) \Bigl( \frac{1}{6} \cdot 0 + \frac{2}{3} \cdot \Bigl(\frac{1}{\sqrt{\sin x}} - \frac{1}{\sqrt{x}} \Bigr) + \frac{1}{6} \cdot \Bigl( \frac{1}{\sqrt{\sin x}} - \frac{1}{\sqrt{x}} \Bigr) \Bigr)\\ &\approx& 0.03504029268746660. \end{eqnarray*} \end{split}\]

This gives the approximate value \(2.03504029268746660\). The exact integral is \(2.034805319207\ldots\)