. 50
( 136 .)


depth. First, the more general problem du/dt + iωu = f ( ωt) can be converted into the
form of (12.39) merely by rescaling the time through t ’ ωt. Second, as noted in the ¬rst
section, any linear partial differential equation with time-independent coef¬cients can be
converted into the form of (12.39) by shifting to a basis of eigenfunctions of the matrix of
the discretized spatial derivatives.
Third, even a nonlinear problem can be cast in the form of (12.39) if the forcing f ( t) is
reinterpreted as a symbol for the nonlinear terms that couple different eigenmodes of the
linear part of the partial differential equation together. This point of view is popular in
weather forecasting because the nonlinear terms are an order of magnitude smaller than
the linear terms (for gravity modes); thus, global weather can be conceptualized as a sea of
independently propagating wave modes with weak wave-wave coupling.
Actually, the nonlinear interactions of the Rossby waves among themselves are rather
strong, but to a ¬rst approximation the gravity waves are linear oscillators forced by the
nonlinear Rossby-Rossby interactions. These depend very weakly on the gravity waves
themselves, so that these Rossby-Rossby couplings may be regarded as external forcing for
the gravity waves.
These ideas are clari¬ed by the following.


This model is described by Lorenz and Krishnamurthy (1987). It is obtained by truncat-
ing an atmospheric model to the most ruthless extreme that allows nontrivial interactions:
three Rossby modes and two gravity waves. If we denote the (large) Rossby amplitudes by
uppercase letters (U, V, W ) and the (small) gravity mode amplitudes by (x, z) “ note that
these are Fourier coef¬cients and not velocities or coordinates “ then the LK Quintet is

’V W + βV z ’ aU
U W ’ βU z ’ aV + F
Vt =
’U V ’ aW (12.40)
Wt =
’z ’ ax
xt =
bU V + x ’ az
zt =

where subscript t denotes time differentiation and where (a, F, b, β) are constant parame-
ters. We shall mostly concentrate on the “inviscid” quintet, which is the special case that

Slow Modes Fast Modes
x 10
0.04 2
0.02 1
-0.04 -3
-0.06 -4
0 20 40 60 80 0 20 40 60 80

Figure 12.2: Example of solution to the “one-way coupled” inviscid LK Quintet (b = 1, β =
0 for the initial conditions [U, V, W, x, z] = [1/10, 0, ’9/100, 0, 0]. Left panel: amplitudes of
the three Rossby waves. (i) U : Thickest curve (never changes sign). V : Medium thickness
(“tent-shaped”). W : Thinnest curve. Right panel: amplitudes of the two gravity waves
with x(t) as the thicker, larger amplitude curve.

the damping coef¬cient a and the forcing F are both zero:

’V W + βV z
Ut =
U W ’ βU z
Vt =
’U V (12.41)
Wt =
xt =
zt = bU V + x

It is often convenient to set β = 0 to obtain the “one-way coupled” model because an exact
analytical solution has been given by Boyd(1994c).
A representative solution is shown in Fig. 12.2.

EXAMPLE FOUR: KdV Solitary Wave
The Korteweg-deVries equation is

ut + u ux + uxxx = 0

The general solution with spatially periodic boundary conditions consists of solitary waves
and quasi-sinusoidal travelling waves. The simplest solitary wave or “soliton” is an exact
nonlinear solution which steadily translates at a constant phase speed c:

usol (x, t) = 12 2 sech( [x ’ ct] + φ); 2

where and φ are constants. Strictly speaking, the soliton is de¬ned only for a spatially
unbounded interval, but its periodic generalization, called a cnoidal wave, is well approx-
imated by (12.43) when ∼ O(1) or larger.

KdV equation: real part of exp(i 5 x)
0 0.5 1 1.5 2 2.5
Slow (thick curve) and fast (thin curve) parts




0 0.5 1 1.5 2 2.5

Figure 12.3: The real part of the Fourier coef¬cient of exp(i5x) as a function of time for a
Korteweg-deVries (KdV) solution which is the sum of a solitary wave of unit phase speed
and a small amplitude sine wave of wave number k = 5. The amplitude has been scaled to
unity for the soliton contribution to this Fourier coef¬cient.

If the initial conditions are slightly perturbed from a soliton, then the ¬‚ow will evolve
as a solitary wave plus small amplitude sine waves which satisfy the linear KdV equation:

usine (x, t) = ± cos(k[x ’ k 2 t]) (12.44)

More general solutions with multiple solitons and large amplitude non-soliton waves are
also possible, but for our purposes, it will suf¬ce to consider a small amplitude solitary
wave perturbed by very small sine waves. The interaction between the soliton and the rest
of the ¬‚ow, although not zero, can be neglected to a ¬rst approximation.
Fig. 12.3 illustrates the time evolution of a typical Fourier coef¬cient. The slow compo-
nent, shown separately in the lower panel, is the solitary wave. The fast component, which
controls the maximum time step for an explicit method, is the perturbing sine wave. The
crucial point is that because the frequency of a sine wave of wavenumber k grows rapidly
as k 3 , these small waves are extremely fast-moving.
If the goal is to resolve the fast sine waves, then one might as well use an explicit
method. If K is the maximum wavenumber, then one must use „ ∼ O(1/K 3 ) to accurately
track its propagation.
If the goal is only to understand the evolution of the solitary waves ” historically true
of most KdV numerical studies ” then one can use an implicit or semi-implicit scheme

with a much longer time step. For example, the fourth order Runge-Kutta method, which
has a rather large stability limit, is stable only when

„ < 2.8/K 3 (12.45)

However, aK (t) varies due to the solitary wave as cos(Kct), which has a temporal period
of 2π/(K c). Thus, the largest timestep which is stable for RK4 is 2.2K 2 /c timesteps per
temporal period. For c = 1 and K = 16, this is 570 time steps per period, and the number
grows quadratically with K. Implicit algorithms allow one to accurately model the solitary
wave with a much longer time step than 1/570 of a temporal period.
There is actually a close relationship between this example and the forced linear oscil-
lator because the Fourier-Galerkin ODEs in time for the KdV equation are

’ ik 3 ak = ’ i (12.46)
ak’m am
dt m=’∞

This is identical in form to our second example (12.39) with ω = ’k 3 and with the nonlinear
terms in the in¬nite series serving as the forcing for the k-th Fourier coef¬cient. Of course,
one must be careful not to push the analogy too far because (12.46) is a coupled system, not
merely a single isolated ODE. However, for large k, the nonlinear self-interaction of ak can
be neglected and the forced linear oscillator is a good ¬rst approximation.

12.9 Numerically-Induced Slow Manifolds:
The Stiffness of Chebyshev Differentiation
As noted in the CFL Stability Limit Rule-of-Thumb in Chapter 9, the maximum timestep
is the same order of magnitude as the time scale for advection or diffusion or wave prop-
agation or whatever across the smallest distance h between two grid points. For wave
propagation and for diffusion, the limits are

„max = d ν h2
[Waves] [Diffusion] (12.47)
„max = d

where d, d are O(1) constants and where cmax is the speed of the fastest-moving waves
and ν is the diffusion or viscosity coef¬cient.
For Chebyshev and Legendre methods, these stability limits are very bad news because
these use highly nonuniform grids with the smallest h ” right next to the endpoints ”
being proportional to 1/N 2 . In contrast, an evenly spaced grid has h = 2/N everywhere.
It follows that explicit time-marching schemes are subject to very short time step limits
when applied to ODEs derived from Chebyshev or Legendre spatial discretizations. For
this reason, implicit time-marching schemes are very widely used with Chebyshev and
Legendre spectral methods. The longest chapter in the book (Chap. 15) will focus on
ef¬ciently solving the boundary value problems, such as those from implicit algorithms.
An even more pressing concern is: Are implicit schemes suf¬ciently accurate to resolve
everything that needs to be resolved when used with Chebyshev spatial discretizations?
An example is useful: the one-dimensional diffusion equation with homogeneous Dirichlet
boundary conditions is

ut = uxx + f (x); u(’1) = u(1) = 0; u(x, 0) = Q(x)

where for simplicity f (x) is independent of time.
We can exploit symmetry by splitting u = S(x, t) + A(x, t) where S and A are u™s sym-
metric and antisymmetric parts with respect to x = 0. We can de¬ne symmetric cardinal
functions as a basis for S via

φj (x) ≡ Cj (x) + C’j (x) (12.49)

where the Cj (x; 2N ) are the usual Chebyshev cardinal functions on the (2N + 2)-point
Lobatto grid, numbered so that the cardinal functions which are one at the endpoints are
C±(N +1) .
The equation for the symmetric part is then reduced by the Chebyshev cardinal function
spatial discretization to

St = D S + f

where the elements of the N -dimensional column vector S(t) are the time-dependent val-
ues of S(x, t) at each of the positive collocation points, excluding the endpoint x = 1, which
is omitted so that S(1, t) = 0 for all t to satisfy the boundary conditions, and f is similarly
the grid point values of f (x) for x ≥ 0. The elements of the square matrix D, which is the
discrete representation of the second derivative operator, are

Dij ≡ φj,xx (xi ), (12.51)
i, j = 1, 2, . . . , N

where the boundary conditions have been implicitly included by omitting the cardinal
function and grid point associated with x = 1.
If we compute the eigenvectors and eigenvectors of D, we can switch to an eigenvector
basis. In terms of this, the second derivative matrix is diagonal and (12.50) reduces to the
uncoupled set of equations

= ’»j bj + gj (12.52)
which is identical to (12.5) except that we have changed the sign of the eigenvalue so that »j
is positive for diffusion; the {bj (t)} and {gj }are the spectral coef¬cients in the eigenvector
basis. The solution is
S(t) = f + bj (0) exp(’ »j t) ej

Table 12.1: Eigenvalues of the Second Derivative with Homogeneous Dirichlet Boundary
Conditions Including Both Symmetric and Antisymmetric Modes

Method Eigenvalues »j , j = 1, 2, ..., N max(»)
j 2 π2 / 4 N 2 π2 / 4

(N + 1)2 sin2 N2
2d order Finite Difference 2(N +1)
0.048 N 4


. 50
( 136 .)