„

’ 4un + 3un’1 ’ 4 un’2 + 1 un’3

25 n+1

12 u

= F (un+1 , tn+1 ) [BDF4]

3 4

(12.27)

„

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 )

n+1

[CN] (12.30)

u = u+

1 ’ „ »/2 1 ’ „ »/2

2

12.5. SEMI-IMPLICIT METHODS 229

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

log

µ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

(12.32)

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)

+

2

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

CHAPTER 12. IMPLICIT SCHEMES & THE SLOW MANIFOLD

230

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

acceptable.

Semi-implicit algorithms have therefore become almost universal in weather forecast-

ing and climate models.

12.6 Speed-Reduction Rule-of-Thumb

Rule-of-Thumb 12 (IMPLICIT SCHEME FORCES PHYSICS SLOWDOWN)

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.

Rule-of-Thumb 13 (EXPLICIT-DEMANDING)

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)

12.7. SLOW MANIFOLD: METEOROLOGY 231

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

q

„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

scheme.

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

CHAPTER 12. IMPLICIT SCHEMES & THE SLOW MANIFOLD

232

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

meteorology.

“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.

EXAMPLE ONE: “TRIVIAL” Slow Manifold: STEADY-STATE

12.8. SLOW MANIFOLD: DEFINITION & EXAMPLES 233

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.

EXAMPLE TWO: FORCED LINEAR OSCILLATOR

(12.39)

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