. 56
( 136 .)


un Advection

Waves & Pressure

un+2/3 Viscosity/Diffusion

Figure 13.4: Splitting by process: Advection, adjustment of ¬elds (linear wave dynamics
and pressure), and viscosity and diffusion are split into three separate, one-after-the-other

’ f v n+2/3 = ’φx
+ f un+2/3 = ’φy
[Pressure Step alias (13.49b)
n+2/3 n+2/3 n+2/3
“Adjustment of Fields”] (13.49c)
φt + ux + vy =0

un+1 = µ un+1 (13.50a)
v n+1 [Diffusion Step] (13.50b)
vt =µ
φn+1 = φn+2/3 (13.50c)
In each case, the solution at the previous fractional step provides the initial condition for
the equations that model the next physical process.
If we wish, we can apply an additional splitting to decompose the diffusion step into
a sequence of one-dimensional diffusion problems. However, it is usually a bad idea to

attempt dimensional splitting for the “adjustment of ¬elds” step because the Coriolis force
and wave dynamics strongly couple the coordinates. Splitting advection into separate x,
y and z advection is unnecessary because this subproblem is normally integrated by an
explicit method.
As noted earlier, the great advantage of splitting is that we can solve each fractional
step using a different numerical method. The advective step is almost always advanced
explicitly and the adjustment of ¬elds is integrated implicitly. In engineering problems,
diffusion is usually marched in time by an implicit algorithm, but Marchuk (1974) barely
discusses diffusion because µ in his numerical weather prediction model is so small that
he can use an explicit method.
It is not even necessary to be consistent by using spectral or pseudospectral algorithms
for all sub-problems. Lambiotte et al. (1982) compute only the advection and pressure steps
pseudospectrally. Since µ is small for their high Reynolds ¬‚ow, they employ ¬nite differences
to integrate the diffusive step.
Chapter 14

Semi-Lagrangian Advection &
Other Integrating Factor Methods

“Go with the ¬‚ow”

“ New Age proverb

“... there is a mistaken belief that semi-Lagrangian schemes are only good for smooth

“ Staniforth & Cot´ (1991), pg. 2211

14.1 Introduction: Concept of an Integrating Factor
An “integrating factor” is a transformation of the unknowns in a differential equation. The
most familiar example yields the solution of the general linear ¬rst order ODE, which is

ux ’ q(x) u = f (x) (14.1)

Let u(x) be replaced by the new unknown v(x) where

u(x) ≡ I(x) v(x) (14.2)

where the integrating factor is
I(x) ≡ exp (14.3)
q(y) dy

Substitution gives the transformed problem

I vx + q(x) I v ’ q(x) I v = f (x) (14.4)

using Ix = q(x)I. Dividing by I(x) gives
’ (14.5)
vx = exp q(y) dy f (x)


Thus, the integrating factor has reduced the differential equation to a quadrature, i. e., v(x)
is the inde¬nite integral of the right-hand side of (14.5):
x z
’ (14.6)
v =A+ exp q(y) dy f (z) dz

where A = u(0).
Integrating factors come in all sorts of ¬‚avors, like ice cream, but we shall concentrate
on two particular classes. The ¬rst is when a time-dependent PDE can be written in the
ut = L u + M (u)
where L is a linear differential operator with spatial derivatives only and where M (u) is
an operator, possibly nonlinear, that contains everything else in the equation including the
forcing terms. The new unknown de¬ned by
u ≡ exp(Lt) v (14.8)
solves the simpli¬ed equation
vt = exp(’Lt) M (v)
Just as for the ¬rst order ODE, the integrating factor has purged the transformed equation
of the differential operator L except through the exponential. This is sometimes a good
idea ” and sometimes not.
The second class of integrating factors is to remove translation/advection by a change-
of-coordinates to a moving coordinate system. Formally, this change-of-coordinate can
be expressed in terms of operators, too. However, it is more convenient to perform the
change-of-coordinate directly.
The simplest case is to shift into a frame of reference which is moving at a constant
speed with respect to the original frame. For example,
ut + cux = f (x, t), u(0, t) = g(t), u(x, 0) = Q(x)
can be transformed by the shift
ξ ≡x’ct ” u(x, t) ≡ v(ξ, t) (14.11)
x=ξ+ct &
vt = f (ξ + ct, t), v(’ct, t) = g(t), v(ξ, 0) = Q(ξ)
For the special case f = 0, this is clearly an improvement on the original solution because
then vt = 0 which implies that v is a function of ξ only, that is, u depends on (x, t) only as
the combination x ’ ct.
Semi-Lagrangian time-marching schemes remove nonlinear advection through a sim-
ilar but more complicated transformation. Such algorithms have so many good proper-
ties that they have become standard on two important geophysical models: the NCAR
Community Climate Model (CCM) and the operational forecasting code of the European
Centre for Medium Range Forecasting, which provides ¬ve-day global predications to an
18-nation consortium.
However, integrating factor schemes can also be a trap. The crucial distinction is that
integrating factors are useful only when the integrating factor removes an important part
of the dynamics. We will explain this through a couple of bad examples in the next section.
(This section can be skipped without loss of continuity.)
In the remainder of the chapter, we describe semi-Lagrangian schemes. We begin with
second order schemes, which are the easiest to describe, and then explain a number of

14.2 Misuse of Integrating Factor Methods

Our ¬rst counterexample is the ODE which was used earlier to explain the concept of
the slow manifold:

ut + »u = f ( t)

Assume »/ >> 1, i. e., that the ratio of the “fast” time scale to the slow time scale is
large, and also that the real part of » is positive so that the homogeneous solution decays
exponentially with time. We shall further assume that the goal is to compute the slow
manifold of this equation, that is, the part which varies only on the same slow O(1/ ) time
scale as the forcing f .
Unfortunately, the time step for an explicit scheme is restricted by the fast time scale of
the homogeneous solution exp(’»t). Implicit methods and Nonlinear Galerkin methods
both remove this restriction, but have their own disadvantages. Is an integrating factor a
viable alternative?
De¬ne a new unknown by

u = exp(’» t)v(t)

The ODE becomes

vt = exp(» t)f ( t)

Hurrah! The term »u has been removed and with it the usual CFL time-step limitation. To
show this, let us analyze the Forward Euler method as applied to (14.15).
This scheme with time step „ and with superscript n denoting the time level is

v n+1 ’ v n
” v n+1 = v n + „ exp(» n „ )f ( n „ ) (14.16)
= exp(» n „ )f ( n „ )

From an arbitrary starting value v 0 and for the special case f ( t) = F , a constant,

= v 0 + „ F {1 + exp(» „ ) + exp(2 » „ ) + . . . + exp([n ’ 1] » „ )}
1 ’ exp(n » „ )
= v0 + „ F (14.17)
1 ’ exp(» „ )

Because of the integrating factor transformation, v n is stable for arbitrarily large time steps
„ . To be sure, the solution is blowing up exponentially, but since the transformed forcing
is growing exponentially, a similar growth in v is okay.
The same solution in terms of the original unknown u is, using v 0 ≡ u0 ,

exp(’n » „ ) ’ 1
un = exp(’n »„ )u0 + „ F (14.18)
1 ’ exp(» „ )

In the limit that „ ’ 0 and n ’ ∞ such that t = n „ grows large, the Forward Eu-
ler/integrating factor solution simpli¬es to


un „F
1 ’ (1 + » „ + O(»2 „ 2 ))
≈ {1 + O(»„ )} (14.19)

where F/» is the slow manifold for this problem. Thus, for suf¬ciently small time step, the
scheme does successfully recover the slow part of the solution as t ’ ∞.
However, the whole point of the integrating factor is to use a long time step without the
bother of implicit or Nonlinear Galerkin methods. Since the stability limit of the Forward
Euler scheme without integrating factor is „max = 1/», it is interesting to see how well the
integrating factor variant does for this „ . Unfortunately, substituting „ = „max into (14.18)
F 1 F
un ∼ t ’ ∞, „ = 1/» (14.20)
= 0.58 ,
» ’1 + exp(1) »
The good news is that the large time solution is proportional to the true slow manifold,
F/»; the bad news is that the amplitude of the computational solution is barely half that of
the true solution even when „ is no larger than the largest time step that would be stable
in the absence of the integrating factor. When „ >> 1/», which is of course what we really
want if the integrating factor scheme is to compete with implicit schemes,
un ∼ t ’ ∞, »„ >> 1 (14.21)
» „ exp(’» „ ),
which is exponentially small compared to the correct answer.
Thus, the integrating factor method is a dismal failure. If the time step „ = r „max
where „max = 1/» is the stability limit for the Forward Euler method (without integrating
factor), then the Forward Euler/integrating factor algorithm gives a slow manifold which
is too small by a factor of r exp(’r). In contrast, implicit methods like Crank-Nicholson
and Backwards-Euler and also Nonlinear Galerkin schemes are accurate for the slow part
of the solution even for large r.
A second bad example is the Fourier pseudospectral solution to the KdV equation
ut + uux + uxxx = 0
when the goal is to model the formation and dynamics of solitary waves. For a soliton of
phase speed c, the coef¬cient of exp(ikx) varies with time as exp(’ikct) which grows only
linearly with k. However, the equation also admits small amplitude waves whose phase
speeds grows as k 3 , so the stability limit for an explicit time-marching scheme is inversely
proportional to kmax . The KdV equation is very “stiff” in the sense that one must use a
very short time step for stability reasons compared to the slow time scale of the dynamics
(i. e., the solitons) that one is actually interested in.
The remedy ” well, sort of a remedy ” is to use the integrating factor exp(’ ‚xxx t) to
remove the third derivative from the transformed equation:
vt + v exp(’‚xxx t) vx = 0


. 56
( 136 .)