(2.41)

u(x, t = 0) = Q(x)

Although the diffusion equation is a partial differential equation, it is one-dimensional and

linear. The particular initial condition

Q(x) = x(π ’ x) (2.42)

satis¬es the boundary conditions and is a polynomial which can be exactly represented by

the sum of the lowest three Chebyshev polynomials. Nevertheless, u(x, t) is singular.

The reason is that the solution to the problem as posed is actually the restriction to

the interval x ∈ [0, π] of the diffusion equation on the in¬nite interval, subject to an initial

2.8. LOCATION OF SINGULARITIES 39

condition which is spatially periodic and antisymmetric with respect to both the origin and

x = π. We can create such an initial condition by either (i) expanding the initial Q(x) as a

sine series or (ii) de¬ning it directly as

∀x ∈ [’∞, ∞]

P (x) = sign(sin(x))Q(x) (2.43)

where the sign function equals one when its argument is positive and is equal to minus

one when its argument is negative. Fig 2.13 shows the initial condition (Eq. 2.42) and its

second derivative. The latter has jump discontinuities at x = ±0, π, 2 π, . . . .

At t = 0, these discontinuities cause no problems for a Chebyshev expansion because

the Chebyshev series is restricted to x ∈ [0, π] (using Chebyshev polynomials with argu-

ment y ≡ (2/π)(x ’ π/2)). On this interval, the initial second derivative is just the constant

’2. For t > 0 but very small, diffusion smooths the step function discontinuities in uxx ,

replacing the jumps by very narrow boundary layers. As t ’ 0+, the layers become in-

¬nitely thin, and thus a Chebyshev approximation for any ¬xed truncation N must converge

slowly for suf¬ciently small t.

Fortunately, this pathology is often not fatal in practice because these diffusive bound-

ary layers widen rapidly so that the evolution for later times can be easily tracked with a

Chebyshev spectral method for small or moderate N . Indeed, many scientists have hap-

pily solved the diffusion equation, graphing the answer only at longish time intervals, and

missed the narrow transient boundary layers with no ill effect.

EXAMPLE: One-Dimensional Wave Equation

(2.44)

utt = uxx , u(0) = u(π) = 0

Q Qxx

2 2

1 1

0 0

-1 -1

-2 -2

-5 0 5 -5 0 5

u(x,t=1/10) uxx(x,t=1/10)

2 2

1 1

0 0

-1 -1

-2 -2

-5 0 5 -5 0 5

Figure 2.13: A representative solution to the diffusion equation. Upper panels: Initial con-

dition (left) and its second derivative (right). Bottom panels: Solution at t = 1/10 (left) and

its second derivative (right).

CHAPTER 2. CHEBYSHEV & FOURIER SERIES

40

f fxx

2 2

1 1

0 0

-1 -1

-2 -2

0 2 4 6 0 2 4 6

u(x,t=0.5) uxx(x,t=0.5)

2 2

1 1

0 0

-1 -1

-2 -2

0 2 4 6 0 2 4 6

Figure 2.14: Wave equation. Upper panels: Initial condition (left) and its second derivative

(right). Bottom panels: Solution at t = 1/2 (left) and its second derivative (right). The

singularities, which are initially at x = 0, π, 2π, etc., propagate both left and right.

(2.45)

u(x, t = 0) = f (x), ut (x, t = 0) = g(x)

The particular initial condition

f (x) = x(π ’ x), (2.46)

g(x) = 0

is the same as for the diffusion example, and also yields a singular solution

Again, the reason is that the problem as posed is actually the restriction to the interval

x ∈ [0, π] of the same equation on the in¬nite interval with the initial condition:

P (x) = sign(sin(x))f (x) ∀x ∈ [’∞, ∞] (2.47)

Fig 2.14 shows the initial condition and its second derivative. The latter has jump disconti-

nuities at x = ±0, π, 2 π, . . . .

The general solution to the wave equation (ut (x, t = 0) = 0) is

u(x, t) = (1/2){ f (x ’ t) + f (x + t) } (2.48)

Instead of diffusing away, the jumps in the second derivative, initially at t = mπ where m is

an any integer, propagate both left and right. Thus, at the later time illustrated in the lower

panels of Fig. 2.14, the exact solution has a second derivative which is discontinuous at two

points in the interior of x ∈ [0, 1], even though the initial condition was super-smooth: a

quadratic polynomial.

The singularities of the diffusion and wave equations are similar to those of Poisson™s

equation in that they, too, are located at the corners, at least initially. However, the corners

are now in space-time, that is, at the corners of a domain in the x ’ t plane.

2.9. FACE: INTEGRATION-BY-PARTS BOUND 41

In contrast to the singularities of an elliptic equation, the corner singularities of a hy-

perbolic or parabolic equation are usually unphysical. The reason is that the choice of the

initial time is usually arbitrary. In weather forecasting, for example, the forecast is usually

begun at midnight Greenwich mean time. Why should the singularities be located at the

boundary of the atmosphere at this time and no other? Would physical singularities move

if the forecast were begun an hour earlier or later? Boyd and Flyer(1999) discuss how initial

conditions can be slightly modi¬ed to satisfy “compatibility” conditions and thus eliminate

the space-time corner singularities.

These three examples should chill the blood of any aspiring computational scientist.

Albert Einstein noted ”Subtle is the Lord”, which is the great man™s way of saying that your

computer program is always trying to kill you. Fear is a very good thing for a numerical

analyst. One may wish mournfully for faster silicon, but the only absolutely fatal disease

to a scientist is a de¬ciency of thinking.

With thinking, spectral methods do very well indeed, even for nasty problems. In the

rest of the book, we explain how.

2.9 FACE: Integration-by-Parts Bound on Fourier Coef¬cients

The coef¬cients of the complex-exponential form of a Fourier series are:

π

1

(2.49)

cn = f (x) exp(’in x) dx

2π ’π

If we integrate by parts, repeatedly integrating the exponential and differentiating f (x),

we obtain without approximation after (J + 1) steps:

J j+1

1 i

f (j) (π) ’ f (j) (’π)

j+n

cn = (’1)

2π n

j=0

J+1 π

1 i

’ f (J+1) (x) exp(’in x ) dx (2.50)

+

2π n ’π

where f J+1 (x) denotes the (J + 1)-st derivative of f (x). (Note that these integrations-by-

parts assume suf¬cient differentiability of f (x).) The integral can be bounded by the length

of the integration interval times the maximum of the integrand, which is the maximum of

the absolute value of | f J+1 |. In the limit n ’ ∞ for ¬xed J, it follows that the integral is

O(1/nJ+1 ); ignoring it gives what Lyness (1971, 1984) has dubbed the “Fourier Asymptotic

Coef¬cient Expansion” (FACE):

J j+1

1 i

f (j) (π) ’ f (j) (’π) + O(n’(J+1) ), n ’ ∞, ¬xed J (2.51)

cn ∼ j+n

(’1)

2π n

j=0

This expansion is an asymptotic series, diverging in the limit of ¬xed n, J ’ ∞.

By taking real and imaginary parts and recalling that the cosine and sine coef¬cients

are related to the complex coef¬cients via (for real-valued f (x))

bn = ’ 2 cn ∀n>0 (2.52)

an = 2 (cn );

CHAPTER 2. CHEBYSHEV & FOURIER SERIES

42

one obtains an equivalent form of the FACE as (n > 0)

J

f (2j+1) (π) ’ f (2j+1) (’π)

1

+ O(n’(2J+4) ), n ’ ∞, ¬xed J (2.53)

an ∼ (’1)n+j

n2j+2

π j=0

J

f (2j) (π) ’ f (2j) (’π)

1

+ O(n’(2J+3) ), n ’ ∞, ¬xed J

bn ∼ (’1)n+1+j (2.54)

n2j+1

π j=0

If a function has discontinuities in the interior of the interval x ∈ [’π, π], then the inte-

grals can be split into segments and integration-by-parts applied to derive similar asymp-

totic expansions.

Although sometimes used to accelerate the calculation of Fourier coef¬cients for func-

tions whose series lack exponential convergence, the FACE is most important as the proof

of the following theorem.

Theorem 4 (INTEGRATION-BY-PARTS COEFFICIENT BOUND) If

1.

f (π) = f (’π), f (1) (π) = f (1) (’π), ..., f (k’2) (π) = f (k’2) (’π) (2.55)

2. f (k) (x) is integrable

then the coef¬cients of the Fourier series

∞ ∞

(2.56)

f (x) = a0 + an cos(nx) + bn sin(nx)

n=1 n=1

have the upper bounds

| an |¤ F/nk ; | bn |¤ F/nk (2.57)

for some suf¬ciently large constant F , which is independent of n.

An equivalent way of stating the theorem is that, if the two conditions above are satis¬ed, then

the algebraic index of convergence is at least as large as k.

Notes: (a) f (k) denotes the k-th derivative of f(x) (b) The integrability of f (k) requires that f (x),

f (1) (x), . . . , f (k’2) (x) must be continuous.

PROOF: Under the conditions of the theorem, the ¬rst (k ’ 1) terms in the FACE are

zero. For suf¬ciently large n, the lowest nonzero term in the series will dominate, implying

that the Fourier coef¬cients are asymptotically O(nk ). Q. E. D.

A few remarks are in order. First, the usual way to exploit the theorem is to integrate-

by-parts as many times as possible until lack of smoothness or a mismatch in boundary

values forces us to stop. At that point, the number k that appears in the theorem is a lower

bound on the algebraic convergence index de¬ned earlier.

However, the index of convergence need not be an integer. For example, if f (x) is a

periodic function with a cube root singularity for real x, such as sin1/3 (x), one integration

by parts shows that the Fourier coef¬cients decrease at least as fast as O(1/n). Because

the second derivative of a cube root singularity has a branch point proportional to x’5/3 ,

2.9. FACE: INTEGRATION-BY-PARTS BOUND 43

which is not integrable, we cannot integrate-by-parts a second time. However, one can

show by other arguments that the integral

π

1

f (1) (x) exp(’inx) dx (2.58)

n ’π

is actually O(1/n4/3 ). (It would take us too far a¬eld to give the detailed proof of this.) The

theorem gives the largest integer which is less than or equal to the actual algebraic index

of convergence. For the sawtooth function and the half-wave recti¬er function, the index

is an integer and then the theorem gives it precisely.

Second, for pure sine or pure cosine series, every other integration-by-parts is triv-

ial. (Note that the FACE for the cosine coef¬cients proceeds only in even powers of n and

odd derivatives of f (x) while the sine coef¬cient expansion involves only odd powers of

n but even derivatives of f (x).) For example, if f (x) = f (’x), then the function is sym-

metric about x = 0 and all its sine coef¬cients are zero. The even order boundary terms

for the cosine coef¬cients vanish independent of whether f (2k) (π) = f (2k) (’π); as if to re-

inforce the point, the symmetry condition ensures that all the even derivatives are equal

at the boundaries anyway. Thus, for a cosine series, it is suf¬cient to examine whether

f (2k’1) (π) = f (2k’1) (’π). Similarly, for a sine series, boundary matching of the odd order