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.

EXAMPLE THREE: LK QUINTET

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

=

Ut

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

CHAPTER 12. IMPLICIT SCHEMES & THE SLOW MANIFOLD

234

Slow Modes Fast Modes

-3

x 10

0.1

5

0.08

4

0.06

3

0.04 2

0.02 1

0

0

-1

-0.02

-2

-0.04 -3

-0.06 -4

-5

-0.08

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 =

’z

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

(12.42)

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

(12.43)

c=4

12.8. SLOW MANIFOLD: DEFINITION & EXAMPLES 235

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)

1

0.5

0

-0.5

-1

0 0.5 1 1.5 2 2.5

Slow (thick curve) and fast (thin curve) parts

1

0.5

0

-0.5

-1

0 0.5 1 1.5 2 2.5

time

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

CHAPTER 12. IMPLICIT SCHEMES & THE SLOW MANIFOLD

236

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

∞

dak

’ 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

Matrices

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

h

„max = d ν h2

[Waves] [Diffusion] (12.47)

„max = d

cmax

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

(12.48)

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

12.9. NUMERICALLY-INDUCED SLOW MANIFOLDS 237

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

(12.50)

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

dbj

= ’»j bj + gj (12.52)

dt

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

N

(12.53)

S(t) = f + bj (0) exp(’ »j t) ej

j=1

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

Exact

jπ

(N + 1)2 sin2 N2

2d order Finite Difference 2(N +1)

0.048 N 4