. 71
( 136 .)


or otherwise

EDT (L) ≡ max (|u(y)|) (17.2)

We are forced to use the phrase “estimate” because our de¬nition of the domain truncation
error is merely the maximum value of the exact solution outside the limits of the truncated
interval. The error in the approximate solution to a differential equation is another matter,
and may be much larger than EDT (L). However, the Assumption of Equal Errors asserts
that max |u(y) ’ uN (y)| ∼ O(EDT (L)) for most real-world problems, even though one can
contrive examples for which this is not true. In the rest of the chapter, we shall use EDT (L)
as our sole measure of the numerical error induced by truncating the computational inter-
There are a couple of subtleties even in this most obvious and unsubtle method. The
¬rst is that if L is ¬xed while the number of basis functions N tends to ∞, the error will not
converge to zero, but rather to EDT (L). It follows that to improve accuracy, both the size of
the domain L and the number of basis functions N must increase simultaneously.

This increase in L moves the poles and branch points of u(y) closer to the real y-axis
after [’L, L] has been rescaled to [’1, 1] for the application of a Chebyshev polynomial
basis (or whatever). This movement-of-the-singularities reduces the rate of convergence
of the Chebyshev series. Thus, when N is doubled, the logarithm of the total error is not
doubled, but is increased by a smaller amount because of the increase in L with N . As N
and L jointly increase, the error decreases to zero at a subgeometric rate in N . (If L could
be ¬xed, then the error would decrease geometrically in N as for a ¬nite interval, but this,
alas, is only a dream when the problem is posed on an unbounded domain.)
The second subtlety is that it is more ef¬cient to use an ordinary Fourier series instead
of Chebyshev polynomials, even though the former is prone to the wild oscillations of
Gibbs™ Phenomenon at y = ±L. The reason is that the amplitude of the Gibbs™ overshoot is
0.09|u(L) ’ u(’L)|/N (Boyd, 1988e), that is to say, the size of the wiggles is proportional to
but smaller than the magnitude of u(y) at the ends of the interval, i. e., the domain trunca-
tion error EDT (L). It follows that the effects of Gibbs™ Phenomenon are always smaller (by
O(1/N )) than the domain truncation error, and thus negligible. The Fourier basis is better
because for a given N , it has a higher density of points (by a factor of π/2) at the center
of the interval than the Chebyshev collocation grid, which wastefully clusters points near
y = ±L where u(y) is very small.
Domain truncation is a good because of its simplicity. However, it can improved further
by combining it with a mapping as shown by Weideman and Cloot (1990), and explained
later in this chapter. Furthermore, domain truncation is more sensitive to the choice of
the domain size L than the rational Chebyshev functions are to the choice of their map
parameter L (Boyd, 1982a, 1987a, 1987b).

17.2.2 Domain Truncation for Slowly-Decaying Functions
If u(y) decays as an inverse power of |y|, instead of exponentially, it is often still possible
to obtain high accuracy from domain truncation by altering the boundary conditions to
match the analytic, asymptotic form of the slowly-decaying solution. For example, Corral
and Jim´ nez(1995) developed an analytic irrotational correction to apply domain domain
to a ¬‚ow which is unbounded in one coordinate. Rennich and Lele(1997) have applied a
three-dimensional Fourier series to a vortical ¬‚ow which is unbounded in two coordinates,
extending Corral and Jim´ nez™s method to an additional dimension. The analytical bound-
ary condition makes it possible to use a computational domain of modest size even though
the decay of the irrotational part of the vorticity is rather slow, which would otherwise
demand a huge domain.
Domain truncation-with-matching is not the only option for slowly-decaying functions.
Algorithms that combine in¬nite interval basis sets with special asymptotic-mimicing basis
functions are the themes of Sec. 13 in this chapter (illustrated by the J0 Bessel function) and
Sec. 4 in Chapter 19 (wave scattering).

17.2.3 Domain Truncation for Time-Dependent Wave Propagation:
Sponge Layers
When waves evolve far from boundaries, it is common to idealize the geometry as an
in¬nite, boundary-free domain. Standard in¬nite interval basis sets have troubles because
the waves do not decay for large |y|; instead, some waves simply propagate out of the
region of interest. Domain truncation is a popular alternative. However, the usual Dirichlet
boundary conditions are equivalent to rigid boundaries, which would re¬‚ect the waves
back towards y = 0 like a mirror. The physics for large |y| must therefore be modi¬ed to
absorb the waves.

For simple wave equations (or for wave equations which asymptote to linear, con-
stant coef¬cient differential equations for large |y|), this modi¬cation can be accomplished
merely by altering the boundary conditions. For example, if the waves have the asymp-
totic spatial structure exp(±iky) for some wavenumber k where the plus sign describes
outward-propagating waves for positive y, then the right boundary condition can be mod-
i¬ed to
’ iku = 0 at y = L (17.3)

and similarly at the left boundary. This Robin condition ¬lters the inward-propagating
wave, thus allowing waves to pass through the boundary from small y to y > L without
Unfortunately, this simple strategy fails if the far ¬eld waves are a mixture of wavenum-
bers or vary non-sinusoidally. It is possible to generalize this “propagate-through-the-
boundary” idea; there has been much recent work on so-called “one-way wave equations”,
but it would take us too far a¬eld to describe it here.
A second option is to add an arti¬cial damping which is negligible for small y (where
the real physics happens) but is allowed to increase exponentially fast for large y. The
ends of the computational domain then soak up the waves like a sponge soaking up water;
this is called the “sponge layer” method. Any old boundary conditions at y = ±L will
do because the waves that re¬‚ect from the computational boundaries will never make it
back to vicinity of the origin. If, Berg, Christiansen and Skovgaard(1987) and Zampieri
and Tagliani(1997) are spectral examples.
The Fourier basis, used by If et al., is superior to a Legendre or Chebyshev basis because
the high resolution of the orthogonal polynomials near y = ±L is completely wasted when
these regions are “sponges”.
There is one peril in “sponge layers”: If the arti¬cial viscosity increases too rapidly,
accuracy may be awful. If the damping coef¬cient is discontinuous, then the wave equation
no longer has analytic coef¬cients and the spectral series will converge only at an algebraic
rate. However, even if the viscosity is an analytic function, the sponge layer may fail to
be a good absorber if the damping increases too rapidly with |y|. The reason is that if the
damping varies more rapidly than the wave, the wave will be partially re¬‚ected ” any
rapid variation in the wave medium or wave equation coef¬cients, whether dissipative or
not, will cause strong partial re¬‚ection. The physically-interesting region around y = 0 will
then be corrupted by waves that have re¬‚ected from the sponge layers.
Fortunately, it turns out that if is de¬ned to be the ratio of the length scale of the wave1
to the length scale of the dissipation, then the re¬‚ection coef¬cient is exponentially small
in 1/ (Boyd, 1998b). It follows that if both the length scale of the arti¬cial viscosity and the
domain size are systematically increased with the number of spectral coef¬cients N , then
the amplitude of re¬‚ected waves at y = 0 will decrease exponentially fast with N . Thus, if
implemented carefully, the sponge layer method is a spectrally-accurate strategy.

17.3 Whittaker Cardinal or “Sinc” Functions
Sir Edmund Whittaker (1915) showed that functions which (i) are analytic for all ¬nite real
y and (ii) which decay exponentially as |y| ’ ∞ along the real axis could be approximated
as sums of what are now called “Whittaker cardinal” or “sinc” functions. The convergence
is exponential but subgeometric with an exponential convergence index usually equal to 0.5.
1 The wave scale is 1/k, i. e., the local wavelength divided by 2π.

(That is, the error typically decreases proportional to exp(’qN 1/2 ) for some positive con-
stant q.)
f (y) ≈ (17.4)
f (yj ) Cj (y)


Cj (y; h) ≡ sinc[(y ’ jh)/h] (17.5)

sinc(y) (17.6)

and where the points of the associated collocation grid are

yj ≡ j h (17.7)

The coef¬cients of the expansion are simply the values of f (y) on the grid because the basis
functions are also cardinal functions with the property

Ci (yj ; h) = δij

What is most striking about (17.4) is that we must choose two independent parameters
to determine the approximation: (i) the truncation N and (ii) the scaling factor h. On a ¬nite
interval, only one parameter ” the truncation N ” would suf¬ce. This complication is not
special to sinc series: all spectral methods on an in¬nite or semi-in¬nite interval require us
to specify a scaling parameter of some sort in addition to N .
Spectral methods on an in¬nite or semi-in¬nite interval have other subtleties, too. A
truncated Chebyshev series is exact on a ¬nite interval only when f (x) is a polynomial;
what makes it useful is that one can approximate any function by a polynomial of suf-
¬ciently high degree. The analogous statement for sinc expansions is that the Whittaker
cardinal series is exact for so-called “band-limited” functions; sinc series are effective for
general functions that decay exponentially as |y| ’ ∞ because one can approximate ar-
bitrary functions in this class with arbitrary accuracy by a band-limited function of suf¬-
ciently large bandwidth.
“Band-limited” functions are de¬ned in terms of Fourier integrals, which are yet an-
other way of representing a function on an in¬nite interval:

F (ω) e’iωy dω
f (y) = √ [Fourier integral] (17.9)
2π ’∞

The Fourier transform F (ω) is analogous to the function a(n) that gives the coef¬cients of
an ordinary Fourier or Chebyshev series except that a(n) is meaningful only for integer n
whereas F (ω) is a continuous function of ω.
In the same way that series coef¬cients a(n) are computed by an integral which is the
product of f (t) with the basis functions,

F (ω) = √ f (y) eiωy dy [Inverse Transform] (17.10)
2π ’∞

An exponentially convergent numerical scheme to compute either the forward or inverse
transform is to replace the in¬nite limits of integration by ±W for some suf¬ciently large
W and then apply the rectangle rule or trapezoidal rule. We do not offer a separate chapter

or even a section on Fourier integrals because this algorithm is equivalent to a Fourier series
pseudospectral method with a very large period.
Communications engineers have introduced a special name for a truncated Fourier in-
tegral because electronic devices have maximum and minimum frequencies that physically
truncate the ω-integration.

A function fW (y) is said to be BAND-LIMITED with BANDWIDTH W if it can be repre-
sented as the TRUNCATED Fourier integral
F (ω)e’iωy dω
fW (y) = √ (17.11)
2π ’W

In the world of Fourier integrals, band-limited functions are the analogues of polyno-
mials. For any degree N , polynomials are “entire functions”, that is to say, have no singu-
larities for any ¬nite x in the complex x-plane. Even so, a polynomial may be an extremely
accurate approximation to a function that has poles or branch points. The resolution of this
apparent contradiction is that the polynomial is a good representation of f (x) only on an
interval which is free from singularities; the in¬nities must be elsewhere. In a similar way,
“band-limited” functions are always “entire functions”, too. This does not alter the fact
that any function f (y), even if it has singularities for complex y, can be approximated for
real y to within any desired absolute error by a band-limited function of suf¬ciently large
Eq. (17.4) and (17.8) show that the sinc series interpolates to f (y) at every point of
the grid. However, this hardly guarantees a good approximation; recall the Runge phe-
nomenon for polynomials. For band-limited functions, success is assured by the following.

If fW (y) is a BAND-LIMITED function of bandwidth W , then the Whittaker cardinal series is
an EXACT representation of the function if h ¤ π/W where h is the grid spacing:

πj Wy
sinc (17.12)
fW (y) = f
W π

PROOF: Stenger (1981) or Dym and McKean (1972).
As noted above, band-limited functions are rather special. However, under rather mild
conditions on f (y), one can prove that F (ω) will decrease exponentially with |ω|, so the
error in writing f (y) ≈ fW (y) decreases exponentially with W , too. Thus, (17.12) really
implies that we can use Whittaker cardinal functions to approximate any function f (y)
that has a well-behaved Fourier transform.

De¬nition 38 (SINC EXPANSION OF A FUNCTION) An approximation to a function f (y)
of the form
y ’ jh
f (hj) sinc (17.13)
f (y) = + EW (h) + EDT (N h)

is called the SINC EXPANSION of f (y) where h is the grid spacing.

The error EW (h) is called the “grid-spacing” or “bandwidth” error. It arises because the grid
points are spaced a ¬nite distance h apart. Equivalently, since the sinc expansion is exact even with
a non-zero h if f (y) is band-limited with a suf¬ciently large W , we may say that this error is caused
by (implicitly) approximating f (y) by a band-limited function fW (y) with W = π/h.
The error EDT (N h) is the “grid-span” or “domain truncation” error. It arises because we must
truncate the in¬nite series to limits of ±N/2 for some ¬nite N . This is equivalent to interpolating
f (y) on a grid of (N + 1) points that span only a portion of the interval [’∞, ∞]. This error is
O(f [hN/2]).

Figure 17.1: The —™s mark the evenly spaced interpolation points used with sinc functions
to solve problems on y ∈ [’∞, ∞]. It can be shown that for most situations, best results

are obtained by choosing a grid spacing h proportional to N where N is the number of
grid points. This implies that interval spanned by the grid points ” the “grid span” ”


. 71
( 136 .)