. 57
( 136 .)


The stability limit of the modi¬ed equation is then controlled by nonlinear term, which
involves only a ¬rst derivative and demands a time step inversely proportional to kmax
itself rather than to its cube.
The exponential-of-derivatives operator looks forbidding, but it can be interpreted in
the usual way as the exponential of the matrix that discretizes the third spatial derivative,
which in turn can be interpreted as multiplication of the k-th eigenvector of the derivative
matrix by the exponential of (’t) times the k-th eigenvalue. In a Fourier basis, this is
trivial because the basis functions, exp(ikx), are the eigenvalues of the derivative operator.
It follows that for any function v(x, t) with Fourier coef¬cients vk (t),

exp(i k 3 t) vk exp(ikx) (14.24)
exp(’‚xxx t) v(x, t) =

To put it another way, the integrating factor method for the KdV equation is equivalent to
transforming the unknowns of the ODE system which comes from the usual Galerkin or
pseudospectral discretization. In place of the Fourier coef¬cients ak (t) of u(x, t), use the
modi¬ed Fourier coef¬cients (of v(x, t)) given by

vk (t) ≡ exp(’ik 3 t) ak (t) (14.25)

Unfortunately, this is really a terrible idea when the prime target is solitary waves. The
simplest counterexample is the special case that the exact solution is a solitary wave of
phase speed c. The exact Fourier coef¬cients ” before the integrating factor ” are

ak = constant exp(’ikct) (14.26)

and thus vary rather slowly in the sense that the frequency of the highest wavenumber
in the truncation grows only linearly with kmax . In contrast, the time dependence of the
transformed coef¬cients is

vk (t) = constant exp(’ik 3 t ’ ikct) (14.27)

which varies rapidly in the sense that the highest frequency after transformation grows as
3 3
kmax . To accurately compute vkmax (t), a time step proportional to kmax is necessary ”
in other words, a time step as short as the explicit stability limit for the original, untrans-
formed KdV equation.
In discussing implicit methods, we have stressed repeatedly that their usefulness is
SITUATION-DEPENDENT. Implicit methods are good for computing KdV solitons because
they misrepresent only the high frequency transient while accurately computing the slow-
moving solitary waves. Integrating factor methods, in contrast, are bad for solitons because
the integrating factor extracts a high frequency motion that is not actually present in the
low frequency motion (i. e., the soliton) that we are computing. Thus, the usefulness
of integrating factor methods is also situation-dependent, but often inversely to implicit
Suppose our goal is to compute KdV solutions which are not solitons but rather a soup
of weakly interacting waves such that the Fourier components are approximately

ak = Ak ( t) exp(ik 3 t) (14.28)

where the Ak are amplitudes that vary slowly under weak wave-wave interactions. Now
because the integrating factor method successfully removes the high frequency oscillations,
the transformed coef¬cients are

vk (t) = Ak ( t)

These vary slowly with time, allowing a long time step. In contrast, implicit and Nonlinear
Galerkin methods as well as standard explicit schemes all require a very short time step
to be accurate for the high frequency motion. Weakly coupled wave systems (“resonant
triads” and their generalization) have been widely analyzed by theorists (Craik, 1985).
Thus, it is quite wrong to conclude that integrating factor methods are counter-productive
in numerical applications. Sometimes they are very valuable. The difference between the
two KdV cases is that in the ¬rst (soliton) case, the integrating factor extracts a high fre-
quency dependence that was not actually present in the slow solutions (bad!) whereas in
the second (weakly interacting sine waves) case, the integrating factor extracts high fre-
quency motion that is an essential part of the desired solution (good!).

14.3 Semi-Lagrangian Advective Schemes: Introduction
A “semi-Lagrangian” (SL) time-marching algorithm is an integrating-factor method in
which the integrating factor is an advection operator that shifts into a coordinate system
that travels with the ¬‚uid. In so doing, it borrows ideas from both the method of char-
acteristics and the Lagrangian coordinate formulation of ¬‚uid mechanics. The “total” or
“particle-following” derivative of any quantity “q”, for simplicity given in two space di-
mensions, is
Dq ‚q ‚q ‚q
≡ (14.30)
+u +v
Dt ‚t ‚x ‚y

In the moving coordinate system, the total derivative becomes the ordinary partial time
derivative, and the advective terms disappear, buried in the change-of-coordinate.
The idea of tracing ¬‚ow by following the motion of individual blobs or particles of
¬‚uid is at least two centuries old. Particle-following coordinates are called “Lagrangian”
coordinates in ¬‚uid mechanics. However, exclusive use of a Lagrangian coordinate system
is usually a disaster-wrapped-in-a-catastrophe for numerical schemes because the parti-
cle trajectories become chaotic and wildly mixed in a short period of time where “chaos”
is meant literally in the sense of dynamical systems theory even for many laminar, non-
turbulent ¬‚ows.1 Semi-Lagrangian (SL) algorithms avoid this problem (and earn the mod-
i¬er “semi”) by reinitializing the Lagrangian coordinate system after each time step. The
usual choice is that the tracked particles are those which arrive at the points of a regular
Fourier or Chebyshev spatial grid at time tn+1 . Since the departure points for these trajec-
tories at t = tn are generally not on the regular grid, we follow a different set of N particles
on each time interval of length „ where „ is the time step. In contrast, in a pure Lagrangian
coordinate system, one follows the evolution of the same N particles forever.
There are two advantages to an SL method:
1 Trajectories of ¬‚uid blobs are often chaotic even when the advecting velocity ¬eld is laminar.

t4 t4
t1 t1
t0 t0
x x
Figure 14.1: Left panel: Fluid mechanics in Lagrangian coordinates tracks the same set
of particles throughout all time. Right panel: Semi-Lagrangian methods. A new set of
particles to track is chosen at the beginning of each time step.

1. If the scheme is “multiply upstream”, a term de¬ned later, SL will be both stable
and accurate for time steps far greater than usual explicit limit. To be precise, the
advective time step limit disappears, and the maximum time step is then controlled
entirely by non-advective processes like viscosity and wave propagation.
2. SL methods are very good at resolving fronts and shock waves even in the absence
of arti¬cial viscosity.
Because of these virtues, SL time-marching is now employed in the NCAR Community
Climate Model and the European Centre for Medium-Range Forecasting weather predic-
tion model. Neither is a toy model, as too often used to test algorithms; instead, both are
global, three-dimensional spherical harmonic models with lots of physics.
Simple semi-Lagrangian algorithms have at least ¬ve disadvantages:
1. Dif¬cult to parallelize (for “multiply upstream” schemes).
2. Not strictly conservative of mass, energy, etc.
3. Inherent computational dissipation (McCalpin, 1988)
4. Instability and inaccuracy for ¬‚ows over mountains.
5. More expensive per time step than most competing schemes because of the need for
• “off-grid” interpolation and
• ¬xed point iteration.
None of these drawbacks, however, has proved fatal. In particular, conservative, shape-
preserving and low dissipation schemes have been recently developed, and the mountain
¬‚ow or “orographic” problems have been solved, too.
Semi-Lagrangian method are designed to be especially good at advection. However,
the algorithms can be applied to the full Navier-Stokes equations without dif¬culty and
can be combined with other ideas. For example, in meteorology it is common to combine
semi-Lagrangian treatment of advection with a semi-implicit treatment of the viscous and
Coriolis force terms (“SL/SI” schemes).
Semi-Lagrangian time-stepping can be combined with any strategy for discretization of
the spatial coordinates.

14.4 A Brief Digression on the One-Dimensional and Two-
Dimensional Advection Equations and the Method of
The Navier-Stokes equations reduce to the so-called “advection equations” if we neglect
all processes except advection. Restricting ourselves to two-dimensions for simplicity, the
pure advection equations are
Du Dv
=0 & =0
Dt Dt
or written in full, where (u, v) are the x and y velocities,
‚u ‚u ‚u
+u +v = 0
‚t ‚x ‚y
‚v ‚v ‚v
+u +v = 0
‚t ‚x ‚y

Figure 14.2: Contour plot in the space-time plane of the solution u(x, t) to the one-
dimensional advection equation, ut + uux = 0. The initial condition is sinusoidal; negative
contours are dashed. Where contours intersect, the solution becomes multi-valued and
ux ’ ∞. The ¬‚ow is said to have “formed a front” or experienced “wave breaking”.

The reason that this drastically simpli¬ed system is interesting is that it can be solved ex-
actly by the “method of characteristics”. Semi-Lagrangian time-marching schemes, even
with all the complications of pressure and viscosity, etc., merely mimic the method of char-
The key observation is that the advection equations state that the velocities are con-
served following the motion. In other words, if we track the fate of a particular blob of
¬‚uid, perhaps marked with a dab of red dye, the velocity of that blob will remain constant
forever. Consequently, the trajectory of that blob and of all blobs is a straight line whose
slope in each spatial dimension is proportional to its velocity with respect to that coordi-
nate. (Fig. 14.2 is a one-dimensional illustration.) If the blob originates at x = ξ and y = ·
at t = 0, then its trajectory is

x = ξ + Q(ξ, ·) t, & y = · + S(ξ, ·) t

where the initial conditions are

u(x, y, t = 0) = Q(x, y), & v(x, y, t = 0) = S(x, y)

The exact solution for the velocities is

u(x, y, t) = Q(ξ[x, y, t], ·[x, y, t]), & v(x, y, t) = S(ξ[x, y, t], ·[x, y, t])

In other words, the initial conditions furnish the exact solution for all time if the arguments
of the functions (Q, S) are replaced by the “characteristic coordinates” or simply “charac-
teristics”, (ξ, ·).

Unfortunately, Eq.(14.33) gives explicit solutions only for (x, y) as a functions of the
characteristics. To evaluate the solutions at a particular point (x, y, t), we need it the other
way around. Thus, it is necessary to solve (14.33) to determine the departure point of a
given blob where the “departure point“ or “foot of the trajectory” is the location (ξ, ·) of
the blob at t = 0. However, the equations that must be solved for the characteristics are
algebraic equations, which are much simpler than differential equations.
Semi-Lagrangian time-marching schemes also require computing the departure point
of a blob by solving a system of algebraic equations. However, this is fairly easy because
the blob moves only a short distance during a time interval of length „ where „ is the time
step. The change of coordinates (14.33) does not trivialize the dynamics of the full Navier-
Stokes equations, as it does for the pure advection equations, but it does remove advection
just as thoroughly.

14.5 The Simplest SL Scheme: Three-Level, Second Order
Semi-Implicit Method
A semi-Lagrangian scheme by itself is not very useful in meteorology because it only re-
moves the advective time step limit, which is much less severe than the restriction imposed
by the propagation of gravity waves. When combined with a semi-implicit treatment of the
non-advective terms, however, SL becomes very powerful as ¬rst demonstrated by Robert
(1981, 1982).
Let the problem be
+ G(u, t) = F (u, t)
where D/Dt denotes the total, particle-following derivative including the advective terms,
G is the sum of all terms which are to be treated implicitly while F represents all the terms
which can be treated explicitly. The basic three-level, second order SL/SI algorithm is

u + ’ u’ 1
+ (G+ + G’ ) = F 0 (14.37)
2„ 2
where the plus and minus signs denote the quantities evaluated at t = tn±1 and super-
script 0 denotes a quantity evaluated at t = tn . In form, this looks like a standard Crank-
Nicholson/explicit method except that (i) the advective terms are missing (ii) the CN-like
step is over a time interval of 2„ instead of „ where „ is the time step.
The magic is in the coordinate transform. For simplicity, assume one space dimension.
Let ±(x) denote the displacement, over the time interval t ∈ [tn , tn+1 ] of a particle that
arrives at point x at t = tn+1 . Then for any vector appearing in (14.37),

≡ u(x, tn + „ )
≡ u(x ’ ±, tn )
u’ ≡ u(x ’ 2 ±, tn ’ „ )

When the spatial dependence is discretized, u, x, and ±(x) all turn into vectors. The dis-
placement is de¬ned to be the solution of the algebraic equation

±j = „ u(xj ’ ±j , tn ) (14.39)

where j is the integer index of the grid points.

The novelty of Eq.(14.39) is that when we move back one time level, we must simulta-
neously shift our spatial position by the displacement ±j . The predicted value u(xj , tn+1 )
is combined not with u at the same spatial location at the two previous time levels, but


. 57
( 136 .)