Spectrally accurate integration - Fundamentals of Numerical Computation

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 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 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)

Julia

MATLAB

Python

Example 9.6.1

f(t) = π * sqrt( cospi(t)^2 + sinpi(t)^2 / 4 );
n = 4:4:48
perim = zeros(size(n))
for (k, n) in enumerate(n)
    h = 2 / n
    t = @. h * (0:n-1) - 1
    perim[k] = h * sum(f.(t))
end
err = @. abs(perim - perim[end])    # use last value as "exact"
@ptconf formatters=ft_printf(["%d", "%.15f", "%.2e"], 1:3)
@pt :header=["n", "perimeter", "error estimate"] [n perim err][1:end-1, :]

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

Example 9.6.1

f = @(t) pi * sqrt( cos(pi*t).^2 + sin(pi*t).^2 / 4 );
N = (4:4:48)';
C = zeros(size(N));
for k = 1:length(N)
    h = 2 / N(k);
    t = h * (0:N(k)-1);
    C(k) = h * sum(f(t));
end
format long
disp(table(N, C, variableNames=["number of nodes", "perimeter"]))

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

Example 9.6.1

f = lambda t: pi * sqrt(cos(pi * t) ** 2 + sin(pi * t) ** 2 / 4)
N = arange(4, 48, 6)
results = PrettyTable(["N", "perimeter estimate"])
for k in range(N.size):
    h = 2 / N[k]
    t = h * arange(N[k]) - 1
    results.add_row([N[k], h * sum(f(t))])

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(x_k) \ell_k(x).

(9.6.5)

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(x_k) \ell_k(x) \, d x\\ &= \sum_{k=0}^n f(x_k) \int_{-1}^1 \ell_k(x) \, d x = \sum_{k=0}^n w_k f(x_k), \end{align*}

(9.6.6)

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.7)

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)

Julia

MATLAB

Python

Clenshaw–Curtis integration

Clenshaw-Curtis integration

ccint.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
function [I,x] = ccint(f,n)
% CCINT  Clenshaw-Curtis numerical integration.
% Input:
%   f     integrand (function)
%   n     one less than the number of nodes (even integer)
% Output:
%   I     estimate of integral(f,-1,1)
%   x     evaluation nodes of f (vector)

% Find Chebyshev extreme nodes.
theta = pi*(0:n)'/n;
x = -cos(theta);

% Compute the C-C weights.
c = zeros(1,n+1); 
c([1 n+1]) = 1/(n^2-1); 
theta = theta(2:n); 
v = ones(n-1,1);
for k = 1:n/2-1
  v = v - 2*cos(2*k*theta)/(4*k^2-1);
end
v = v - cos(n*theta)/(n^2-1);
c(2:n) = 2*v/n;

% Evaluate integrand and integral.
I = c*f(x);   % use vector inner product

Clenshaw–Curtis integration

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.8)

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.10)

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 9.6.1

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.13)

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.10) is trivially satisfied. Otherwise, the error formula (9.1.8) implies

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

(9.6.14)

where Φ 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.15)

where $\Psi(x)\in \mathcal{P}_{n-1}$ is unknown. The optimality requirement (9.6.10) 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.16)

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.17)

for all $q \in {\mathcal{P}}_{n-1}$ . Hence satisfaction of (9.6.17) implies satisfaction of (9.6.10). But by the orthogonality property of Legendre polynomials, satisfaction of (9.6.17) is guaranteed if $\Phi(x)=cP_n(x)$ for a constant $c$ . Thus Φ 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)

Julia

MATLAB

Python

Gauss–Legendre integration

Gauss-Legendre integration

glint.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
function [I,x] = glint(f,n)
% GLINT  Gauss-Legendre numerical integration.
% Input:
%   f     integrand (function)
%   n     number of nodes (integer)
% Output:
%   I     estimate of integral(f,-1,1)
%   x     evaluation nodes of f (vector)

% Nodes and weights are found via a tridiagonal eigenvalue problem.
beta = 0.5./sqrt(1-(2*(1:n-1)).^(-2));
T = diag(beta,1) + diag(beta,-1);
[V,D] = eig(T);
x = diag(D); [x,idx] = sort(x);    % nodes
c = 2*V(1,idx).^2;                 % weights

% Evaluate the integrand and compute the integral.
I = c*f(x);      % vector inner product

Gauss–Legendre integration

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)

Julia

MATLAB

Python

Example 9.6.3

First consider the integral

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

(2)

f(x)= 1 / (1 + 4x^2);
exact = atan(2);

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

using Plots
n = 8:4:96
err = zeros(length(n), 2)
for (k, n) in enumerate(n)
  err[k, 1] = abs(exact - FNC.ccint(f, n)[1])
  err[k, 2] = abs(exact - FNC.glint(f, n)[1])
end

err[iszero.(err)] .= NaN    # remove from log-scale plot
plot(n, err, m=:o, label=["CC" "GL"],
    xaxis=("number of nodes"),  yaxis=(:log10, "error", [1e-16, 1]), 
    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).

(3)

f(x) = 1 / (1 + 16x^2);
exact = atan(4) / 2;

n = 8:4:96
err = zeros(length(n), 2)
for (k,n) in enumerate(n)
  err[k, 1] = abs(exact - FNC.ccint(f, n)[1])
  err[k, 2] = abs(exact - FNC.glint(f, n)[1])
end

err[iszero.(err)] .= NaN    # remove from log-scale plot
plot(n, err, m=:o, label=["CC" "GL"],
    xaxis=("number of nodes"),  yaxis=(:log10, "error", [1e-16, 1]), 
    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$ .

tol = 10 .^(-2.0:-2:-14)
n = zeros(size(tol))  
errAdapt = zeros(size(tol))
for (k, tol) in enumerate(tol)
  Q, t = FNC.intadapt(f, -1, 1, tol)
  errAdapt[k] = abs(exact - Q)
  n[k] = length(t)
end

errAdapt[iszero.(errAdapt)] .= NaN
plot!(n, errAdapt, m=:o, label="intadapt")
plot!(n, n.^(-4), l=:dash, label="4th order",
        xaxis=(:log10),  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.

Example 9.6.3

First consider the integral

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

(1)

f = @(x) 1 ./ (1 + 4*x.^2);
exact = atan(2);

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

n = (8:4:96)';
errCC = zeros(size(n));
errGL = zeros(size(n));
for k = 1:length(n)
  errCC(k) = exact - ccint(f, n(k));
  errGL(k) = exact - glint(f, n(k));
end
clf,  semilogy(n, abs([errCC errGL]), 'o-')
xlabel('number of nodes'),  ylabel('error')
title('Spectral integration')   
legend('Clenshaw–Curtis', 'Gauss–Legendre')

(The missing points 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).

(2)

f = @(x) 1 ./ (1 + 16*x.^2);
exact = atan(4) / 2;

n = (8:4:96)';
errCC = zeros(size(n));
errGL = zeros(size(n));
for k = 1:length(n)
  errCC(k) = exact - ccint(f, n(k));
  errGL(k) = exact - glint(f, n(k));
end
clf,  semilogy(n, abs([errCC errGL]), 'o-')
xlabel('number of nodes'),  ylabel('error')
title('Spectral integration')   
legend('Clenshaw–Curtis', 'Gauss–Legendre')

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

tol = 10 .^ (-2:-2:-14)';
n = zeros(size(tol));  
errAdapt = zeros(size(tol));
for k = 1:length(n)
  [Q, t] = intadapt(f, -1, 1, tol(k));
  errAdapt(k) = exact - Q;
  n(k) = length(t);
end
hold on;  semilogy(n, abs(errAdapt), 'o-')
plot(n, n.^(-4), 'k--')        % 4th order error
set(gca, 'xscale', 'log')     
legend('ccint', 'glint', 'intadapt', '4th order')  
title(('Spectral vs 4th order'));

Example 9.6.3

First consider the integral

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

(1)

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).

(2)

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.

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");

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¶

✍ Suppose $f$ is periodic on $[-1,1]$ . Show that the result of applying the trapezoid formula on $2n+1$ points is identical to (9.6.2).
✍ For each integral, use Gauss–Legendre integration with $n=2$ to write out the terms $w_1f(t_1)$ and $w_2f(t_2)$ explicitly.
(a) $\displaystyle\int_{-1}^1 e^{-x}\, dx = 2 \sinh(1) \qquad$ (b) $\displaystyle\int_{-1}^1 e^{-x^2} \qquad$ (c) $\displaystyle\int_{-1}^1 (2+x)^{-1}$

⌨ 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)$

(a) ✍ (See also Exercise 9.3.5.) Using the change of variable
$z = \phi(x) = a + (b-a)\frac{(x+1)}{2},$
(9.6.18)
show that
$\int_{a}^b f(z) \, d z= \frac{b-a}{2} \int_{-1}^{1} f( \phi(x) ) \, d x .$
(9.6.19)
(b) ⌨ Rewrite Function 9.6.1 and Function 9.6.2 to accept additional inputs for $a$ and $b$ and compute integrals over $[a,b]$ .
(c) ⌨ Repeat the steps of Exercise 2 for the integral
$\int_{\pi/2}^{\pi} x^2 \sin 8x \, d x = -\frac{3 \pi^2}{32}.$
(9.6.20)
⌨ A particle moves in a circular path with angular velocity given by $\omega(\theta)=\sin(\exp(\sin \theta))$ . The time it takes to complete one full orbit is
$T = \int_0^T dt = \int_{0}^{2\pi} \frac{d\theta}{d\theta / dt } = \int_{0}^{2\pi} \frac{d\theta}{\omega(\theta) }.$
(9.6.21)
Use Function 5.6.1 to find the period of the motion. Show results for different values of $n$ to establish convergence to at least 12 digits.
✍ Prove the claim about linearity of the Gauss–Legendre integration formula alluded to in the derivation of Theorem 9.6.1. Namely, show that condition (9.6.10) is true if and only if
$\int_{-1}^1 x^j\,dx = \sum_{k=1}^n w_k x_k^j$
(9.6.22)
for all $j=0,\ldots,2n-1$ .

Footnotes¶

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

References¶

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

Preface

Trigonometric interpolation

Preface

Improper integrals