or otherwise

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

(ii)

|y|≥L

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-

val.

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.

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

340

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

e

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-

e

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.

17.3. WHITTAKER CARDINAL OR “SINC” FUNCTIONS 341

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

‚u

’ iku = 0 at y = L (17.3)

‚y

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

re¬‚ection.

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

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

342

(That is, the error typically decreases proportional to exp(’qN 1/2 ) for some positive con-

stant q.)

N/2

f (y) ≈ (17.4)

f (yj ) Cj (y)

j=’N/2

where

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

sin(πy)

≡

sinc(y) (17.6)

πy

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

(17.8)

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:

∞

1

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,

∞

1

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

17.3. WHITTAKER CARDINAL OR “SINC” FUNCTIONS 343

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.

De¬nition 37 (BAND-LIMITED FUNCTIONS)

A function fW (y) is said to be BAND-LIMITED with BANDWIDTH W if it can be repre-

sented as the TRUNCATED Fourier integral

W

1

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

bandwidth.

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.

Theorem 32 (SHANNON-WHITTAKER SAMPLING THEOREM)

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

’j

sinc (17.12)

fW (y) = f

W π

j=’∞

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

N/2

y ’ jh

f (hj) sinc (17.13)

f (y) = + EW (h) + EDT (N h)

h

j=’N/2

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

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

344

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