0.303 N 4

Chebyshev: Gauss grid accurate for j < N (2/π)

0.102 N 4

Legendre: Gauss grid accurate for j < N (2/π)

CHAPTER 12. IMPLICIT SCHEMES & THE SLOW MANIFOLD

238

5

Circles: Chebyshev

10

Fast Modes

4

10 Slow Manifold

3

10

Solid: exact eigenvalues

2

10

1

10

0

10

0 5 10 15

j

Figure 12.4: The exact eigenvalues of the second derivative operator with homogeneous

Dirichlet boundary conditions on x ∈ [’1, 1] are shown as the heavy solid line; the cor-

responding eigenvalues of the Chebyshev collocation matrix are the circles-with-thin-line.

Roughly (2/π) ≈ 0.64 of the Chebyshev eigenvalues ” roughly 10 of the 16 eigenmodes

illustrated ” are good approximations. Eigenvalues 11 to 16 are way off, and there is no

point in accurately tracking the time evolution of these modes of the Chebyshev pseu-

dospectral matrix because they have no counterpart in the solution of the diffusion equa-

tion. The slow manifold is composed of the modes left of the vertical dividing line whose

Chebyshev eigenvalues do approximate modes of the diffusion equation. The largest nu-

merical eigenvalue is about 56,000 whereas the largest accurately approximated eigenvalue

is only about 890, a ratio of roughly 64, and one that increases quadratically with N . Only

symmetric eigenmodes are shown.

The exact solution to the diffusion equation is of the same form except that the eigen-

vectors are trigonometric functions and the eigenvalues are different as given by the ¬rst

line of Table 12.1. The good news of the table is: the ¬rst (2/π)N eigenvalues are well-

approximated by any of the three spectral methods listed. This implies that we can follow

the time evolution of as many exact eigenmodes as we wish by choosing N suf¬ciently

large.

The bad news is: the bad eigenvalues are really bad because they grow as O(N 4 ) even

though the true eigenvalues grow only as O(N 2 ). In the parlance of ODE theory, the

Chebyshev-discretized diffusion equation is very “stiff”, that is, it has an enormous range

in the magnitude of its eigenvalues. The stability limit (set by the largest eigenvalue) is

very small compared to the time scale of the mode of smallest eigenvalue.

The slow manifold is the span of the eigenmodes which are accurately approximated

by the discretization, which turns out to be the ¬rst (2/π)N modes for this problem. The

fast-and-we-don™t-give-a-darn manifold is the subspace in con¬guration space which is

spanned by all the eigenmodes that are are poorly resolved by N Chebyshev polynomials

(Fig. 12.4).

However, the physics is also relevant: if only the steady-state solution matters, than the

slow manifold is the steady-state, and the time evolution of all modes can be distorted by

a long time step without error. The physics and numerics must be considered jointly in

identifying the slow manifold.

12.10. INITIALIZATION 239

10

8

6

4

2

0

-2

-4

-6

-8

-10

0 1 3 4 5

2

TIME (days)

Figure 12.5: The time evolution of the north-south wind at 50.40 N. on the Greenwich

meridian for the barotropic forecast of Richardson (1922). The dotted curve is the forecast

made from the raw data without initialization. The solid curve is not the result of a running

time average, but rather is the forecast from the initial conditions obtained by adjusting the

raw data onto the slow manifold. [Redrawn from Fig. 8c of Lynch(1992) with permission

of the author and the American Meteorological Society.]

12.10 Initialization

Meteorological jargon for the statement that “the large scale dynamics is near the slow

manifold” is to say that the ¬‚ow is “in quasi-geostrophic balance”. Unfortunately, raw

observational data is corrupted by random instrumentation and sampling errors which,

precisely because they are random, do not respect quasi-geostrophic balance. Initializ-

ing a forecast model with raw data generates a forecast with unrealistically large gravity

wave oscillations superimposed on the slowly-evolving large-scale dynamics (Daley, 1991,

Lynch, 1992, Fig. 12.5). The remedy is to massage the initial conditions by using an algo-

rithm in the following category.

De¬nition 28 (SLOW MANIFOLD INITIALIZATION)

A numerical algorithm for adjusting unprocessed initial values from a point in con¬guration space

off the slow manifold to the nearest point on the slow manifold is called a SLOW MANIFOLD

INITIALIZATION or simply (especially in meteorology) an INITIALIZATION procedure.

Fig. 12.6 schematically illustrates initialization onto the slow manifold.

Meteorological simulations of large-scale weather systems can be subdivided into fore-

casting and climate modelling. For the latter, the model is run for many months of time

until it “spins up” to statistical equilibrium. The time-averaged temperature and winds

then de¬ne climate of the model. The choice of initial conditions is irrelevant because the

¬‚ow will forget them after many months of integration; indeed, independence of the initial

conditions is an essential part of the very de¬nition of climate.

CHAPTER 12. IMPLICIT SCHEMES & THE SLOW MANIFOLD

240

Raw Data

Adjusted

Initial Conditions

Figure 12.6: Schematic of an initialization onto the slow manifold, which is shown as the

cross-hatched surface in a three-dimensional con¬guration space. (A real con¬guration

space might have millions of dimensions where each coordinate is a spectral coef¬cient or

grid point value.)

A slow manifold initialization is unnecessary for a climate model. The inherent dissipa-

tion of the model gradually purges the fast transients, like white blood cells eating invading

germs. Similarly, the backwards Euler scheme, which damps all frequencies but especially

high frequencies, will pull the ¬‚ow onto the slow manifold for almost any physical system,

given long enough to work. The BE algorithm, viscosity, and other physical or computa-

tional damping mechanisms are a kind of “Poor Man™s Slow Manifold Initialization”.

When the ¬‚ow must be slow from the very beginning, however, stronger measures are

needed. Meteorologists have tried a variety of schemes including:

1. Modifying the partial differential equations so that the approximate equations sup-

port only “slow” motions (quasi-geostrophic and various “balanced” systems, re-

viewed in Daley, 1991).

2. “Dynamic initialization”, which is integrating a few steps forward in time, followed

by an equal number back to t = 0, using special algorithms designed to damp high

frequencies, rather than accurately integrate in time (Daley, 1991, Chap. 11).

3. Time-frequency ¬ltering via Laplace transforms (Lynch and Huang, 1992).

4. Method of multiple scales [“Normal mode initialization”] (Machenhauer, 1977, Baer,

1977, Baer and Tribbia, 1977).

Without assessing the relative merits of these schemes for operational forecasting, we shall

describe the method of multiple scales, or “normal mode initialization” as meteorologists

call it, because it provides the theoretical underpinning for the Nonlinear Galerkin meth-

ods.

12.11. THE METHOD OF MULTIPLE SCALES(BAER-TRIBBIA) 241

12.11 The Method of Multiple Scales:

The Baer-Tribbia Series

When a partial differential equation with both fast and slow modes is discretized through

Galerkin™s method using the normal modes as the spectral basis, the result is a system of

coupled ordinary differential equations in time that can be written

St + i ωS . — S = fS (S, F )

Ft + i ωF . — F = fF (S, F )

after partitioning the solution u into its slow and fast modes, S and F , respectively. (In

meteorology, S is the vector of amplitudes of the Rossby modes while F is the amplitude

of the gravity waves.) The symbol .— is used to denote elementwise multiplication of two

vectors, that is, ωS . — S is a vector whose j-th element is the product of the frequency of

the j-th mode, i. e., the j-th element of ωS , with the j-th element of S. The obvious way

to integrate on the slow manifold is to simply ignore the fast modes entirely and solve the

reduced system

St + i ωS . — S = fS (S, 0) (12.54)

Machenhauer (1977) pointed out that this is a rather crude approximation. The non-

linear interaction of the slow modes amongst themselves will create a forcing for the fast

modes. If the time scale for the slow modes is O( ) compared to the fast frequencies, then

the second line of Eq. (12.54) implies that

(fF (S, 0))j

Fj ≈ ’i (12.55)

ωj

to within a relative error of O( ). This forced motion in the fast modes varies slowly with

time because it is forced only by the slow modes. Machenhauer showed that initializa-

tion schemes for weather forecasting were signi¬cantly improved by computing the initial

values of the gravity modes from Eq.(12.55). If the gravity modes are initialized to zero,

then the forecast must include (thrice-accursed!) high frequency gravity waves to cancel

(at t = 0) the slow particular solution given by Eq.(12.55), and the forecast will resemble

the wiggly dotted curve in Fig. 12.5.

Baer(1977) and Baer and Tribbia (1977) observed that Machenhauer™s approximation

was really just the lowest order of a singular perturbation expansion in , derivable by the

so-called “method of multiple scales”. The underlying idea can be illustrated by computing

the perturbation series for a linear, uncoupled ODE:

du

(12.56)

= »u + f ( t)

dt

which is just Eq.(12.5) again except that we make one crucial assumption: that the forcing

varies slowly with time as measured by the small parameter which appears inside its

argument. (Implicitly, » ∼ O(1), but this is not crucial; the effective perturbation parameter

will turn out to be the ratio /» which is the ratio of the “fast” and “slow” times scales for

the problem.)

It is convenient to introduce the slow time variable

T≡ t (12.57)

CHAPTER 12. IMPLICIT SCHEMES & THE SLOW MANIFOLD

242

The ODE becomes

du

(12.58)

= » u + f (T )

dT

It is then clear that lowest nontrivial approximation is the same as Machenhauer™s: ne-

glecting the time derivative so that there is an approximate balance between the two terms

on the right-hand-side. This approximation, however, is sensible only for motion on the

slow manifold or within O( ) of the slow manifold; the homogeneous solution is the fast-

evolving function exp(»t), which is poorly approximated by neglecting the time derivative.

De¬ning

N

uj j

(12.59)

u(T ; ) =

j=0

we obtain

= ’f (T )/»

u0 (T )

du0 /dT

u1 (T ) =

»

df /dT

=’

»2

du1 /dT

2

u (T ) =

»

d f /dT 2

2

=’

»3

and in general