. 94
( 136 .)


2. x ∈ [0, ∞]: The transformation of the integral is (Boyd, 1987b, and Chapter 17), using
the same map as yields the T L rational Chebyshev basis:
x = L cot2 (t/2), L is a user-choosable constant (19.54)
∞ π
f (L cot2 (t/2)) 2L (19.55)
I = f (x)dx = dt
[1 ’ cos(t)]2
0 0

The quadrature is
≡ wj f (L cot2 (tj /2)) (19.56)
≡ πj/(N + 1), (19.57)
tj j = 1, 2, . . . , N
[1 ’ cos(mπ)]
sin(tj ) 2
≡ 2L (19.58)
wj sin(mtj )
[1 ’ cos(tj )]2 N + 1 m

3. x ∈ [’∞, ∞]: The transformation of the integral is (Boyd, 1987b, and Chapter 17),
using the same map as yields the T B rational Chebyshev basis:

L is a user-choosable constant (19.59)
x = L cot(t),
∞ π
I = f (x)dx = f (L cot(t)) 2L dt
sin2 (t)
0 0

The quadrature is
≡ (19.61)
IN wj f (L cot(tj ))
≡ πj/(N + 1), (19.62)
tj j = 0, 1, . . . , N + 1
Lπ/[sin2 (tj )(N + 1)], 0 < j < (N + 1)
≡ (19.63)
Lπ/[sin2 (tj )(2N + 2)], j = 0 and j = N + 1

Because the transformed integrals have integrands which are antisymmetric with respect
to t = 0 for the ¬nite and semi-in¬nite intervals, Boyd expanded the product of f with the
metric factor in terms of sine cardinal functions. Thus, the endpoint are omitted for these
two cases so that each uses only the N interior points of the (N + 2)-point Lobatto grid.
Alternatively, one may expand only f as a Chebyshev or Fourier series and then evalu-
ate integrals of the products of the cardinal functions with the metric factor, which is sin(t)
for the ¬nite interval. This requires two additional evaluations of f , but raises the degree
of f for which the approximation is exact only by one. The quadrature weights given by
Fraser and Wilson (1966) as their Eq. (4.3).
Clenshaw and Curtis (1960) and Gentleman(1972a, b, c) prefer to compute the Cheby-
shev series for f (x) ¬rst by a Fast Fourier transform and then integrate term-by-term
through a recurrence formula. This is an O(N log2 (N )) procedure whereas the cost of eval-
uating the weights through Boyd™s method is O(N 2 ). Since these weights need only be
computed once for a given N , however, we recommend Boyd™s procedure.

19.8.4 Integration of Periodic Functions and the Trapezoidal Rule
Every student of numerical analysis learns the trapezoidal rule, and learns that the error is
O(1/N 2 ) where N is the number of grid points. Indeed, this estimate of only second order
accuracy is correct for non-periodic functions. However, when the integrand f (x) is periodic
and in¬nitely differentiable for all real x and the integration interval coincides in length with
the spatial period P , the trapezoidal rule becomes spectrally accurate.

trapezoidal rule approximation to the integral of a function f (x) on the interval x ∈ [’P/2, P/2]
± 
 f (’P/2) + f (P/2) 
N ’1
f (x)dx ≡ IN = h (19.64)
+ f (’P/2 + hk)
 
’P/2 k=’(N ’1)

h ≡ P/N is the grid spacing. If f (x) has the Fourier series
∞ ∞
1 2π 2π
f (x) ≡ ±0 + x ∈ [’P/2, P/2], (19.65)
±n cos n x + βn sin n x ,
2 P P
n=1 n=1

then the error in the trapezoidal rule is given without approximation as
EN ≡ f (x)dx ’ IN = ’P {±N + ±2N + ±3N + . . . } (19.66)

If f (x) is non-periodic, then the Fourier coef¬cients decrease as O(1/N 2 ) and the trapezoidal
rule is only second-order accurate, as taught in elementary numerical analysis texts.
If f (x) is periodic with period P and analytic on the interval x ∈ [’P/2, P/2], then the Fourier
coef¬cients are exponentially small in N and the trapezoidal rule is spectrally accurate.

Proof: The error series follows trivially from Theorem 19 because the integral is also the
constant in the Fourier series. The rest of the theorem follows from the Fourier convergence
theory for non-periodic and periodic functions described in Chap. 2. Q. E. D.
It is important to note that since the trapezoidal rule uses an evenly spaced grid, the
Clenshaw-Curtis adaptivity still applies.

19.8.5 In¬nite Intervals and the Trapezoidal Rule
If f (x) is a function that decays exponentially fast with |x|, then integrals on the in¬nite
interval can be performed with spectral accuracy by applying the trapezoidal rule on a
large but ¬nite interval through “domain truncation” (Sec 17.2). The approximation is the
same as Eq.(19.64),

N ’1

f (’P/2) + f (P/2)
f (x)dx ≈ IN ≡ h (19.67)
+ f (’P/2 + hk)
’∞ k=1


|I ’ IN | ∼ max (|f (x)|) + O(±N ) (19.68)

where ±N is again the N -th Fourier coef¬cient when f (x) is expanded as a Fourier series
on the interval x ∈ [’P/2, P/2]. Note that because f (x) is exponentially small at the ends
of the (large!) interval x ∈ [’P/2, P/2], it is quite unnecessary for f (x) to be periodic, and
the discontinuity in the Fourier series at the ends of the interval will be O(f (P/2)) and
therefore exponentially small if we have chosen a suf¬ciently large P .
Because of its simplicity, the trapezoidal rule is usually preferable to the T B-integration
formula (19.61). However, the latter can be sometimes applied to integrands that decay
only as an inverse power of x whereas the trapezoidal rule is always restricted to expo-
nentially decaying f . Its accuracy is also more sensitive to the choice of P than is the
T B-integration rule to the choice of map parameter L (Boyd, 1982a). However, because
numerical integration is fairly cheap task, this sensitivity of the trapezoidal rule is rarely a

19.8.6 Singular Integrands
Integrals with integrable singularities can be computed to spectral accuracy by making a
change of coordinates. For example, suppose that the function is singular at both ends of
the interval x ∈ [’1, 1]. Make the change of variable

x = tanh(z/L) (19.69)

Error in Tanh Rule
(1-x 2)-3/4 dx

sum of (1-x 2)-3/4


sum of sech(zk)

0 100 200 300 400

Figure 19.4: An illustration of the exponential convergence of the tanh-rule where N is the
number of grid points. The integrand (1 ’ x2 )’3/4 is singular at both endpoints. After the
change of coordinates to x = tanh(x/L) where L is arbitrarily chosen to be one here, the
integral is converted to a nonsingular integral on z ∈ [’∞, ∞]. The “tanh rule” is simply
the trapezoidal rule, applied on a large but ¬nite interval in z. The tanh rule is sensitive to
roundoff. When the integrand is expressed in terms of sech(z), one can obtain better than
ten decimal place accuracy (thick, solid curve). However, when the integrand is computed
in its original form (1’x2 )’3/4 , the tanh rule gives only three correct digits (dashed curve).
The dif¬culty with (1’x2 )’3/4 is that near the endpoints where x2 ≈ 1, this is the fractional
power of the small difference of two much larger terms, generating large roundoff error.
(Note that the tanh rule has an exponentially high density of points, i. e., lots of points, in
these endpoint regions of strong cancellation.)

where L is a user-choosable scaling factor. Noting that dx/dz = (d/dz)tanh(z/L) = (1/L)sech2 (z/L),
it follows that

f (tanh(z/L)) sech2 (z/L)dz
I≡ (19.70)
f (x) dx =
’1 ’∞

If f (x) blows up more slowly than a ¬rst order pole as x ’ 1, that is, f (x) grows no faster
than (1 ’ x)’± for some ± < 1 and similarly at the other endpoint, then (i) the original
integral in x is ¬nite and (ii) the transformed integral converges because the sech2 factor,
which arises from the change of coordinates, decays faster than f (tanh(z/L)) grows as
|z| ’ ∞.

To recover an exponential rate of convergence, the transformed integral can be approx-
imated by any spectrally-accurate quadrature for the in¬nite interval. The simplest is the

trapezoidal rule, which becomes

h f (tanh(’P/(2L)) + f (tanh(P/(2L))
sech2 (P/(2L))

L 2
N ’1
f (tanh(’P/(2L) + hk/L)) sech2 (’P/(2L) + hk/L) (19.71)

where h ≡ P/N is the grid spacing in the new coordinate z and where P is chosen suf¬-
ciently large that the integrand is negligible at z = P/(2L). This formula is usually called
the “tanh” rule for obvious reasons.
When the integral is singular at a point x = a in the interior of the integration interval,
one should split the integral into two and apply the tanh rule separately to each piece. If
there is only a single singularity, one can replace the tanh mapping by a similar exponential
mapping onto a semi-in¬nite integral.

19.8.7 Sets and Solitaries
Sometimes, one needs to repeatedly integrate a large set of integrals which fall into a com-
mon class. For this case, the most ef¬cient procedure is to estimate how large N must be by
computing a few members of the class for several different N . One may then use the small-
est acceptable N for all the remaining integrals. Since most computations are performed
for ¬xed N , Gaussian quadrature is the cheapest way to integrate a whole set of integrals.
When one needs to integrate a single, solitary integral, then an adaptive procedure
like the Curtis-Clenshaw strategy is more ef¬cient. The non-Gaussian spectral quadrature
schemes described above are then very effective.
Chapter 20

Symbolic Calculations: Spectral
Methods in Algebraic
Manipulation Languages

“[A computer] can arrange and combine its numerical quantities exactly as if they were
letters or any other general sysmbols; and in fact it might be possible to bring out its
results in algebraic notation, were provisions made accordingly.”

“ Augusta Ada Byron, Countess of Lovelace (1844)

“Originally, the calculation [of the Bondi metric problem in general relativity] had
required something like six months by hand ... the computer calculation, using


. 94
( 136 .)