. 49
( 136 .)


= F (un+1 , tn+1 ) [BDF3] (12.26)

’ 4un + 3un’1 ’ 4 un’2 + 1 un’3
25 n+1
12 u
= F (un+1 , tn+1 ) [BDF4]
3 4

un+1 ’ un F (un+1 , tn+1 ) + F (un , tn+1 )
[Crank-Nicholson (CN)] (12.28)
„ 2
For the special case of the linear, scalar-valued ODE, du/dt = »u + f (t), it is easy to
solve the implicit formulas:

un f (tn+1 )
un+1 = [BE] (12.29)
1 ’ „» 1 ’ „»

1 + „ »/2 n „ f (tn+1 ) + f (tn )
[CN] (12.30)
u = u+
1 ’ „ »/2 1 ’ „ »/2

If » is negative and real, then the solutions to the difference equations decay geometrically
to the steady state as exp(µt) where

µBE log(1 + „ »)
» „»
1+„ »/2
µCN 1’„ »/2
=’ (12.31)
» „»
Since 1 + „ » must always be smaller than unity (since » < 0), it follows that the logarithms
in (12.31) are always negative so that the ratios of µ/» are always positive. In other words,
the numerical solutions decay, as does the true solution, for all „ ∈ [0, ∞]. The conditional
stability of the Forward Euler scheme (where the condition is „ < 1/|»|) has been replaced,
for these two implicit schemes, by unconditional stability.
It is left to the reader to demonstrate “ if not already convinced by earlier numerical
analysis texts or courses “ that the unconditional stability extends to » anywhere in the left
half-plane including the imaginary axes, which is called the property of being “A-stable”.
Higher order schemes are not so well-behaved; BDF 3 and BDF 4, for example, are un-
stable for small imaginary ». However, these algorithms are still sometimes used because
weak dissipation ( (») < 0) is suf¬cient to stabilize them. Fornberg(1996, Appendix G)
gives stability diagrams for a wide variety of schemes.

12.5 Semi-Implicit Methods
For nonlinear equations, an implicit algorithm has a high cost-per-timestep because one
must solve a nonlinear boundary value problem at every time step. For the Navier-Stokes
equations and for numerical weather forecasting and climate modelling, it is therefore com-
mon to use an algorithm in the following class.

De¬nition 26 (SEMI-IMPLICIT) If some terms in a differential equation are approximated by
an explicit scheme while others are approximated implicitly, then the resulting time-marching algo-
rithm is said to be SEMI-IMPLICIT.
For example, for the nonlinear PDE

ut = F (u, x, t) + L(u, x, t)

where F and L denote the nonlinear and linear parts of the differential equation, a popular semi-
implicit scheme is (“[AB3CN]”)

23 4 5
F (un , tn ) ’ F un’1 , tn’1 +
un+1 = un + „ F un’2 , tn’2
12 3 12

L(un+1 , tn+1 ) + L(un , tn ) (12.33)
In general, the implicit and explicit schemes may be of the same or different order.

It may seem illogical to treat some terms explicitly while other terms implicitly, and
even more bizarre to use schemes of different orders, as in the AB3CN™s mix of third order
Adams-Bashforth with the Crank-Nicholson scheme, which is only second order. How-
ever, there are major advantages.
First, because the nonlinear term is treated explicitly, it is only necessary to solve a lin-
ear boundary value problem at each time step. Second, the viscous terms, which involve

second derivatives, impose much stiffer time step requirements („ ∼ O(1/N 4 ) for a Cheby-
shev basis) than the advective terms, which involve only ¬rst derivatives and impose a
timestep limit proportional to 1/N 2 . In other words, the semi-implicit algorithm stabilizes
the most dangerous terms. Third, in weather forecasting and other ¬‚uid dynamics, ad-
vection is crucial. It is important to use a high order time-marching scheme with a short
timestep to accurately compute frontogenesis, turbulent cascades, advection of storm sys-
tems and so on. There is little advantage to treating the nonlinear terms implicitly because
a timestep longer than the explicit advective stability limit would be too inaccurate to be
Semi-implicit algorithms have therefore become almost universal in weather forecast-
ing and climate models.

12.6 Speed-Reduction Rule-of-Thumb
Implicit schemes obtain their stability by slowing down the time evolution of the numerical
solution so that explicit stability criteria are satis¬ed. In other words, if the implicit algorithm
approximates the true eigenvalue » by µ where the homogeneous solution evolves as exp(»t), then

„ | µ(„ ) | < O(1) (12.34)

where „ is the time step. (Recall that the stability condition for the Forward Euler method is „ |»| < 1
and for RK4 is „ » < 2.8.)

The Rule-of-Thumb has not been demonstrated for all possible implicit algorithms, but
no counterexamples have been published. The Rule-of-Thumb is also not quite a theorem
because it ignores the distinction between logarithmic growth and a constant.
For example, it is easy to show from (12.31) that

„ |µBE | = | log(1 + „ »)|
1 + „ »/2
„ |µCN | = log (12.35)
1 ’ „ »/2

Thus, the slowed-down numerical eigenvalue µ does not quite satisfy a stability require-
ment of the form

„ |µ| < q (12.36)

for a constant q for all „ but only for q which is allowed to grow logarithmically with the
time step. However, when „ |»| = 1000, that is, when the time step is 1000 times larger than
would be stable for the Forward Euler method, q < 8, that is, the homogeneous solution to
du/dt = »u + f (t) in the implicit difference schemes are evolving roughly 8/1000 times as
fast as for the exact solution.
Obviously, only a lunatic would use so long a time step if tracking evolution on the
time scale of 1/» is physically signi¬cant. This can be formalized as the following.

Suppose that the stability requirement for an explicit time-marching scheme for a time-dependent
PDE is, for some constant q,

„ |»max | < q (12.37)

where, after discretization using N degrees of freedom for the spatial derivatives, »max (N ) is the
eigenvalue of the discretization matrix which is largest in magnitude. If it is important to follow
wave propagation, damping or other time behavior on the time scale of exp(i»max t), then one
should use an EXPLICIT scheme with a time step which is SMALL compared to
„limit ≡ (12.38)
|»max |

The ¬rst point in justifying this Rule-of-Thumb is the observation that if „ ∼ „limit ,
that is, if the time step is near the stability limit, then any explicit or implicit algorithm
will have large errors in following decay or propagation on a time scale of O(1/|»max |). (It
may be highly accurate for slower components of the solution, but the Rule only applies
when these slow components are not enough, and one needs to accurately resolve the
fastest component, too.) It follows that we must use a time step which is small compared to
the stability limit to accurately track the dynamics of the fast component that imposes the
stability limit.
Given that the time step is so small, we may use either an implicit or explicit method.
Since implicit methods are usually much more expensive than explicit schemes for the
same time step, however, the second point is that it is usually cheaper to use an explicit
There are some exceptions. First, implicit methods require solving a Boundary Value
Problem (BVP) at every time step, but if the basis is Fourier or spherical harmonics and
the operator of the BVP is constant coef¬cient or a Laplace operator, solving the BVP may
be trivial. In that case, the Crank-Nicholson scheme, which is second order with a smaller
proportionality constant than most explicit schemes of the same order, is quite attractive.
Second, a Chebyshev basis usually generates some eigenvalues which are poor approxi-
mations to those of the corresponding differential operator. For the Laplacian, for example,
some of the Chebyshev eigenvalues are O(N 4 ), implying a time step which is no larger than
q/N 4 for an explicit scheme, while the ¬rst N true eigenvalues are bounded by N 2 . In this
case, stability is limited by garbage modes, that is, by large spurious eigenvalues whose
magnitude is unrelated to anything physical. It is sensible to use an implicit method with
a time step which is O(1/N 2 ) because then all the modes whose eigenvalues are good ap-
proximations to those of the differential operator will be integrated accurately in time; the
modes which are advanced with poor time accuracy are nonsense anyway.
Even so, it is clear that implicit methods must be used with discretion. In meteorology,
where implicit methods have been wisely employed for many years, this has led to the
concept of the “slow manifold”.

12.7 Slow Manifold: Meteorology
Weather forecasting is done by solving the “primitive equations”, which ¬lter sound waves
so that only two important classes of waves remain. Rossby waves (and also nonlinear
advection, which has roughly the same time scale) are responsible for the travelling high-
pressure and low-pressure centers that are the large-scale weather. Gravity waves are gen-
erated by small-scale phenomena like thunderstorms and have frequencies an order of
magnitude higher than Rossby waves or the advective time scale.
Observations supply initial conditions for the numerical code. Rossby waves are well-
resolved by the radiosonde and satellite network, but gravity waves are not. However,
it has been known for many years that the amplitude of gravity waves, when averaged
over the globe, is very small compared to Rossby waves. If we imagine a con¬guration

space of 3N dimensions where each coordinate is the amplitude of a wave mode in a nu-
merical model, then the instantaneous state of the atmosphere should lie close to the N -
dimensional subspace of pure Rossby wave motion with all the gravitational modes equal
to zero.
It follows that a good global forecasting model does not need to accurately track move-
ments in the whole phase space, but only ¬‚ow on the “slow manifold”. Implicit methods
are natural tools to integrate on the slow manifold. Their stability allows a long time step.
On the slow manifold “ but only on the slow manifold “ the long time step does not com-
promise the accuracy of the numerical solution.
Gravity waves limit the timestep for explicit algorithms to less than ten minutes. How-
ever, with an implicit or semi-implicit algorithm, one can make good ¬ve-day forecasts
with a timestep of an hour. Because of the errors introduced by subgridscale turbulence,
cumulus convection, photochemistry and so on, the total forecast error cannot be signi¬-
cantly reduced by using a shorter timestep.
Because the advective limit is roughly an hour and a timestep longer than an hour
would damage accuracy anyway, it is usual in forecasting to treat the nonlinear terms ex-
plicitly. The linear terms associated with gravity waves are treated implicitly, but in a
spherical harmonics basis, the (linear) BVP is easy to solve as explained in Chapter 18. The
semi-implicit algorithm costs only a little more per timestep than an explicit scheme, but
permits a time step six times longer.
Because the goal is to forecast only slow motion, implicit methods make sense for global
“Global” means a model whose goal is to forecast weather over a large area such as the
continental United States. In small areas (such as the rain shadow of a thunderhead), local
“fast” dynamics such as thunderstorms may be very important. (In densely-populated
areas, local Doppler radar and high-density instrument arrays can provide enough “fast”
data to make mesoscale modelling feasible locally, though not globally.) It follows that
so-called “mesoscale” models, which attempt to forecast severe weather for a few hours
ahead, must resolve the fast motion, not just the slow manifold, and therefore usually
employ explicit time-marching.
The concept of a “slow manifold” and the usefulness of implicit time-stepping methods
are situational.

12.8 Slow Manifold: General De¬nition and
Physical Examples
Many scienti¬c problems have multiple time scales. A major branch of perturbation the-
ory, the “method of multiple scales”, exploits these different time scales. Implicit time-
marching schemes are a numerical strategy for doing the same thing.
As explained in the previous section, however, implicit methods are only sensible when
the goal is to track motion on the “slow manifold”.

De¬nition 27 (SLOW MANIFOLD)
In the phase space of a numerical model which allows motion on both a fast time scale and a slow
time scale, the SLOW MANIFOLD is the submanifold where the time evolution is on the longer
time scale only.


If a partial differential equation (and its discretization) have a steady state, then this is
the ultimate slow manifold. If we wish to integrate in time to allow the ¬‚ow to ¬nd its own
way to the steady state, we may use inaccurate time-marching schemes because only the
¬nal steady state matters.


du/dt + iu = f ( t), << 1

The homogeneous solution of the oscillator, exp(’it), is the “fast” motion with a time scale
of unity; the O(1/ ) time scale of the forcing is the “slow” scale. The “slow manifold” is
that particular solution which varies only on the slow time scale; the general solution to
the ODE contains fast oscillations also.
Although the forced, linear oscillator seems too simple to be useful, it contains hidden


. 49
( 136 .)