Stable Quadrature Schemes

Stable Quadrature Schemes#

Composite Quadratures#

When dealing with functions \(f\) that may lack smoothness, it can be advantageous to divide the integration interval \([a,b]\) into \(m\) smaller, equally sized subintervals \(\{[x_{i-1},x_i]\}_{i=1}^m\). Each subinterval has a length of \(h:=(b-a)/m\). Within each subinterval, we apply a low-degree interpolatory quadrature formula.

For example, using the trapezium rule on \([x_{i-1},x_i]\), we find

\[\begin{split} \begin{align} \mathcal{C}_{1,m}(f) &= \frac{b-a}m\, \Bigl[ \Bigl( \frac{1}{2} f(x_0) + \frac{1}{2} f(x_1) \Bigr) + \Bigl(\frac{1}{2} f(x_1) + \frac{1}{2} f(x_2) \Bigr) + \cdots + \Bigl( \frac{1}{2} f(x_{m-1}) + \frac{1}{2} f(x_m) \Bigr) \Bigr] \notag\\ &= h \Bigl[ \frac{1}{2} f(x_0) + \sum_{i=1}^{m-1} f(x_i) + \frac{1}{2} f(x_m) \Bigr]. \end{align} \end{split}\]

This scheme involves \(m+1\) evaluations of \(f\). To use Simpson’s rule, we introduce an additional midpoint node \(x_{i-1/2}:=(x_{i-1}+x_i)/2\) in each subinterval. The composite Simpson’s rule is then given by:

\[\begin{split} \begin{align} \mathcal{C}_{2,m}(f) &= \frac{b-a}{m} \, \Bigl[ \frac{1}{6} f(x_0) + \frac{2}{3} f(x_{1/2}) + \frac{1}{6} f(x_1) + \frac{1}{6} f(x_1) + \frac{2}{3} f(x_{3/2}) + \frac{1}{6} f(x_2)\\ & \phantom{= \frac{b-a}{m} \, \Bigl[} + \cdots + \frac{1}{6} f(x_{m-1}) + \frac{2}{3} f(x_{m-1/2}) + \frac{1}{6} f(x_m)\Bigr]\\ &= \frac{h}{6} \Bigl[ f(x_0) + 4 \sum_{i=1}^{m} f(x_{i-1/2}) + 2 \sum_{i=1}^{m-1} f(x_i) + f(x_m) \Bigr]. \end{align} \end{split}\]

It involves \(2m+1\) evaluations of \(f\).

Error bounds for these methods can be obtained directly from those of the underlying methods. For \(\mathcal{C}_{1,m}\) we have, in the subinterval \([x_{i-1},x_i]\),

\[ \Bigl| \int_{x_{i-1}}^{x_i} f(x) dx - \tfrac{h}{2} (f(x_{i-1}) + f(x_i)) \Bigr| \le \frac{1}{2!}\,\max_{\xi\in[x_{i-1},x_i]}\,|f''(\xi)| \,\int_{x_{i-1}}^{x_i} |\omega_2(x)| dx = \frac{h^3}{12}\,\max_{\xi\in[x_{i-1},x_i]}\,|f''(\xi)|, \]

where \(\omega_2(x)=(x-x_{i-1})(x-x_i)\). Adding these up, we have in \([a,b]=[x_0,x_m]\),

\[ \bigl|\mathcal{I}(f)-\mathcal{C}_{1,m}(f)\bigr| \le \frac{mh^3}{12}\,\max_i\Bigl\{\max_{\xi\in[x_{i-1},x_i]}\,|f''(\xi)|\Bigr\} \leq (b-a)\,\frac{h^2}{12}\,\max_{\xi\in[a,b]}\,|f''(\xi)|. \]

Gaussian Quadratures#

The degree of exactness of an interpolatory quadrature \(\mathcal{Q}_n\) is at least \(n\). By considering \(\omega_{n+1}^2(x)=(x-x_0)^2\cdots(x-x_n)^2 \in \mathcal{P}_{2n+2}\), we find that the degree of exactness of any interpolatory quadrature cannot be equal to or larger than \(2n+2\). Gauss showed that by judicious choice of the nodes \(\{x_k\}\), one can indeed obtain a quadrature exact for all \(f \in \mathcal{P}_{2n+1}\). The resulting Gaussian quadrature rests on the following:

Theorem 9 (Gaussian quadrature)

The interpolatory quadrature formula \(\mathcal{Q}_n\) achieves a degree of exactness \(2n+1\) if and only if its nodes \(x_0, \ldots, x_n\) are such that the polynomial

\[ \omega_{n+1}(x) = \prod_{i=0}^n (x - x_i), \]

meets the condition \(\int_a^b \omega_{n+1}(x) \, p(x) \, dx = 0\) for all \(p \in \mathcal{P}_n\).

Thus, the nodes \(\{x_k\}_{k=0}^n\) are chosen so that the error polynomial \(\omega_{n+1} \in \mathcal{P}_{n+1}\) is in the orthogonal complement relative to the \((n+1)\)-dimensional space \(\mathcal{P}_n\) within the \((n+2)\)-dimensional space \(\mathcal{P}_{n+1}\) with respect to the inner product \(\langle p, q \rangle = \int_a^b p(x) \, q(x) \, dx\). A proof of the theorem is given in the optional materials section.

Definition 4

An interpolatory quadrature method that achieves a degree of exactness \(2n+1\) is termed a Gaussian quadrature method.

In short,

\[ \mathcal{Q}_n \text{ Gaussian} \qquad \Leftrightarrow \qquad \langle \omega_{n+1}, p \rangle = 0 \quad \forall \, p \in \mathcal{P}_n. \]

Once the nodes are determined, the coefficients can be computed as usual. Gaussian quadrature weights are always positive, unlike the coefficients \(\sigma_k\) of a general quadrature formula. This can be illustrated by considering

\[ L_k^2(x) = \prod_{i = 0 \atop i \neq k}^n \frac{(x-x_i)^2}{(x_k-x_j)^2} \in \mathcal{P}_{2n}. \]

With the degree of exactness of \(2n+1\) and \(L_k^2(x_j) = \delta_{jk}\), we find that:

\[ 0 < \mathcal{I}(L_k^2) = \mathcal{Q}_n(L_k^2) = \sigma_k. \]

In summary, Gaussian quadrature combines two crucial properties: stability and speed of convergence.

While we only scratch the surface of Gaussian quadrature, it’s close relationship to orthogonal polynomials makes it a fascinating topic for further study, e.g. see [Süli and Mayers, 2003] and [Plato, 2003] for more details.

Python skills#

Composite Newton-Cotes quadrature rules#

Scipy provides an implementation of the composite trapezoidal and Simpson’s rules, Numpy only has the former.

import numpy as np
from scipy import integrate

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

# Define the interval
a = 0  # Lower limit
b = 1  # Upper limit
n = 100  # Number of sample points

# Generate sample points
x = np.linspace(a, b, n)
y = my_function(x)

# Trapezoidal Rule
trapz_result_np = np.trapz(y, x)
trapz_result_scipy = integrate.trapz(y, x)

# Simpson's Rule
simps_result = integrate.simps(y, x)

# Display the results
print("Trapezoidal Rule (numpy):", trapz_result_np)
print("Trapezoidal Rule (scipy):", trapz_result_scipy)
print("Simpson's Rule:", simps_result)

Gaussian quadrature#

The code below compares exact integration with SciPy’s Gaussian quadrature implementation integrate.quad, which uses an adaptive composite formula. In addition, the Gaussian quadrature function integrate.fixed_quad with a specified number of quadrature points is evaluated.

import scipy.integrate as integrate

# Define a function to integrate and its exact integral
def my_function(x):
    p = 9
    return x**p

def integral_of_my_function(a,b):
    p = 9
    return (b**(p+1) - a**(p+1)) / (p+1)

# Define the interval
a = 1  # Lower limit
b = 5  # Upper limit

# Using exact integration formula
result_exact = integral_of_my_function(a,b)
print("Exact integration:", result_exact)

# Using scipy.integrate.quad, which uses an (adaptive) composite Gaussian quadrature
result_quad, _ = integrate.quad(my_function, a, b)
print("Gaussian Quadrature (scipy.integrate.quad):", result_quad)

# Using scipy.integrate.fixed_quad with different orders
for n in [2, 3, 4, 5]:  # Different number of quadrature points
    result_fixed_quad, _ = integrate.fixed_quad(my_function, a, b, n=n)
    print(f"Gaussian Quadrature of order {n} (fixed_quad):", result_fixed_quad)

Self-check questions#

Question

Show, for \(f \in C^3[a,b]\),

\[ \bigl|\mathcal{I}(f)-\mathcal{C}_{2,m}(f)\bigr| \le (b-a)\,\frac{h^3}{192}\,\max_{\xi\in[a,b]}\,|f'''(\xi)|. \]

Question

Let \(f(x) = \exp(-x^2)\). We wish to compare several methods to approximate

\[ \mathcal{I}(f) := \int_0^1 f(x) dx \]

where only three evaluations of \(f\) are allowed.

Writing your formulas clearly, compute approximations to \(\mathcal{I}(f)\) using the:

rectangle rule, i.e. composite quadrature with \(n=0\) and \(m=3\);
trapezium rule, i.e. composite quadrature with \(n=1\) and \(m=2\);
closed Newton-Cotes with \(n=2\).

Compare with the exact answer \(\mathcal{I}(f) = (\sqrt\pi/2)\,\textrm{erf}\,1=0.746824132812\cdots\).

Question

Let \([a,b] = [-1,1]\).

Given the nodes \(\{t_k\}=\bigl\{1/2,\bigl(1\pm\sqrt{3/5}\bigr)/2\bigr\}\), compute the coefficients \(\{\sigma_k\}\) of the interpolatory quadrature formula.
Use your result to compute the interpolatory quadrature \(\mathcal{Q}_2'(f)\) for \(f(x)=\exp(-x^2)\) and compare with the results of the methods of the previous question.

Question

Find nodes \(\{t_0,t_1\}\subset[0,1]\) and coefficients \(\{\sigma_0,\sigma_1\}\) such that the resulting quadrature formula has degree of exactness \(3\), i.e.

\[ \sigma_0 f(t_0) + \sigma_1 f(t_1) = \int_0^1 f(t) \>\mathrm{d}t, \]

for every \(f \in \mathcal{P}_3\). Hint: consider polynomials in \(t-\frac{1}{2}\) instead of \(t\).

Question

Similarly, find nodes \(\{t_0,t_1,t_2\}\subset[0,1]\) and coefficients \(\{\sigma_0,\sigma_1,\sigma_2\}\) such that

\[ \sigma_0 f(t_0) + \sigma_1 f(t_1) + \sigma_2 f(t_2) = \int_0^1 f(t) dt \]

for every \(f \in \mathcal{P}_5\). Hint: symmetry!

Question

Decompose the interval \([a,b]\) into \(m\) subintervals of equal length \(h := (b-a)/m\). Let \(\mathcal{C}_{h,1}\) be the composite trapezium rule defined on these subintervals and let \(\mathcal{C}_{h,2}\) be the composite Simpson’s rule defined on these subintervals (thus for the calculation of \(\mathcal{C}_{h,2}(f)\) the function \(f\) is evaluated with a step-size of \(h/2\)). Show that

\[ \mathcal{C}_{h,2}(f) = \frac{4}{3} \mathcal{C}_{h/2,1}(f) - \frac{1}{3} \mathcal{C}_{h,1}(f). \]

Here \(\mathcal{C}_{h/2,1}\) is the trapezium rule with subintervals of the length \(h/2\).

Stable Quadrature Schemes

Contents

Stable Quadrature Schemes#

Composite Quadratures#

Gaussian Quadratures#

Python skills#

Composite Newton-Cotes quadrature rules#

Gaussian quadrature#

Self-check questions#

Optional material#