from numpy import *
from scipy import linalg
from scipy.linalg import norm
from matplotlib.pyplot import *
from prettytable import PrettyTable
from timeit import default_timer as timer
import sys
sys.path.append('fncbook/')
import fncbook as FNC

# This (optional) block is for improving the display of plots.
# from IPython.display import set_matplotlib_formats
# set_matplotlib_formats("svg","pdf")
# %config InlineBackend.figure_format = 'svg'
rcParams["figure.figsize"] = [7, 4]
rcParams["lines.linewidth"] = 2
rcParams["lines.markersize"] = 4
rcParams['animation.html'] = "jshtml"  # or try "html5"

9.6.Spectrally accurate integration¶

In Numerical integration we derived methods of order 2, 4, and higher for numerical integration. They all use a formula

\int_{-1}^1 f(x)\, dx \approx \sum_{k=0}^n w_k f(t_k)

(9.6.1)

for a collection of nodes $t_0,\ldots,t_n$ in $[-1,1]$ and weights $w_0,\ldots,w_n$ . (Throughout this section we use $[-1,1]$ as the domain of the integral; for a general interval $[a,b]$ , see Exercise 9.6.4.) The nodes and weights are independent of the integrand $f(x)$ and determine the implementation and properties of the formula.

The process for deriving a specific method was to interpolate the integrand, then integrate the interpolant. Piecewise linear interpolation at equally spaced nodes, for instance, produces the trapezoid formula. When the integrand is approximated by a spectrally accurate global function, the integration formulas are also spectrally accurate.

9.6.1Periodic functions¶

For a function periodic on $[-1,1]$ , the most natural interpolant is the trigonometric polynomial (9.5.4). However, from (9.5.3) one finds that

\int_{-1}^1 \sum_{k=-n}^n y_k\tau_k(x)\, dx = \sum_{k=-n}^n y_k\left[ \int_{-1}^1 \tau_k(x)\, dx\right] = \frac{2}{2n+1} \sum_{k=-n}^n y_k.

(9.6.2)

In Exercise 9.6.1 you are asked to verify that this result is identical to the value of the trapezoid formula on $2n+1$ nodes.

Example 9.6.1 (Perimeter of an ellipse)

We use the trapezoidal integration formula to compute the perimeter of an ellipse with semi-axes 1 and 1/2. Parameterizing the ellipse as $x=\cos \pi t$ , $y=\frac{1}{2}\sin \pi t$ leads to the arc-length integral

\int_{-1}^1 \pi\sqrt{ \cos^2(\pi t) + \tfrac{1}{4}\sin^2(\pi t)}\,dt.

(9.6.3)

f = lambda t: pi * sqrt(cos(pi * t) ** 2 + sin(pi * t) ** 2 / 4)
N = arange(4, 48, 6)
perim = zeros(N.size)
for k in range(N.size):
    h = 2 / N[k]
    t = h * arange(N[k]) - 1
    perim[k] = h * sum(f(t))
err = abs(perim - perim[-1])    # use the last value as reference
results = PrettyTable()
results.add_column("N", N)
results.add_column("perimeter", perim)
results.add_column("error", err)
results

The approximations gain about one digit of accuracy for each constant increment of $n$ , which is consistent with spectral convergence.

9.6.2Clenshaw–Curtis integration¶

Suppose $f$ is smooth but not periodic. If we use a global polynomial interpolating $f$ at the Chebyshev second-kind points from (9.3.2),

%:label: chebextremerepeat t_k = - \cos\left(\frac{k \pi}{n}\right), \qquad k=0,\ldots,n,

(9.6.4)

and integrate the resulting polynomial interpolant, the method should have spectral accuracy for a smooth integrand. The resulting algorithm is known as Clenshaw–Curtis integration.

Having specified the nodes in (9.6.1), all that remains is to find the weights. The Lagrange form of the interpolating polynomial is

p(x) = \sum_{k=0}^{n} f(t_k) \ell_k(x).

From this,

\begin{align*} I = \int_{-1}^1 f(x)\, dx \approx \int_{-1}^1 p(x) \, d x &= \int_{-1}^1 \sum_{k=0}^n f(t_k) \ell_k(x) \, d x\\ &= \sum_{k=0}^n f(t_k) \int_{-1}^1 \ell_k(x) \, d x = \sum_{k=0}^n w_k f(t_k), \end{align*}

where $w_k = \int_{-1}^1 \ell_k(x)\,dx.$ For even values of $n$ , the result is

w_k = \begin{cases} \dfrac{1}{n^2-1}, & \text{$k=0$ or $k=n$},\\[3mm] \dfrac{4}{n} \displaystyle \sum_{j=0}^{n/2} \frac{\cos ( 2 \pi j k / n)}{\gamma_j (1-4 j^2) }, & k=1,\dots,n-1, \end{cases} \quad \gamma_j = \begin{cases} 2, & j=0 \text{ or } n/2,\\ 1, & j=1,2,\dots,n/2-1. \end{cases}

(9.6.5)

There are different formulas for odd values of $n$ . Note that the weights also depend on $n$ ; e.g. $w_2$ for $n=4$ is not the same as $w_2$ for $n=10$ . Also note that the interpolant itself never needs to be computed.

Function 9.6.1 performs Clenshaw–Curtis integration for even values of $n$ .^[1]

Algorithm 9.6.1 (ccint)

ccint.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def ccint(f, n):
    """
    ccint(f, n)

    Perform Clenshaw-Curtis integration for the function f on n+1 nodes in [-1,1]. Return
    integral and a vector of the nodes used. Note: n must be even.
    """
    # Find Chebyshev extreme nodes.
    theta = np.linspace(0, np.pi, n + 1)
    x = -np.cos(theta)

    # Compute the C-C weights.
    c = np.zeros(n + 1)
    c[[0, n]] = 1 / (n**2 - 1)
    v = np.ones(n - 1)
    for k in range(1, int(n / 2)):
        v -= 2 * np.cos(2 * k * theta[1:-1]) / (4 * k**2 - 1)
    v -= np.cos(n * theta[1:-1]) / (n**2 - 1)
    c[1:-1] = 2 * v / n

    # Evaluate integrand and integral.
    I = np.dot(c, f(x))  # use vector inner product
    return I, x

9.6.3Gauss–Legendre integration¶

Let us reconsider the generic numerical integration formula (9.6.1),

\int_{-1}^1 f(x)\, dx \approx \sum_{k=1}^n w_k f(t_k) = Q_{n}[f],

(9.6.6)

where $Q_n[f]$ stands for the application of the formula to function $f$ . (We start the sum from $k=1$ instead of $k=0$ for notational convenience in what follows.)

The interpolation approach spurred us to use Chebyshev nodes. But it’s not clear that these are ideal nodes for the specific application of finding an integral. Instead, we can define formula as the integral of a polynomial interpolant, but with the weights and nodes chosen to satisfy an optimality criterion. As usual, we denote the set of all polynomials of degree at most $m$ by $\mathcal{P}_m$ .

Since there are $n$ nodes and $n$ weights available to choose, it seems plausible to expect $m=2n-1$ , and this intuition turns out to be correct. Hence, the goal is now to find nodes $t_k$ and weights $w_k$ such that

\int_{-1}^1 p(x)\,dx = Q_{n}[p] = \sum_{k=1}^n w_k p(t_k), \qquad p \in \mathcal{P}_{2n-1}.

(9.6.8)

If these conditions are satisfied, the resulting method is called Gauss–Legendre integration or simply Gaussian integration. Because the integration formula is linear, i.e., $Q_n[\alpha p + q] = \alpha Q_n[p] + Q_n[q]$ , it is sufficient to show that $Q_n$ gets the exact value for the monomials $1,x,x^2,\ldots,x^{2n-1}.$

Generalizing the process above to general $n$ would be daunting, as the conditions on the nodes and weights are nonlinear. Fortunately, a more elegant approach is possible.

Proof

Choose an arbitrary $p\in\mathcal{P}_{2n-1}$ , and let $\hat{p}_n(x)$ be the lowest-degree interpolating polynomial for $p$ using the as-yet unknown nodes $t_1,\dots,t_n$ . By definition,

Q_n[p] = \int_{-1}^1 \hat{p}_n(x)\, dx.

(9.6.11)

Since $\hat{p}_n(x)$ has degree at most $n-1$ , it is exactly equal to $p$ if $p\in\mathcal{P}_{n-1}$ , and (9.6.8) is trivially satisfied. Otherwise, the error formula (9.1.6) implies

p(x) - \hat{p}_n(x) = \frac{p^{(n)}(\xi(x))}{n!} \Phi(x),

(9.6.12)

where $\Phi$ is the error indicator function $\prod_k (t-t_k).$ Trivially, the left-hand side is a polynomial in $\mathcal{P}_{2n-1}$ of degree at least $n$ , so the right-hand side must be too. Thus, we can write

p(x) - \hat{p}_n(x) = \Psi(x) \Phi(x),

(9.6.13)

where $\Psi(x)\in \mathcal{P}_{n-1}$ is unknown. The optimality requirement (9.6.8) becomes

0 = \int_{-1}^1 p(x)\,dx - Q_{n}[p] = \int_{-1}^1 \bigl[p(x) - \hat{p}_n(x)\bigr]\,dx = \int_{-1}^1 \Psi(x) \Phi(x) \, dx.

(9.6.14)

Given that $\Psi(x)\in \mathcal{P}_{n-1}$ , we can ensure that this condition is satisfied if

\int_{-1}^1 q(x)\Phi(x) \,dx = 0

(9.6.15)

for all $q \in {\mathcal{P}}_{n-1}$ . Hence, satisfaction of (9.6.15) implies satisfaction of (9.6.8). But by the orthogonality property of Legendre polynomials, satisfaction of (9.6.15) is guaranteed if $\Phi(x)=cP_n(x)$ for a constant $c$ . Thus, $\Phi$ and $P_n$ have the same roots.

From Theorem 9.4.3 we know that the roots of $P_n$ are distinct and all within $(-1,1)$ . (Indeed, it would be strange to have the integral of a function depend on some of its values outside the integration interval!) While there is no explicit formula for the roots, there are fast algorithms to compute them and the integration weights on demand. Function 9.6.2 uses one of the oldest methods, practical up to $n=100$ or so.

Algorithm 9.6.2 (glint)

glint.py

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def glint(f, n):
    """
    glint(f, n)

    Perform Gauss-Legendre integration for the function f on n nodes in (-1,1). Return
    integral and a vector of the nodes used.
    """
    # Nodes and weights are found via a tridiagonal eigenvalue problem.
    beta = 0.5 / np.sqrt(1 - (2.0 * np.arange(1, n)) ** (-2))
    T = np.diag(beta, 1) + np.diag(beta, -1)
    ev, V = eig(T)
    ev = np.real_if_close(ev)
    p = np.argsort(ev)
    x = ev[p]  # nodes
    c = 2 * V[0, p] ** 2  # weights

    # Evaluate the integrand and compute the integral.
    I = np.dot(c, f(x))  # vector inner product
    return I, x

9.6.4Convergence¶

Both Clenshaw–Curtis and Gauss–Legendre integration are spectrally accurate. The Clenshaw–Curtis method on $n+1$ points has degree $n$ , whereas the Gauss–Legendre method with $n$ points has degree ${2n-1}$ . For this reason, it is possible for Gauss–Legendre to converge at a rate that is “twice as fast,” i.e., with roughly the square of the error of Clenshaw–Curtis. But the full story is not simple.

Example 9.6.3 (Comparing spectral integration methods)

First consider the integral

\int_{-1}^1 \frac{1}{1+4x^2} \, dx = \arctan(2).

(9.6.16)

f = lambda x: 1 / (1 + 4 * x**2)
exact = arctan(2)

We compare the two spectral integration methods for a range of $n$ values.

N = range(8, 100, 4)
errCC = zeros(len(N))
errGL = zeros(len(N))
for k, n in enumerate(N):
    errCC[k] = exact - FNC.ccint(f, n)[0]
    errGL[k] = exact - FNC.glint(f, n)[0]

semilogy(N, abs(errCC), "-o", label="Clenshaw–Curtis")
semilogy(N, abs(errGL), "-o", label="Gauss–Legendre")
xlabel("number of nodes"),  ylabel("error"),  ylim(1e-16, 0.01)
legend(),  title("Spectral integration");

(The missing dots are where the error is exactly zero.) Gauss–Legendre does converge faster here, but at something less than twice the rate.

Now we try a more sharply peaked integrand:

\int_{-1}^1 \frac{1}{1+16x^2} \, dx = \frac{1}{2}\arctan(4).

(9.6.17)

f = lambda x: 1 / (1 + 16 * x**2)
exact = atan(4) / 2

N = range(8, 100, 4)
errCC = zeros(len(N))
errGL = zeros(len(N))
for k, n in enumerate(N):
    errCC[k] = exact - FNC.ccint(f, n)[0]
    errGL[k] = exact - FNC.glint(f, n)[0]

semilogy(N, abs(errCC), "-o", label="Clenshaw–Curtis")
semilogy(N, abs(errGL), "-o", label="Gauss–Legendre")
xlabel("number of nodes"),  ylabel("error"),  ylim(1e-16, 0.1)
legend(),  title("Spectral integration");

The two are very close until about $n=40$ , when the Clenshaw–Curtis method slows down.

Now let’s compare the spectral performance to that of our earlier adaptive method in intadapt. We will specify varying error tolerances and record the error as well as the total number of evaluations of $f$ .

loglog(N, abs(errCC), "-o", label="ccint")
loglog(N, abs(errGL), "-o", label="glint")

tol_ = 1 / 10 ** arange(2, 15)
n = zeros(tol_.size)
errAdapt = zeros(tol_.size)
for k, tol in enumerate(tol_):
    Q, t = FNC.intadapt(f, -1, 1, tol)
    errAdapt[k] = exact - Q
    n[k] = t.size

loglog(n, abs(errAdapt), "-o", label="intadapt")
loglog(n, 1 / (n**4), "--", label="4th order")
xlabel("number of nodes"),  ylabel("error"),  ylim(1e-16, 1)
legend(),  title("Spectral vs 4th order");

At the core of intadapt is a fourth-order formula, and the results track that rate closely. For all but the most relaxed error tolerances, both spectral methods are far more efficient than the low-order counterpart. For other integrands, particularly those that vary nonuniformly across the interval, the adaptive method might be more competitive.

The difference in convergence between Clenshaw–Curtis and Gauss–Legendre is dwarfed by the difference between spectral and algebraic convergence. It is possible, though, to encounter integrands for which adaptivity is critical. Choosing a method is highly problem-dependent, but a rule of thumb is that for large error tolerances, an adaptive low-order method is likely to be a good choice, while for high accuracy, the spectral methods often dominate.

9.6.5Exercises¶

Exercise 9.6.3

⌨ For each integral, compute approximations using Function 9.6.1 and Function 9.6.2 with $n=4,6,8,\ldots,40$ . Plot the errors of both methods together as functions of $n$ on a semi-log scale.

(a) $\displaystyle\int_{-1}^1 e^{-4x}\, dx = \sinh(4)/2$

(b) $\displaystyle\int_{-1}^1 e^{-9x^2} = \sqrt{\pi}\, \operatorname{erf}(3)/3$

(c) $\displaystyle\int_{-1}^1 \operatorname{sech}(x) \, dx = 2 \tan^{-1} [ \sinh (1) ]$

(d) $\displaystyle\int_{-1}^1 \frac{1}{1+9x^2}\, dx = \frac{2}{3} \tan^{-1}(3)$

Footnotes¶

This function is modeled after the function clencurt.m of Trefethen (2000).
↩

References¶

Trefethen, L. N. (2000). Spectral Methods in MATLAB. SIAM.