FOURTH-ORDER TIME STEPPING FOR STIFF PDES
∗
ALY-KHAN KASSAM
†
AND LLOYD N. TREFETHEN
‡
Abstract. A modi fication of the ETDRK4 (Exponential Time Differencing fourth-order Runge-
Kutta) method for solving stiff nonlinear PDEs is presented that solves the problem of numerical
instability in the scheme as proposed by Cox and Matthews and generalizes the method to non-
diagonal operators. A comparison is made of the performance of this modified ETD s cheme against
the competing methods of implicit–explicit differencing, integrating factors, time-splitting, and Forn-
berg and Driscoll’s “sliders” for the KdV, Kuramoto-Sivashinsky, Burgers, and Allen-Cahn equations
in one space dimension. Implementation of the method is illustrated by short Matlab programs for
two of the equations. It is found that for these applications with fixed time steps, the modified ETD
scheme is the best.
Key words. ET D, exponential time- differencing, KdV, Kuramoto-Sivashinsky, Burgers, Allen-
Cahn, implicit–explicit, split step, integrating factor
AMS subject classifications.
1. Introduction. Many time-dependent partial differential equations (PDE s)
combine low-order nonlinear terms with higher-order linear terms. Ex amples include
the Allen-Cahn, Burgers, Cahn-Hilliard, Fisher-KPP, Fitzhugh-Naguno, Gray-Scott,
Hodgkin-Huxley, Kuramoto- Sivashinsky, Navier-Sto kes, nonlinear Schr¨odinger, and
Swift-Hohenberg equations. To obtain acc urate numerical solutions of such prob-
lems, it is desirable to use high-order approximations in space and time. Yet because
of the difficulties introduced by the combination of nonlinearity and stiffness, most
computations here tofore have been limited to second order in time.
Our subject in this paper is fourth o rder time-differencing. We shall write the
PDE in the form
u
t
= Lu + N (u, t),(1.1)
where L and N ar e linear and nonlinea r operators, respectively. Once we discretise
the spatial part of the PDE we get a system of ODEs,
u
t
= Lu + N(u, t).(1.2)
There seem to be five principal competing methods for solving problems of this kind,
which we will abbreviate by IMEX, SS, IF, SL, and E TD. Of course these are not
the only schemes that are being used. Noteworthy schemes that we ignore ar e the
exp onential Runge–Kutta schemes [26] and deferred correction [37] or semi–implicit
deferred correction [6, 45].
IMEX = Implicit–explicit. These are a well studied family of schemes that have
an established history in the solution of stiff PDE. Early work loo king at some
stability issues dates to the beginning of the 1980’s [61]. Schemes have been prop osed
for specific examples, e.g. the KdV equatio n [13] and the Navier-Stokes equations
[11, 32, 34], as well as certain classes of problems, for example re action-diffusion
∗
This work was supported by the Engineering and Physical Sciences Research Council (UK) and
by MathWorks, Inc.
†
Oxford University Computing Laboratory, Wolfson Bldg., Parks Road, Oxford OX1 3QD, UK
(akk@comlab.ox.ac.uk)
‡
same address (LNT@comlab.ox.ac.uk).
1
2 A.-K. KASSAM AND L. N. TREFETHEN
problems [51] and atmospheric modelling problems [62]. An overview of the stability
properties and derivations of implicit–ex plicit schemes can be found in [2].
Implicit–explicit schemes consist of using an explicit multi-step formula, for ex-
ample the second order Adams–Bashforth formula, to advance the nonlinear part of
the problem and an implicit scheme, for example the second order Adams–Moulton
formula, to advance the linear part. Other kinds of formulations also exist; for de-
velopments based around Runge-Kutta rather tha n Adams-Bashforth formulae, for
example, see again work by Ascher et al. [3] as well as very recent work by Ca lvo et al.
[10] and Kennedy and Carpenter [33]. In this report, we use a scheme known either
as AB4BD4 (in [14]) or SBDF4 (in [2]), which consists of combining a fourth order
Adams–Bashforth and a fourth order backward differentiation scheme. The formula
for this scheme is
u
n+1
= (25 − 12hL)
−1
( 48u
n
− 36u
n−1
+ 16u
n−2
− 3u
n−3
+ 48hN
n
(1.3)
− 72N
n−1
+ 48N
n−2
− 12N
n−3
) (AB4BD4).
SS = Split Step. The idea of split step methods seems to have originated with
Bagrinovskii and Godunov in the late 1950’s [4] and to have been independently
developed by Strang for the construction of finite difference schemes [57] (the simplest
of these is often called ‘Strang Splitting’). The idea has been widely used in modelling
Hamiltonian dynamics, with the Ha miltonian of a system split into its potential and
kinetic energy pa rts. Some early work on this was done by Ruth [50]. Yoshida
[63] developed a technique to produce split-step methods of arbitrary even order.
McLachlan and Atela [41] studied the accuracy of such schemes and McLachlan [42]
made some further comparisons of different symplectic and non-symplectic schemes.
Overviews of these methods can be found in Sanz-Serna and Calvo [53] and Boyd [7],
and a recent discussion of the relative merits of operator splitting in general can be
found in a paper by Schatzman [54].
In essence, with the split step method, we want to write the solution as a com-
position of linear and nonlinear steps
u(t) ≈ exp(c
1
tL)F (d
1
tN) exp(c
2
tL)F (d
2
tN)· · · exp(c
k
tL)F (d
k
tN)u(0),(1.4)
where c
i
and d
i
are rea l numbers and represent fractional time steps (though we
use product notation, the nonlinea r substeps are nonlinea r). Generating split step
methods becomes a process of generating the appropriate sets of real numbers, {c
i
}
and {d
i
}, such that this product matches the exact evolution operator to high order.
The time–stepping for such a scheme can be e ither a multistep or a Runge-Kutta
formula. We use a fourth order Runge-Kutta formula for the time–stepping in this
exp eriment.
IF = Integrating factor. Techniques that multiply both sides of a differential
equation by some integrating factor and then make a relevant change of variable a re
well known in the theory of ODEs (see for example [38]). A similar method has been
developed for the study of PDEs. The idea is to make a change of variable that allows
us to solve for the linear part exactly, and then use a numerical scheme of our choosing
to solve the transformed, nonlinear equation. This technique has been used for PDEs
by Milewski and Tabak [44], Maday et al. [40], Smith and Waleffe [55, 56], Fornberg
FOURTH-ORDER TIME-STEPPING FOR STIFF PDE 3
and Driscoll [20], Trefethen [60], Boyd [7] and Cox and Matthews [14]. Starting with
our generic discretised PDE, we define
v = e
−Lt
u.(1.5)
The term e
−Lt
is known as the integrating factor. In many applications we can work in
Fourier space and render L diagonal, so that scalars rather than matrices ar e involved.
Differentiating (1.5) gives
v
t
= −e
−Lt
Lu + e
−Lt
u
t
.(1.6)
Now, multiplying (1.2 ) by the integrating factor gives
e
−Lt
u
t
− e
−Lt
Lu
|
{z }
v
t
= e
−Lt
N(u),(1.7)
that is,
v
t
= e
−Lt
N(e
Lt
v).(1.8)
This has the effect of ameliorating the stiff linear par t of the PDE, and we can use
a time–stepping method of our choice (for example a fourth order Runge-K utta for-
mula) to advance the transformed equation. In practice, one doesn’t use the equation
as it is written in (1.8), but rather replaces actual time, t, with the timestep, ∆t,
and incrementally updates the formula from one timestep to the next. This greatly
improves the stability.
In both the split step method and the integrating factor method, we use a fourth
order Runge–Kutta method for the time–stepping. The fourth order Runge–Kutta
algorithm that we used to perform the time integration for this method was:
a = hf(v
n
, t
n
),(1.9)
b = hf(v
n
+ a/2, t
n
+ h/2),
c = hf(v
n
+ b/2, t
n
+ h/2),
d = hf(v
n
+ c, t
n
+ h),
v
n+1
= v
n
+
1
6
(a + 2b + 2c + d), (Fourth order RK).
where h is the time step and f is the nonlinear functional on the RHS of (1.8). Fo r
the split step method, we simply replace f in (1.9) with F from (1.4).
SL = Sliders. In a recent paper [20], Fornberg and Driscoll describ e a clever
extension of the implicit–explicit conce pt described above. In addition to splitting
the problem into a linear and a nonlinear part, they als o split the linear part (after
transformation to Fourie r space) into three regions: low, medium and high wavenum-
bers. The slider method involves using a different numerical scheme in each region.
The advantage of this method is tha t one can combine high order methods for the
low wave numbers with high-stability methods for the higher wave numbers. We can
summarise one version of this method with the following table:
4 A.-K. KASSAM AND L. N. TREFETHEN
Low |k| Medium |k| High |k|
AB4/AB4 AB4/ AM6 AB4/AM2*
Here k is the wavenumber, AB4 denotes the fourth order Adams-Bashforth formula,
AM6 denotes the sixth order Adams-Moulton for mula, and AM2* denotes a modified
second o rder Adams-Moulton formula specified by
u
n+1
= u
n
+
h
2
(
3
2
Lu
n+1
+
1
2
Lu
n−1
),(1.10)
where h is the timestep.
Unfortunately, this scheme is stable only for purely dispersive equations. In order
to generalise the concept, Driscoll has developed a very similar idea using Runge-
Kutta time–stepping [17]. Again, the idea is to make use of different schemes for
‘fast’ a nd ‘slow’ modes. In this c ase, he uses the fourth-order Runge–Kutta formula
to deal with the slow, nonlinear modes, and an implicit–explicit third order Runge–
Kutta method to advance the ‘fast’ linear modes. This is the method that we explore
in this paper.
ETD = Exponential Time Differencing. This method is the main focus of this
paper, and we will describe it in Section 2.
***
One might imag ine that extensive comparisons would have been carried out of the
behavior of these methods for various PDE s such as those listed in the first paragraph,
but this is not so. One reason is that SL and ETD are quite new; but even the other
three methods have not been co mpared as systematically as one might expect.
Our aim in beginning this project was to make such a co mparison. However,
we soon found tha t further development of the ETD schemes was first needed. As
partly rec ognized by their originators Cox and Matthews, these methods as originally
proposed encounter certain problems associated with eig e nvalues equal to or close
to zero, especially when the matrix L is not diagonal. If these problems are not
addressed, ETD schemes prove unsuccessful for some PDE applications.
In Section 2 we propose a modification of the ETD schemes that solves these
numerical problems. The key idea is to make use of complex analysis and evaluate
certain coefficient matrices or scalars via contour integrals in the complex plane.
Other modifications would very po ssibly also achieve the same purpose, but so fa r as
we know, this is the first fully practical ETD method for general use.
In Section 3 we summarize the results of experimental comparison of the five
fourth-order methods listed above for four PDEs: the Burgers, Korteweg-de Vries,
Allen-Cahn, and Kuramoto-Sivashinsky equations. We find that the ETD scheme
outperforms the o thers. We believe it is the best method currently in existence for stiff
PDEs, at least in one space dimension. In making such a bold statement, however, we
should add the caveat that we are only considering fixed time steps. Our ETD methods
do not e xtend cheaply to variable time-stepping; an IMEX scheme, for example, is a
more natural candidate for s uch problems.
Sections 4 and 5 illustr ate the methods in a little more detail for a diagonal
example (Kuramoto-Sivashinsky) and a non- diagonal e xample (Allen-Cahn). They
also provide brief Matlab codes for use by our re aders as templates.
2. A modified ETD scheme. Low order ETD s chemes arose originally in the
field of computational electrodynamics [59]. They have been independently derived
FOURTH-ORDER TIME-STEPPING FOR STIFF PDE 5
several times [5, 12, 14, 21, 46, 48] – indeed Arieh Ise rles has pointed o ut to us that in
the ODE context, related ideas g o back as far as Filon in 1928 [18, 30] – but the most
comprehensive treatment, and in particular the fourth order ETDRK4 formula, is in
the paper by Cox and Matthews [14], and it is from this paper that we take details of
the scheme. Cox and Matthews argue that ETD schemes outper form IMEX schemes
because they tre at transient solutions better (where the linear term dominates), and
outperform IF schemes because they treat non-transient solutions better (where the
nonlinear term dominates).
Algebraically, ETD is simila r to the IF method. The difference is that we do not
make a complete change of variable. If we proceed as in the IF approach and apply
the same integrating factor and then integrate over a single time step of length h, we
get
u
n+1
= e
Lh
u
n
+ e
Lh
Z
h
0
e
−Lτ
N(u(t
n
+ τ), t
n
+ τ)dτ.(2.1)
This equation is e xact, s o far, and the various order ETD schemes come fr om how
one approximates the integral. In their paper Cox and Matthews first present a
sequence of recurrence formulae tha t provide higher and higher order approximations
of a multistep type. They pro pose a generating formula
u
n+1
= e
Lh
u
n
+ h
s−1
X
m=0
g
m
m
X
k=0
(−1)
k
m
k
N
n−k
,(2.2)
where s is the order of the scheme. The coefficients g
m
are given by the recurre nce
relation
Lhg
0
= e
Lh
− I,
Lhg
m+1
+ I = g
m
+
1
2
g
m−1
+
1
3
g
m−2
+ . . . +
g
0
m + 1
, m ≥ 0.(2.3)
Cox and Matthews also derive a set of ETD methods based o n Runge-Kutta time
stepping, which they call ETDRK schemes. In this report we consider only the fourth
order scheme of this type, k nown as ETDRK4. According to Cox and Matthews, the
derivation of this scheme is not at all obvious and required a symbolic manipulation
system. Here are the formulae:
a
n
= e
Lh/2
u
n
+ L
−1
(e
Lh/2
− I)N(u
n
, t
n
),
b
n
= e
Lh/2
u
n
+ L
−1
(e
Lh/2
− I)N(a
n
, t
n
+ h/2),
c
n
= e
Lh/2
a
n
+ L
−1
(e
Lh/2
− I)(2N(b
n
, t
n
+ h/2) − N(u
n
, t
n
)),
u
n+1
= e
Lh
u
n
+ h
−2
L
−3
{[−4 − Lh + e
Lh
(4 − 3Lh + (Lh)
2
)]N(u
n
, t
n
)
+ 2[2 + Lh + e
Lh
(−2 + Lh)](N(a
n
, t
n
+ h/2) + N(b
n
, t
n
+ h/2))
+ [−4 − 3Lh − (Lh)
2
+ e
Lh
(4 − Lh)]N (c
n
, t
n
+ h)}.
Unfortunately, in this form, ETDRK4 (and indeed any of the ETD schemes of
order higher than two) suffers from numerical instability. To understand why this is
