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)

∞ π

sin(t)

f (L cot2 (t/2)) 2L (19.55)

I = f (x)dx = dt

[1 ’ cos(t)]2

0 0

The quadrature is

N

≡ wj f (L cot2 (tj /2)) (19.56)

IN

j=1

≡ πj/(N + 1), (19.57)

tj j = 1, 2, . . . , N

N

[1 ’ cos(mπ)]

sin(tj ) 2

≡ 2L (19.58)

wj sin(mtj )

[1 ’ cos(tj )]2 N + 1 m

m=1

19.8. SPECTRALLY-ACCURATE QUADRATURE METHODS 457

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

∞ π

1

(19.60)

I = f (x)dx = f (L cot(t)) 2L dt

sin2 (t)

0 0

The quadrature is

N

≡ (19.61)

IN wj f (L cot(tj ))

j=1

≡ πj/(N + 1), (19.62)

tj j = 0, 1, . . . , N + 1

Lπ/[sin2 (tj )(N + 1)], 0 < j < (N + 1)

≡ (19.63)

wj

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.

Th

Theorem 38 (TRAPEZOIDAL RULE ERROR for PERIODIC INTEGRANDS)

trapezoidal rule approximation to the integral of a function f (x) on the interval x ∈ [’P/2, P/2]

is

±

f (’P/2) + f (P/2)

N ’1

P/2

f (x)dx ≡ IN = h (19.64)

+ f (’P/2 + hk)

2

’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

CHAPTER 19. SPECIAL TRICKS

458

then the error in the trapezoidal rule is given without approximation as

P/2

EN ≡ f (x)dx ’ IN = ’P {±N + ±2N + ±3N + . . . } (19.66)

’P/2

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)

2

’∞ k=1

where

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

|x|>L

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

problem.

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)

19.8. SPECTRALLY-ACCURATE QUADRATURE METHODS 459

Error in Tanh Rule

0

10

1

(1-x 2)-3/4 dx

«

-1

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

k

-5

10

1/2

-10

sum of sech(zk)

10

0 100 200 300 400

N

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

∞

1

1

f (tanh(z/L)) sech2 (z/L)dz

I≡ (19.70)

f (x) dx =

L

’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

CHAPTER 19. SPECIAL TRICKS

460

trapezoidal rule, which becomes

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

sech2 (P/(2L))

≡

IN

L 2

N ’1

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

+

k=1

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