. 51
( 136 .)


Chebyshev: Lobatto grid accurate for j < N (2/π)
0.303 N 4
Chebyshev: Gauss grid accurate for j < N (2/π)
0.102 N 4
Legendre: Gauss grid accurate for j < N (2/π)

Circles: Chebyshev
Fast Modes
10 Slow Manifold


Solid: exact eigenvalues


0 5 10 15

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

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.

0 1 3 4 5
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.

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.

Raw Data

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-

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)

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:

= »u + f ( t)
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)

The ODE becomes
= » u + f (T )
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.

uj j
u(T ; ) =

we obtain

= ’f (T )/»
u0 (T )
du0 /dT
u1 (T ) =
df /dT
du1 /dT
u (T ) =
d f /dT 2
and in general


. 51
( 136 .)