. 11
( 136 .)


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

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

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

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.

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.

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:
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) (’π)
cn = (’1)
2π n
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
2π n

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

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

f (2j+1) (π) ’ f (2j+1) (’π)
+ O(n’(2J+4) ), n ’ ∞, ¬xed J (2.53)
an ∼ (’1)n+j
π j=0

f (2j) (π) ’ f (2j) (’π)
+ O(n’(2J+3) ), n ’ ∞, ¬xed J
bn ∼ (’1)n+1+j (2.54)
π 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.



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

which is not integrable, we cannot integrate-by-parts a second time. However, one can
show by other arguments that the integral
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


. 11
( 136 .)