Partial: Implicit-Explicit Methods For Time-Dependent Differential Equations Ascher, Wetton

SIAM J. NUMER. ANAL.
1995 Society for Industrial and Applied Mathematics

Vol. 32, No. 3, pp. 797--823, June 1995 O06
IMPLICIT-EXPLICIT METHODS FOR TIME-DEPENDENT

Downloaded 12/29/12 to 128.148.252.35. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
PARTIAL DIFFERENTIAL EQUATIONS*

URI M. ASCHER, STEVEN J. RUUTH$, AND BRIAN T.R. WETTON
Abstract. Implicit-explicit (IMEX) schemes have been widely used, especially in conjunction
with spectral methods, for the time integration of spatially discretized partial differential equations
(PDEs) of diffusion-convection type. Typically, an implicit scheme is used for the diffusion term
and an explicit scheme is used for the convection term. Reaction-diffusion problems can also be
approximated in this manner. In this work we systematically analyze the performance of such
schemes, propose improved new schemes, and pay particular attention to their relative performance
in the context of fast multigrid algorithms and of aliasing reduction for spectral methods.
For the prototype linear advection-diffusion equation, a stability analysis for first-, second-,
third-, and fourth-order multistep IMEX schemes is performed. Stable schemes permitting large time
steps for a wide variety of problems and yielding appropriate decay of high frequency error modes are
identified. Numerical experiments demonstrate that weak decay of high frequency modes can lead
to extra iterations on the finest grid when using multigrid computations with finite difference spatial
discretization, and to aliasing when using spectral collocation for spatial discretization. When this
behavior occurs, use of weakly damping schemes such as the popular combination of Crank-Nicolson
with second-order Adams-Bashforth is discouraged and better alternatives are proposed.
Our findings are demonstrated on several examples.
Key words, method of lines, finite differences, spectral methods, aliasing, multigrid, stability
region
AMS subject classifications. 65J15, 65M20
1. Introduction. Various methods have been proposed to integrate dynamical

systems arising from spatially discretized time-dependent partial differential equations
(PDEs). For problems with terms of different types, implicit-explicit (IMEX) schemes
have been often used, especially in conjunction with spectral methods [7], [16]. For
convection-diffusion problems, for example, one would use an explicit scheme for the
convection term and an implicit scheme for the diffusion term. Reaction-diffusion
problems can also be approximated in this manner. In this work we systematically
analyze the performance of such schemes, propose improved new schemes, and pay
particular attention to their relative performance in the context of fast multigrid
algorithms and of aliasing reduction for spectral methods.
Consider a time-dependent PDE in which the spatial derivatives have been dis-
cretized by central finite differences or by some spectral method. This gives rise to
a large system of ordinary differential equations (ODEs) in time which typically has
the form
(1) it- f(u) / ,g(u)
,
where Ilgll is normalized and is a nonnegative parameter. The term f(u) in (1) is
some possibly nonlinear term which we do not want to integrate implicitly. This could
Received by the editors June 1, 1993; accepted for publication December 20, 1993.
Department of Computer Science, University of British Columbia, Vancouver, British Columbia,
V6T 1Z4, Canada (ascher(C)cs.ubc. ca). The research of this author was partially supported by the
Natural Sciences and Engineering Research Council of British Columbia grant OGP0004306.
Department of Mathematics, University of British Columbia, Vancouver, British Columbia,
V6T 1Z2, Canada (ruuth(C)math.ubc. ca). The research of this author was partially supported by the
Natural Sciences and Engineering Research Council of British Columbia postgraduate scholarship.
Department of Mathematics, University of British Columbia, Vancouver, British Columbia,
V6T 1Z2, Canada (wetton(C)math. ubc. ca). The research of this author was partially supported under
Natural Sciences and Engineering Research Council of British Columbia grant OGP0122105.
797
798 U. M. ASCHER, S. J. RUUTH, AND B. T. R. WETTON
be because the Jacobian of f(u) is nonsymmetric and nondefinite and an iterative

solution of the implicit equations is desired, or the Jacobian could be dense, as in
spectral methods, requiring the inversion of a full matrix at each time step. One may
simply wish to integrate f(u) explicitly for ease of implementation. The term ug(u),
however, is a stiff term which should be integrated implicitly to avoid excessively
small time steps. Frequently ug(u) is a linear diffusion term, in which case the implicit
equations form a linear system which is positive definite, symmetric, and sparse. Such
_
systems can be solved efficiently by iterative techniques (e.g., [26], [11]). Thus, for
problems of the form (1) it often makes sense to integrate ug(u) implicitly and f(u)
explicitly, yielding an IMEX scheme.
The most popular IMEX scheme hitherto has been a combination of second-order
Adams-Bashforth for the explicit ("convection") term and Crank-Nicolson for the
implicit ("diffusion") term [16], [7], [1]. Applied to (1) this gives
un+ n 3 1
1)
(2) f(un) -f(u [g(tn+l) nt g(tn)]
where k is the constant discretization step size and u is the numerical approximation
to u(kn). Other IMEX methods have also been considered, e.g., [18], [25], [15].
A wide variety of other applications for IMEX schemes are also possible. For
example, solutions to reaction-diffusion systems arising in chemistry and mathemati-
cal biology can be computed using this technique. For these problems the nonlinear
reaction term can be treated explicitly while the diffusion term is treated implicitly.
Examples of reaction-diffusion systems from a biological standpoint can be found in
[19]. We report about these in a separate paper [21].
Several authors have analyzed specific IMEX schemes. For example, an experi-
mental analysis of several IMEX schemes including (2) was carried out in [2]. Some
stability properties for certain second-order IMEX schemes were determined in [25].
The schemes (17) and (24) below were considered in a Navier--Stokes context in [15].
These studies do not address how to choose the best IMEX scheme for a given sys-
tem (1).
In this paper, we seek efficient IMEX schemes for convection-diffusion-type prob-
lems using a systematic approach. First we derive, in 2, classes of linear multistep
IMEX schemes. Classes of s-step schemes turn out to have optimal order s and to
depend on s parameters.
We then restrict our attention in 3 to a prototype advection-diffusion equation
in one dimension
(3) Ut aUx + ’Uxx

where a and u are constants, > 0, subject to periodic boundary conditions. A
von Neumann analysis (see, e.g., [23]) can then be applied, yielding in effect a safe
diagonalization of (1) and allowing us to consider a scalar test equation
(4) 2 (c + i/)x
as is customary for ODEs (see, e.g., [12]) for values of a,/3 on an ellipse. An analysis
of IMEX schemes for (4) is then performed, seeking methods which allow the largest
stable time steps. In particular, for u >> 1 we seek methods possessing a mild time
step restriction since the system (1) is dominated by the implicitly handled "diffusion
term." For other values of u >_ 0 we seek schemes which have reasonable time-stepping
restrictions, e.g., comparable to the Courant-Friedrichs-Lewy (CFL) condition [23].
IMPLICIT-EXPLICIT METHODS FOR TIME-DEPENDENT PDEs 799
In 4 we perform a variety of numerical experiments using finite difference and

pseudospectral methods on linear and nonlinear problems of convection-diffusion type
in one and two spatial dimensions. These experiments agree with the theory of 3.
Furthermore, we discuss and demonstrate the use of IMEX schemes which yield strong
decay of high frequency spatial modes. This property has important implications for
the efficiency of time-dependent multigrid methods and of pseudospectral methods.
An appropriate IMEX scheme (but not the popular (2)!) can reduce aliasing in pseu-
dospectral methods, in the multigrid context, similar methods can reduce the number
of multigrid cycles needed per time step, in effect acting as smoothers (cf. [14], [10]).
Conclusions and recommendations are summarized in 5. Perhaps the most sur-
prising conclusion is that the most popular IMEX scheme (2) can essentially always
be outperformed by other IMEX schemes. Moreover, a modification of (2) is almost
always at least as good as (2) and at times much better.
2. General linear multistep IMEX schemes. We now derive s-step IMEX
--
schemes for (1), s > 1. Letting k be the discretization step size and u denote the
approximate solution at tn kn, we may write these schemes as
s-1 s-1 s-1
1
tn-t-l" E
j=o
ajtn-j E bJf(un-J)-t-P E cJg(ztn-J)
j=0 j=-
where c-1 0. See [8] for some stability and convergence results. For a smooth
function u(t), expand (5) in a Taylor series about tn nat to obtain the truncation
error. This yields
Applying (1) to the truncation error (6), an order p scheme is obtained provided
s-1
1 + Eaj 0,
j=0
s-1 s-1 s-1
(7)
1- E
j=
j aj E E j=0
by
j=-
cj,
1
j=l
ay E jbj c_
j=l j=l
jcj,
(p 1)! (p- 1)! + (p- 1)!

j’-I j--1 j--1
It is not difficult to prove the following theorem [22].

THEOREM 2.1. For the s-step IMEX scheme (5), the following hold.
(a) The 2p + 1 constraints of the system (7) are linearly independent, provided
p <_ s. Thus, there exist s-step IMEX schemes of order s.
(b) An s-step IMEX scheme cannot have order of accuracy greater than s.
(c) The family of s-step IMEX schemes of order s has s parameters.
We thus restrict all further discussion to s-step IMEX schemes of order s.
3. Analysis for a test advection-diffusion problem. For the prototype prob-
lem (3), using the standard second-order centred approximations D1 and D2 for the
first and second derivatives, respectively, we obtain the corresponding semidiscrete
equations
aDIU + D2U, 1 <_ <_ M.

(Here, a uniform spatial grid with M points has been employed.) Applying a discrete
Fourier transform diagonalizes this system to
2i i3jxi + aixi, j 1,...,

where ai and/i are given by
(8) (oj,/3j) ( 2[cs(2rjh)-l]’a sin(2rjh)
For notational convenience, we write

(9) 2 ix + ax.
We then consider values (a, 3) which lie on the ellipse of Fig. 1. Note that the solution
of
(10)
(9) decays in time,
_ j--0
]X(tn+)]
Applying the general multistep IMEX scheme (5) to
lxn+l ajx
n-j
E bjixn-J
j=0
-nt-
(9) yields
j=--i
which is a linear difference equation with constant coefficients. (When comparing (10)
to (5) note that u is buried in a.) The solutions of the difference equation (10) are of
CjCtX
n-j
the form
X n+l plT + p2T +’’" + psT"ns
where -j is the jth root of the characteristic equation defined by
s-1
(11) (I:)(Z) (1 C_lCtlg)Z -t- E(aj bjik cjak)z s-j-1
j=O
_1-1-1111"] VI’ []--[ IV1 I[TT-TI Ill I[11 ]ii ilT’[ ilil iTT-[TII 1[ ["1 ]1I I[] I1 I1 ill1
a/h
-4v/h
Of
FIG. 1. Half ellipse of possible (c, 1).
and pj is constant for -j simple and a polynomial in n otherwise. Clearly, stability

holds for ITjI<_ 1, -j simple, and ITjI< 1 where Tj is a multiple root.
Remark. Below we will be interested in the entire range of u _> 0. For the purely
hyperbolic case u 0, our choice of the centered difference D1 then leaves out the
important characteristic-based methods. However, it is well known that a one-sided
scheme can be interpreted as a centred scheme plus an artificial diffusion term (e.g.,
[13]), so we do want to pay close attention to
3.1. First-order IMEX schemes. The one-parameter family of first-order IMEX
schemes for (1) can be written as
(12) U n+l
and we restrict 0 <_ /_< 1 to prevent large truncation error.

The choice 0 yields the forward Euler scheme
t n+l
n
t k f(u ) + l/kg(tn).
This scheme is fully explicit, and will not be considered further. The choice
yields the second-order Crank-Nicolson scheme when f 0. Another possibility is to
apply backward Euler to g and forward Euler to f. This choice (/= 1) yields
(3) un+ un kf(un) + l]]gg(u n+l).

IMEX schemes such as (13) which apply a backward differentiation formula (BDF)
discretization [9] to g and which extrapolate f to time step (n + 1) will be referred to
as semiimplicit BDF (SBDF) schemes.
Applied to the test equation (4), the general IMEX scheme (5) gives
X n+l (0, )X n, where (c, )

1 + ka(1 -/) +
The stability region is thus {(c,/): I(c,)l _< 1}.

Note that these one-step schemes trivially accommodate variable time-stepping
and use relatively little storage. The choice 1 is particularly attractive because
strong decay occurs for c large and negative. However, first-order IMEX schemes have
the shortcoming that they are unstable near the nonzero -axis (i.e., for sufficiently
small and a 0 fixed but arbitrary) since
I( 0,)1 11+ ik] v/l+ k2 2 > 1.

Furthermore, at least
second-order spatial discretization is used. We thus consider higher-order methods
for the remainder of this paper.
3.2. Second-order IMEX schemes. Approximating (1) to second order using
IMEX schemes leaves two free parameters. If we center our schemes in time about
time step (n +-) to second order, we obtain the following family:
(14)
1
[( l)un+ ( 1)un--1] 1)
+p [( + )g(n+l)+ (1--c)g(n) + g( --1)].
C C n
Some of these methods are quite familiar. For example, selecting (, c) (, 0) yields
the popular scheme (2),
1
(15) 1 [un+ 3
.
Because it applies Crank-Nicolson for the implicit part and second-order Adams-
Bashforth for the explicit part, this scheme will be referred to
1
This choice gives
CNAB (Crank-
Nicolson, Adams-Bashforth). Below, we show that the best asymptotic decay prop-
erties for 7 occur when c
[9(un+l) + 9(un)+ g(n-l

16
1
Because of the obvious similarity to CNAB, this scheme will be called modified CNAB
(MCNAB). Note that in comparison to CNAB, MCNAB does require the additional
evaluation or storage of g(u n-l).
By setting (7, c) (0, 1) in (14) we obtain another method which has been applied
to spectral applications (e.g., [4])
1 n-1
(16) 2 U f(un) + [g(un+l) + g(un--1)]"
This scheme uses leap frog explicitly and something similar to Crank-Nicolson im-
plicitly (cf. [17]). For this reason, this method shall be referred to as CNLF (Crank-
Nicolson, leap frog).
Finally, setting (7, c) (1, 0) yields

1
(17) 2-- [3un
+1 -4u n + tn-l] 2f(tn) f(tn-l) d- pg(tn+l)
which shall be referred to as SBDF since this scheme is centered about time step
(n + 1). Other authors, e.g. [25], call this scheme extrapolated Gear.
Having determined integration formulae, we direct our attention to obtaining
stability properties for second-order IMEX schemes. The second-degree characteristic
.
polynomial resulting from (14) applied to 2 (c + i/)x is given by
(18)
1
+-k (7+c) z
[Z + iZ( + ) + (1
1 )]z +
1
+ iZ
c
It is easy to verify that at the origin (a,/) (0 0) the roots of (I) are 1 and 27-1
2"y+l
Thus, all of these schemes are stable at the origin, provided u >_ 0.
Because the parameter c in (18) is always multiplied by a we choose c according
to stability properties for lal >> 1/4. For this case, the roots of the characteristic
equation (18) are given approximately by
(7+c ) +(1-7- c)r +

c
0
which yields
7 + c- 1 + V/(1 7) 2c
7-1, 2
27+c
For any (7, c), evaluating
(19)
provides an estimate of the high frequency modal decay for large u. Minimization
over c determines the method with the strongest asymptotic high frequency decay for
a particular 7. This yields (see [22])
(20) c= for 7 > O,

2
1
c>_ for 7=0.
The choice
(21) c 1 7
also proves useful (see [22]). The schemes SBDF, MCNAB, and CNLF satisfy (20),
the schemes SBDF and CNLF satisfy (21), but CNAB satisfies neither. Also, mini-
mization of (19) over 7 and c determines that the SBDF scheme (17) possesses the
strongest asymptotic decay of second-order methods.
Further information can be obtained from the stability contours in the a-
.
plane. These plots are displayed in Figs. 2-5. Figure 2 shows the contours for CNLF.
This method is stable for all > 0, provided k < h Such a time step restriction is
804 U.M. ASCHER, S. J. RUUTH, AND B. T. R. WETTON
I
r]-T]iiTT-rT
O/k -lO/k 0
( 0(,
FIG. 2. Stability contours for CNLF. FIG. 3. Stability contours for CNAB.
-lO/k 0 -lO/k 0
i 0(,
FIG. 4. Stability contours for MCNAB. FIG. 5. Stability contours for SBDF.
high frequency modes can be weak, tending to 1 as a

to other second-order methods such that y
largest stability region among such methods.
-
undesirable since it applies even for large u and small h. Furthermore, the decay of
-. Comparison of CNLF
0 suggests that CNLF produces the
The contours for CNAB are plotted in Fig. 3. This method has a reasonable time
step restriction for larger and small h. It is unstable near the/-axis, however. It
a
-
also suffers from poor decay of high frequency modes, since the decay tends to 1 as
-c. Using MCNAB (7, c) ,
(5, ) the decay tends to which is a significant
improvement. See Fig. 4 for these contours.
The contours for SBDF are displayed in Fig. 5. This method has the mildest
time-stepping restriction when is large and h is small. The decay of high frequency
modes is strong, tending to 0 as a tends to -. This method, however, has the
strictest time step limitation for small lal.
We can now develop a quantitative method for describing time step restrictions
for second-order schemes. Such a method will help select which second-order scheme
to use for a particular problem.
Defining the mesh Reynolds number [20]
ah
=_
10-I
10-2 10-I 10 101
Vi y/h =I/R
FIG. 6. Time step restriction for various where c /.
in Figs. 6-8 we plot the theoretical time step restrictions for various and c. As
can be seen from these figures, increasing ? allows larger stable time steps when
R < 2. The case c 1 --y also has the property that decreasing "y allows larger
time steps for R > 2. Comparison of Figs. 6-8 indicates that the largest time step
can be applied using SBDF for R < 2 and CNLF for R > 2. This result physically
corresponds to selecting SBDF when discrete diffusion dominates, and CNLF when
discrete convection dominates. From this perspective, the popular CNAB is only
competitive when R 2.
Frequently, an important consideration when choosing a second-order scheme is
what the constant of the truncation error is. For example, Crank-Nicolson is known
to have a much smaller truncation error than second-order BDF (see, e.g., [12]),
so we expect CNAB to have a smaller truncation error than SBDF. The scheme
MCNAB is expected to have a truncation error similar to CNAB, however. (Numerical
experiments in 4 support this claim.) Because of its small truncation error and
because it produces stronger high frequency spatial decay than CNAB, MCNAB may
be preferred in certain problems over both CNAB and SBDF when R < 2.
One disadvantage that all second-order IMEX schemes with "7 > 0 (i.e., essentially
all except CNLF) share is that their stability regions do not contain a portion of the
-axis other than the origin. Specifically, it was shown in [25] that when ( 0 one
of the roots -i of (18) satisfies
-
I-i(flk)l 2-1+ /2+ (ilk)4+...>
for flk sufficiently small and > 0. Thus for all k, all second-order schemes such
that -y > 0 are unstable on the nonzero fl-axis. The CNLF scheme is well known to
,
be stable on the fl-axis provided k < h but it is only marginally stable (see, e.g.,
1
" 10 0
ii II
10-2 10 -I 10 0 101 1
Vi ty/l =I/R
FIG 7. Time step restriction for various y where c

l_.y)2
2
101
FIG. 8. Time step restriction for various "y where c O.

IMPLICIT-EXPLICIT METHODS FOR TIME-DEPENDENT PDEs
[23]), providing no damping of high frequency error components anywhere. To obtain

IMEX schemes which are stable for u _> 0 and have strong decay for la[ >>
- 807
we must
consider higher-order schemes.

3.3. Third-order IMEX schemes. Recall from 2 that we have a three-param-
eter family of three-step IMEX schemes of order three. One possible parameterization
yields
(23)
These schemes are centered about time step (n + 7) to third order, provided 0 0.
As for lower-order schemes, the value of 7 should be between 0 and 1 to avoid large
truncation error. Also, the parameter c is multiplied by u, so this parameter should
be adjusted to modify large properties of a scheme.
-
Setting (7, 0, c) (1, 0, 0) yields
(24)
1
un+l 3U n +
1 n--2 3f(un) 3f (tn-- )+ f(u 2) + pg (tn+
k -U )
which is the third-order BDF for the implicit part, and therefore is called third-order
SBDF.
--
The third-degree characteristic polynomial resulting from (23) applied to the test
equation ( + i)x is given by
(25)
(z)= [ ++5+o- ( 72 + c) akl z3
1 7
72 + 37
3
-.y 2 1
+,--+o+ +1+ /
ik + (1- 72 + 2 0) ck] z2
3c
+ +-1+ +2+0 i+ 2
3c+
04) 1ok z
-g+
We now determine which methods produce the strongest asymptotic decay as a
-c. For this case, the roots of the characteristic polynomial (25) are given approxi-
mately by
(26)
72
-t-c
+y
2 ) za+ ( 1-72-3c+]-230h (7-72 ]
z2_
2 3c+-0 z + ---O-c-O.
By determining the solutions, 7"1, T2, and 73, of (26) we may evaluate
to obtain an estimate of the high frequency model decay for large Minimization
over (7, 0, c) determines the method with the strongest asymptotic high frequency
.
decay. Certainly if
(27) D3"o,Oo,co 0
then (7o,0o, co) minimizes the amplification as a -oc. From (26) we satisfy (27)
iff
23
1 -7g 3co + 0o
=0,
3") +3"o -+Co
2
3"g-3"o
(28) 2 +3co-400 =0,
3"( + 3"0 -4- CO
2
--0o
12 co =0.
3"( + 3"0
+co
2
Each term is divided by 2

(g+3"o
+ co) to allow the possibility of satisfying (27) by
+3"0
letting 2 + co) -* +oc. Simplification of the system (28) yields
(370- 1)(7o- 1)
(29) 3" +23’0
0,
(30) 0o 6(7o2 7o) 0,
-
3" +3"o
2
(31)
--0o
12 co 0.
3"g +3"0 CO
2
For (+o
2 + co) to be finite there are two possibilities, (70, 0o, co)(1, 0, 0) and
4
,
(/o,0o,co)=(3, 3, 5) both of which specify third-order SBDF For (3" +3"0
2 +co)
to be infinite, (29) implies Iol << Icol. Using this in (30) implies 10ol << Icol. Applying
these results to (31) results in a contradiction, so 2 (3"o:+3"o
+ co) must be finite.
Because third-order SBDF has the strongest asymptotic decay of third-order
IMEX schemes, special attention is given to its properties throughout the remain-
der of this section. All schemes considered here are stable on a segment of the -axis
including the origin (a, ) (0, 0), as can be verified by an analysis of (25) [22].
Further information about stability in the a-
plane can be obtained by plotting
max{Izl (z) 0}, where q(z) is defined in (25). These stability contours are
displayed in Figs. 9-14.
We begin by examining if it is possible to arrive at a stable scheme for any fixed

7. For a fixed, but arbitrary, 7 and for 101, Icl --, oc we obtain an approximate local
minimum of T,0,c if 0.4661. Using these parameters, the scheme simplifies to
(32)
1
(u+l u
23 4
-f(un) -f(un-l) +
+ .5184ug(u ) + .0650ug(u
2 -1)n-2) +
f(u .4661-g(un+l)
.0494ug(u -2)
which applies third-order Adams-Bashforth to the explicit term. (Note that 7 disap-
pears in (32) through the limiting process.) The stability contours of Fig. 9 suggest
that this method is stable for all u > 0 provided k < 0.62 h. This restriction is more
severe than that for third-order Adams-Bashforth applied to 2 ix because of the
dip in the stability contours when c < 0. As mentioned previously, the -axis is in-
cluded in the absolute stability region. This result demonstrates that for third-order
methods, it is possible to find methods for any 7 which are stable for all u > 0 by
varying 0 and c.
In 3.2, the most interesting second-order schemes were produced by selecting
equal to 0, 3, or 1. We consider schemes for these values of 7 below. The parameters
0 and c are chosen to produce schemes which allow large stable time steps as
For example, the method (7, 0, c) (0,-2.036,-.876) of Fig. 10 is stable for all
u >_ 0 provided k <_ 0.67. h
Similarly (7,0, c) (.5 -1.21 5) of Fig. 12 is stable
for all u >_ 0 if k <_ 0.65-ha
In both these cases substantially larger time steps can be
taken for large or moderate lal than for the method (32). Furthermore, stronger high
frequency decay occurs for these methods than for method (32).
Recall that the strongest decay as a --, -oc occurs for third-order SBDF. The
stability contours for this method are shown in Fig. 13. This plot together with the
zoom-in of Fig. 14 suggest that third-order SBDF is stable for all a < 0 provided
k < 0.62 h The plot of Fig. 13 also indicates that very large time steps can be taken
I.
for large or moderate [a Although applying 7 1 and nonzero 0 and c may allow
somewhat larger stable time steps we focus on 0 c 0 since other choices degrade
the strong asymptotic decay.
For 7 1, the choice 0 c 0 is not recommended. As seen in Fig. 11,
small time steps are needed for moderate or large

lie in the stable region.
,
(7, 0, c) (3,0, 0) results in a small stability region. This plot indicates that very
since the ellipse of Fig. 1 must
3.4. Fourth-order IMEX schemes and a comparison. From 2, the general
-- -
fourth-order, four-step IMEX scheme is a four parameter family of methods. From the
previous section we anticipate that the fourth-order SBDF may have good stability
properties. For/- f(u)+ vg(u), this scheme is given by
+ 3 u -1 4un-2
1
(33) 25un+l 4u
4f(un) 6f(u n-l) + 4f(un-2) f(un-3) + Pg(un+l).

The characteristic polynomial obtained from applying (33) to 2 (a + i)x is
(z) (25) (4 + + (3 + (4)

5+ z +
1
Contour plots similar to those in Figs. 13 and 14 were obtained for the fourth-order
SBDF as well (see [22] for these plots). The scheme is stable for k <_ 0.52 h and
-9/k 0
O
FIG. 9. Stability contours for FIG. 10. Stability contours for
(% 0, c) lim (7, 0, .46610).
0-- (’7, O, c) (0,-2.036,-.876).
-9/k 0 -9/k 0
0(, 0,
(,0,)
FIG. 11. Stability contours
(1/2,0,0),
for FtG. 12.
(., o, ) (-, .,
Stability contours
-.).
for
1/k 1/k
-.2/k 0
FIG. 13. Stability contours for 3 rd FIG. 14. Zoom-in around 0 for
order SBDF, (7, O, c) (1, 0, 0). 3 rd Order SBDF.
larger time steps are permitted as u increases. However, the stability region is smaller
than for the third-order case, so smaller time steps may be required. Furthermore,
the -axis is closer to the boundary of the stability region for fourth-order SBDF,
suggesting that third-order SBDF may dissipate error better for u 0.
Third- and fourth-order SBDF methods are good choices for IMEX schemes for
some problems. For all >_ 0, these methods are stable for a time step restriction
101
100 CNLF // O( CNLF CNLF
&er4
O2orr SBDF
Or
10-1 II I,II II
10-2 10-1 100 101 1
v, / =I/R
FIc. 15. Time step restrictions for various IMEX schemes.
similar to the CFL condition. Greater accuracy and strong high frequency decay also
make these methods very attractive. Nonetheless, for many problems second-order
methods are preferred. Higher-order methods require more storage, and more work
per time step. This extra expense could be justified if greater accuracy permitted
-
larger time steps. Third- and fourth-order schemes, however, have more severe time
step restrictions than second-order schemes. In fact, Fig. 15 shows that larger stable
time steps can be taken with second-order SBDF when R < 9 for the linear advec-
tion-diffusio problem using second-order spatial discretizations. CNLF allows larger
stable time steps than either third- or fourth-order SBDF for R > 1.
Third-order SBDF should be efficient for problems which require strong decay for
Ic >> and a moderate time step restriction for R > 10. It should also be effective
for problems where R >> 1, since a portion of the -axis is within the stability region.
Fourth-order SBDF has a particularly severe time step restriction for the advec-
tion-diffusio problem when R is moderate or small. For example, when R .125,
fourth-order SBDF can only apply one tenth the time step of second-order SBDF,
as can be seen from Fig. 15. This restriction on the step size would appear to limit
fourth-order SBDF to problems where accuracy, and not stability, is the reason for
limiting the time step size.
4. Further considerations and numerical experiments. The previous sec-
tion has dealt with stability properties of IMEX schemes for the one-dimensional
linear constant coefficient advection-diffusion equation. These results provide neces-
sary, but not sufficient conditions for stability for variable coefficient and nonlinear
convection-diffusion problems.
In this section we summarize numerical experiments which verify that our analysis
for the simple, linear problem can be useful for determining which IMEX scheme to
apply to more complicated problems. In addition, strongly damping schemes are

shown to be more efficient in certain spectral collocation and multigrid applications.
In order to calculate starting values for multistep IMEX schemes, we use one-step
(low-order) IMEX schemes with a very small time step, unless otherwise noted.
4.1. Finite difference approximations in one dimension.
Example 1. To examine nonzero viscosity behaviour, we consider the one-dimen-
sional variable coefficient problem
(34) ut + sin(2rx)u.
subject to periodic boundary conditions on the interval [0,1] and initial condition
u(x, O) sin(27rx).
We use centered second-order differences for the spatial derivatives. Note that the
solution is smooth for all u > 0.
To test the theory’s predictions for small mesh Reynolds numbers (22), the model
problem (34) was approximately solved using discretization step sizes h 3
and k 1.8h in time. Use of such step sizes is appropriate only for strongly damped
in space
flow. Utilizing several IMEX schemes, computations to time t 2 are performed for
viscosities u in the range .01 _< u _< .1. These values correspond to mesh Reynolds
numbers R in the range 1.59 >_ R > 0.159.
From Figs. 6, 7, 8, and 15 the theory predicts that for these step sizes SBDF is
stable for a larger viscosity interval than any other scheme. As u decreases, third-order
SBDF followed by CNAB should become unstable. Stability for MCNAB should be
similar to, and somewhat better than, CNAB. Furthermore, fourth-order SBDF and
CNLF should be unstable over the entire interval because the theoretical stability
restriction is violated.
By comparing results with those for h and k .225h using SBDF, the
max norm relative errors for second-order schemes are evaluated. Third- and fourth-
order schemes are compared with computations with these same step sizes but using
third- and fourth-order SBDF. The resultant errors are plotted against u in Fig. 16.
Fourth-order SBDF and CNLF are not included because they are unstable over the
entire interval. The plots of the figure clearly coincide with the results of the theory.
Figure 16 also indicates that SBDF is the only stable method for the above choice
of discretization parameters when 0.0015 < u < 0.0025. (The values of u where the
curves turn upward correspond to the onset of instability for the given values of h and
k.) This agrees with the prediction that SBDF allows the largest stable time steps for
small mesh Reynolds numbers (see 3.2). Although third-order SBDF has a smaller
stability interval, it may be useful in problems where high accuracy is needed since it
produces a smaller error than second-order methods when stable.
To test the theory’s predictions for large mesh Reynolds numbers, example (34)
was solved using discretization step sizes h- in space and k .9h in time. Using
several IMEX schemes, computations to time t- 2 are performed for viscosities u in
the range .001 <_ u <_ .01. These values correspond to mesh Reynolds numbers R in
the range 12.3 > R >_ 1.23.
From Figs. 6, 7, 8, and 15 the theory predicts that for these step sizes only CNLF
is stable over the entire viscosity interval. As u decreases, third-order SBDF followed
by SBDF and finally CNAB should become unstable. Stability for MCNAB should
be similar to CNAB. Furthermore, fourth-order SBDF should be unstable over the
101
10
10-1
10-2
10-3
10-2 10-1
Viscosity
FIG. 16. Large viscosity behaviour for various IMEX methods.
entire interval because the theoretical stability restriction is violated when u < .05
for these step sizes.
By comparing results with those for h and k .225h using SBDF the max
norm relative errors for second-order schemes are evaluated for these computations.
Third- and fourth-order schemes are compared with computations with these same
step sizes but using third- and fourth-order SBDF. The resultant errors are plotted
against u in Fig. 17. Fourth-order SBDF is not plotted because it is unstable over the
entire interval. The plots of the figure support the results of the theory.
Figure 17 also indicates that CNLF is the only stable method for u < 0.002. This
agrees with the prediction that CNLF allows the largest stable time steps for large
mesh Reynolds numbers (see 3.2). CNLF is a particularly attractive choice because
it has the smallest error of second-order methods. Use of SBDF is not recommended
because it has the smallest stable interval and the largest error among any of the
second-order schemes considered.
For the same problem using k .5h, dramatically different results are predicted
because all the methods satisfy their theoretical stability restrictions. Indeed, com-
putations for CNLF, third- and fourth-order SBDF all produce errors which nearly
coincide, since spatial error dominates the solutions. Other second-order methods
appear stable and produce only slightly less accurate results. []
Example 2. To examine the limit u 0, we consider the one-dimensional nonlinear
problem
1
(35) + cos(2rt)(1 + u)u 0
having periodic boundary conditions on the interval [0,1] and initial conditions
u(x, O) sin(27rx).
10
10-I
0-2
Viscosity
FIG. 1T. Small viscosity example for various IMEX methods.
As in the previous two sections, we use second-order centered differences to ap-

proximate ux. For h- 0 and k- .5h, computations are performed to time t- 100.
The time step value k .5h was used to ensure that fourth-order SBDF satisfied the
stability restriction k <_ .52h. Using the fact that the exact solution to this problem
at integer t equals the initial data,, i.e.,
u(x, n) sin(2rx), n 0, 1, 2,...,

we compute the error in the solution at time t n, n 1, 2,..., 100.
The results again agreed with the predictions from 3: All second-order schemes
tested, such that 7 > 0, were unstable, while CNLF, third-order SBDF, and fourth-
order SBDF were all stable. The SBDF schemes are dissipative as well. For fuller
details, see [22]. rn
4.2. Spectral collocation.
Example 3. For our next application we consider the well-known Burgers equation
(36) u + UUx UUxx
subject to periodic boundary conditions on the interval [-1,1] and initial conditions
u(x, O) sin(rx).
A plot of the solution of the Burgers equation for u .01 at several different times
is provided in Fig. 18. This computation was produced using Chebyshev collocation
with 40 basis functions and k 1/160 using SBDF.
The next few subsections discuss Fourier and Chebyshev collocation implementa-
tions for the above model problem. See [71 or [3] for details on these methods.
-.2
-.4
-.b
-.8
.0 .8 .6 -.4 -.2 .2 .4 .6 .8
FIG. 18. Burgers equation for various for v 0.01.
4.2.1. Fourier spectral collocation. Since the problem of Example 3 is pe-

riodic, we expect that Fourier series makes a good basis of trial functions for this
problem. Indeed, since
10u O:u
20x Ox 2
we see that ut is antisymmetric for u antisymmetric. By selecting an initial condition
which is antisymmetric we guarantee that u remains antisymmetric for all t. Since
only these components of the series contribute to the solution we use a Fourier sine
series. We thus approximate u by the series
N
E aj(t)sin(jx).
j=l
Using initial conditions al (0) 1, ak(0) 0, and k 1, the cj(t) are determined by
enforcing the differential equation at the collocation points xj 1 <j <N
i.e.,
OUN
+ u OUN " axe
This scheme is also called a pseudospectrM method since the nonlinear convection
term is evaluated in physical space.
By applying Fourier sine collocation with 40 basis functions and k 1/40, we
solve the model problem at each time step to time t 2. Because the system is small,
the implicit equations are solved in physical space using LU decomposition. In larger
10
10-I
10 -2
10-4
10-5
10-3 1-2
Viscosity
FIG. 19. Fourier spectral collocation for Burgers equation.
systems where greater efficiency is needed these would be solved in Fourier space using
transform methods (see [7]). For CNAB, MCNAB, SBDF, and third-order SBDF the
max norm relative error is plotted against viscosity (see Fig. 19). The "exact solution"
is based on a computation using N 80 modes and k and third-order SBDF.
CNLF was not included because the theoretical stability restriction is violated.
This can be easily seen because the linear advection-diffusion equation has eigenvalues
(ian- n2 2) for eigenfunctions e inx. From these eigenvalues we know that the
stability restriction is k < -, which is violated initially because lu(x, 0)1 1.
As expected, third-order SBDF has the smallest error of any of the methods when
stable. The stable region, however, is smaller than that for CNAB or SBDF. CNAB
and MCNAB are once again very similar in behavior with the modified version being
marginally better. SBDF appears to allow the largest stable time step when 0.01.
Further refinement of the time step to k leaves third-order SBDF as the
method of choice over the entire interval. Such a refinement may be unnecessary in
this example because the error is less than 1% for a step size which is four times
larger.
4.2.2. Chebyshev spectral collocation. Because the solution to the problem
is periodic and antisymmetric we know that u(+l, t) 0 for all t. Using this fact, we
solve (36) subject to the homogeneous Dirichlet boundary conditions u(:l,t) 0,
using a pseudospectral Chebyshev collocation scheme. The Gauss-Chebyshev grid
xi cos[(2i- 1)u/2N], i-1,...,N

is used for collocation points.
Similar to the Fourier case, the max norm relative errors are evaluated for several
IMEX schemes using k 1/40 and N 40. As can be seen from Fig. 20, the results
10-I
20. Chebyshev spectral collocation for Burgers equation.
are qualitatively similar to those of the Fourier case. SBDF performs particularly well
for smaller viscosities. From the theory, it has the widest stability region for large
la] and thus is best able to accommodate the rapidly growing eigenvalues associated
with Chebyshev collocation.
Both Chebyshev and Fourier collocation can be affected by aliasing [7]. We con-
sider aliasing for the Fourier case next.
4.2.3. Aliasing in pseudospectral methods. Aliasing occurs when nonlinear
terms produce frequencies that cannot be represented in the basis, and thus contribute
erroneously to lower frequencies. For instance, in Example 3, the Fourier sine mode
sin(rex) when explicitly evaluated in the convection term produces the contribution
sin(rex)
Osin(mx) m__ sin(2mx).
Ox 2
If 2m is greater than the number of Fourier sine basis functions N, this frequency
cannot be represented correctly and aliasing occurs. We now proceed to demonstrate
that this behavior can plague poorly spatially resolved computations when applying
weakly damping IMEX schemes.
We compute solutions for the model Burgers equation (36) subject to periodic
boundary conditions and the initial conditions of Example 3 modified to read
u(x, 0) 0.98 sin(2x) + HF(x)
where
HF(x) 0.01 sin(61rx) + 0.01 sin(62x).
We use Fourier sine collocation with N 64 basis functions, and integrate to time
t 2 with a viscosity u 0.1.
TABLE
k CNLF ’CNA’B’ MCNAB SB..DF..I .3ra orderSBDF
4--!-1’ ’l ":92 .57’5 "l .0056 .0072’ .0045

192,
.060’ .022 .0013
’.00079.’11 .0010
The value of the approximate solution at time t k is obtained using first-order

SBDF with the same step size. To represent the type of high frequency information
that could be produced during a computation, we add HF(x) to the solution at t k.
This ensures that high frequency information remains after the strongly damping first-
order SBDF step. Subsequent steps are then computed using the relevant second-order
scheme. For third-order SBDF, the value at t 2k is also needed. For the purpose
of demonstrating aliasing effects, this value is computed using CNAB since we wish
to retain most of the added high frequencies.
The max norm relative errors for several IMEX schemes are computed by compar-
1
ing results with those for SBDF using N 128 and k In the third-order case
results are compared with those for third-order SBDF using N 128 and k 24N"
These errors are summarized in Table 1.
For the case k we note that the error for CNAB is far greater than for
MCNAB, SBDF, or third-order SBDF. Using CNAB, a nonaliased computation using
the 2/3’s rule [7] and k results in a relative error of less that 10 -3. This
nonaliased result, along with its aliased counterpart, are plotted in Fig. 21. The
"error" curve in this figure is that of the aliased CNAB. From the figure it is clear
that the main component of the error is proportional to sin(rx). This low frequency
mode is not part of the exact solution and must be an aliasing artifact from high
frequency components. It is interesting to note that even after 128 time steps the
numerical solution is still plagued by high frequency components which have not yet
decayed. A similar study of CNLF reveals that it suffers from aliasing as well.
Resorting to a smaller time step k makes a substantial improvement in the
solution for CNAB and CNLF. Even so, the aliasing error for CNLF is sufficiently
large that further refinement is likely required.
We conclude that use of a strongly damping scheme such as SBDF, MCNAB, or
third-order SBDF gives an inexpensive method to reduce aliasing in poorly resolved
computations. Application of weakly damping schemes like CNAB and especially
CNLF may necessitate undesirably small time steps to produce the high frequency
decay needed to prevent aliasing in an aliased computation. Alternatively, weakly
damping schemes may be used in conjunction with an antialiasing technique such
as the 3/2’s rule or the 2/3’s rule [7]. However, these antialiasing techniques either
increase the expense of the computation or produce a severe loss of high frequency
information at each time step. rn
Remark. It may be argued that our addition of the high frequency components to
the numerical solution is artificial. However, this allows us to study the propagation
of errors related to these effects in simple examples within a controlled environment.
In general, high frequency solution and error components are generated at each time
step by nonlinear terms, forcing terms, and boundary effects (el. [6]). 3
4.3. Time-dependent multigrid in two spatial dimensions.
Example 4. The next model that we consider is the two-dimensional convection-
diffusion problem
(37) ut + (u. V)u

.00050
.00040
.00030
.00020
.00010
-.OOOLO
-.00020
-.00030
.00040
-.00050
-1.0
FI(. 21. Aliased and non-aliased computations for CNAB.
where u (u, v). This model incorporates some of the major ingredients of the
two-dimensional Navier-Stokes equations. Indeed, it is often being solved as part of
projection schemes for incompressible flows (see [13]). We carry out our computations
on the square ft --[0,1] x [0,1] and consider periodic boundary conditions and the initial
conditions
u(x, y, O) sin[2-(x + y)],

v(x, y, O) sin[2r(x + y)].
For spatial discretization we use standard second-order centered differences as
before. For IMEX schemes the convection term, (u. V)u is handled explicitly and
the diffusion term vAu is handled implicitly. This treatment yields a positive definite,
symmetric, sparse linear system to solve at each time step. Such systems are solved
etIiciently using a multigrid algorithm, the components of which are outlined next.
To solve the implicit equations we apply a full multigrid (FMG), full approxima-
tion scheme (FAS) algorithm [5]. Rather than applying the standard algorithm at
each time step, we use the ideas of [6]. Smoothing is accomplished using red-black
Gauss-Seidel. This relaxation technique is chosen because it has a very good smooth-
ing rate. Prolongation is accomplished using bilinear interpolation and restriction
by full weighting. The standard FMG cycle is modified so that the first coarse grid
correction is performed before any fine grid relaxation. Based on the assumption that
the increment between time steps is smooth, fine grid relaxation is most effective af-
ter the smoothness of the increment is accounted for (see [6]). Interpolation for the
FMG step is bilinear. In [6] it is argued that for time-dependent problems higher-
order interpolation only decreases high frequency error, which tends to dissipate in
parabolic systems anyway. This suggestion is utilized here; indeed experiments with
cubic interpolation did not produce any reduction in the number of multigrid itera-
tions to achieve a given tolerance. Another recommendation of [6] is to avoid the final
smoothing pass at each time step in order to reduce aliasing from red-black Gauss-
Seidel relaxation. Aliasing is not a major source of error in Example 4, however, so
this suggestion was not utilized.
For this problem we use a residual test with a tolerance TOL to determine the
number of fine grid iterations to perform at each time step. Next we show that the
choice of time-stepping scheme can affect the number of fine grid iterations to achieve
a given residual tolerance.
The model problem (37) was approximated using step sizes h 1-5-g and k
0.00625 and residual tolerance TOL=0.003. After the first time step, high frequency
information
HF(z, y) 0.005 cos[2r(64x + 63y)]
was added to each of u and v to represent the type of high frequency information that
might be produced during a computation.
For several second-order IMEX schemes, the average number of fine grid iterations
at each time step was computed. The result using V-cycles is graphed in Fig. 22.
Strongly damping schemes such as SBDF and MCNAB require roughly one iteration
per time step. Weakly damping schemes required far more effort to solve the implicit
equations accurately, because lingering high frequency components necessitate more
work on the finest grid. This is evident since CNAB requires more than two iterations
per time step and CNLF required three.
Using a W-cycle improves the relative efficiencies of CNLF and CNAB. Even for
these cases, however, nearly twice the number of fine grid iterations were required
than for more strongly damping schemes such as SBDF and MCNAB. See [22] for a
plot.
Even for a smaller viscosity u 0.02 and a much coarser mesh h
performance of CNLF suffers. It uses about 30% more iterations than SBDF or
,
the
MCNAB to achieve the desired tolerance. This result uses TOL=.009 and is plotted
in Fig. 23.
Thus we conclude that use of a poorly damping IMEX scheme such as CNAB
or especially CNLF can necessitate extra iterations on the finest grid in multigrid
solutions to the implicit equations for small mesh Reynolds numbers R < 2. For large
mesh Reynolds numbers R > 2, this effect was not observed. []
5. Conclusions and recommendations. IMEX scheInes are not a universal

cure for all problems. It is not difficult to imagine situations where a fully implicit
or fully explicit scheme is preferable. However, these schemes can be very effective in
many situations, some of which are depicted in this paper.
It may be important to choose an IMEX scheme carefully. The usual CNAB can
be significantly outperformed by other IMEX schemes. Based on observations in this
and previous sections we provide a few guidelines for selecting IMEX schemes for
convection-diffusion problems.
5.1. Finite differences. For the finite difference (finite element, finite volume)
case, begin by determining an estimate for the mesh Reynolds number R ah where
represents viscosity, a represents characteristic speed, and h represents grid spacing.
For problems where R >> 2 application of CNLF, or a third- or fourth-order
scheme is reasonable. Third- and fourth-order SBDF can be applied to these problems,
Number
of
iterations CNLF
per CNAB
time
step MCNAB
SBDF
10) (1
Second-order method (3’, c).
FIG. 22. Multigrid V-cycle iterations, 0.03, h 128
since a portion of the -axis is included in the stability region even though these
methods were selected primarily for their large u properties. Of the methods we have
considered in detail, CNLF has the mildest time step restriction and is nondissipative
while third-order SBDF is the most dissipative. All other second-order schemes should
be avoided in this case. Explicit schemes should also be considered for these problems.
For R > 2 use of CNLF or third-order SBDF appears appropriate. CNLF has the
mildest time step restriction, but accuracy concerns could make third-order SBDF
competitive. Avoid use of SBDF in this case. For problems of this type, a study to
determine when explicit schemes are competitive would be interesting.
For R 2 the theory predicts that many second-order IMEX schemes have similar
time step restrictions. A study to determine the method with the smallest truncation
error would be useful in this case. For greater accuracy, third-order SBDF appears to
be more useful than fourth-order SBDF, since its time step restriction is less severe. If
strong decay of high frequency spatial modes is a desirable characteristic then CNLF
should be avoided.
For R < 2, use of SBDF permits the largest stable time steps. The modified
CNAB scheme, MCNAB, can also be applied to problems of this type, although its
time step restriction is somewhat stricter, third-order SBDF is recommended when
high accuracy is needed. Numerical experiments in 4.3 demonstrate that in multi-
grid solutions to the implicit equations, application of a strongly damping method is
prudent. MCNAB or SBDF should be useful in such problems. Avoid use of CNLF
in this case.
5.2. Spectral methods. Because the eigenvalues for spectral methods are dif-
ferent than for finite differences as is their meaning (see [24]), we cannot expect the
stability time step restrictions from 3 to hold quantitatively. For this reason, an
822 W. M. ASCHER, S. J. RUUTH, AND B. T. R. WETTON
Number
of
iterations CNLF
per CNAB MCNAB SBDF
time
step
(0,1) (5, 0) (1,) (1,0)

Second-order method (, c).
FIG. 23. Multigrid W-cycle iterations, u 0.02, h
-
32"
analysis of the linear advection-diffusion equation for Chebyshev and Fourier spectral
methods would be interesting. However, this is outside the scope of this paper.
Numerical experiments for the Burgers equation were made for small to moderate
mesh Reynolds numbers. For these problems, CNLF has a very severe time step
restriction. These computations also suggest that SBDF has the mildest time step
restriction of the methods considered. This result is particularly pronounced in the
Chebyshev collocation case.
The large mesh Reynolds number case for spectral collocation was not considered.
In this case, a comparison of the relative eificiencies of IMEX schemes and fully explicit
schemes would be interesting.
Third-order SBDF’appears to be an eicient method for problems where high
accuracy is needed.
In problems where aliasing occurs, a strongly damping scheme, such as SBDF,
MCNAB, or third-order SBDF can be used to inexpensively reduce aliasing. Applica-
tion of a weakly damping scheme such as CNAB or CNLF in poorly spatially resolved
aliased computations should be avoided.
REFERENCES
[1] L. ABIA AND J. SANZ-SERNA, The spectral accuracy of a fully-discrete scheme for a nonlinear
third-order equation, Computing, 44 (1990), pp. 187-196.
[2] C. BASDEVANT, M. DEVILLE, P. HALDENWANG, J. LACROIX, J. OUAZZANI, l:{. PEYRET, P. OR-
LANDI, AND A. PATERA, Spectral and finite difference solutions of the Burgers equation,
Comput. Fluids, 14 (1986), pp. 23-41.
[3] J. BOYD, Chebyshev &: Fourier Spectral Methods, Springer-Verlag, Berlin, New York, 1989.
[4] M. BRACHET, D. MEIRON, S. ORSZAG, B. NICKEL, t{. M. ORE, AND V. FRISCH, Small-scale
structure of the Taylor-Green vortex, J. Fluid Mech., 130 (1983), pp. 411-452.
[5] A. BRANDT, Guide to multigrid development, in Multigrid Methods, W. Hackbusch and U. Trot-
tenberg, eds., Springer-Verlag, Berlin, 1982, pp. 220-312.
[6] A. BRANDT AND J. GREENWALD, Parabolic multigrid revisited, in Multigrid Methods III, W.
Hackbusch and U. Trottenberg, eds., Birkhuser-Verlag, Basel, Boston, 1991, pp. 143-154.
[7] C. CANUTO, M. HUSSAINI, A. QUARTERONI, AND T. ZANG, Spectral Methods in Fluid Dynamics,
Springer-Verlag, Berlin, New York, 1987.
[8] M. CROUZEIX, Une mdthode multipas implicite-explicite pour l’approximation des ’equations
d’dvolution paraboliques, Numer. Math., 35 (1980), pp. 257-276.
[9] C. GEAR, The automatic integration of stiff ordinary differential equations, in Proc. IFIP
Congress 1968, A. Morrell, ed., North Holland Publishing Company, Amsterdam, 1969,
pp. 187-193.
[10] B. GUSTAFSSON AND P. LOTSTEDT, Fourier analysis of multigrid methods for general systems
of PDEs, Math. Comp., 60 (1993), pp. 473-493.
[11] W. HACKBUSCH, Multi-Grid Methods and Applications, Springer-Verlag, Berlin, New York,
1985.
[12] E. HAIRER, S. NORSETT, AND G. WANNER, Solving Ordinary Differential Equations I, Springer-
Verlag, Berlin, New York, 1987.
[13] C. HIISCH, Numerical Computation of Internal and External Flows, Wiley, New York, 1988.
[14] A. JAMESON, Computational transonics, Comm. Pure Appl. Math., XLI (1988), pp. 507-549.
[15] G. KARNIADAKIS, M. ISRAELI, AND S. ORSZAG, High-order splitting methods for the incom-
pressible Navier-Stokes equations, J. Comput. Phys., 97 (1991), pp. 414-443.
[16] J. KIM AND P. MOIN, Application of a fractional-step method to incompressible Navier-Stokes
equations, J. Comput. Phys., 59 (1985), pp. 308-323.
[17] H. KREISS, Numerical methods for solving time-dependent problems for partial differential
equations, Lecture Notes, Univ. de Montreal, 1978.
[18] P. MEEK AND J. NORBURY, Two-stage, two-level finite difference schemes for nonlinear
parabolic equations, IMA J. Numer. Anal., 2 (1982), pp. 335-356.
[19] J. MURRAY, Mathematical Biology, Springer-Verlag, Berlin, New York, 1989.
[20] R. PEYRET AND T. TAYLOR, Computational Methods for Fluid Flow, Springer-Verlag, Berlin,
New York, 1983.
[21] S. RUUTH, Implicit-explicit methods for reaction-diffusion problems, J. Math. Biol., submitted.
[22] Implicit-Explicit Methods for Time-Dependent PDE’s, master’s thesis, Ins. Appl. Math.,
Univ. of British Columbia, Vancouver, 1993.
[23] J. STRIKWERDA, Finite Difference Schemes and Partial Differential Equations, Brooks/Cole,
Pacific Grove, CA, 1989.
[24] L. TREFETHEN, Lax-stability vs. eigenvalue stability of spectral methods, in Numerical Methods
for Fluid Dynamics III, K. Morton and M. Bains, eds., Clarendon Press, Oxford, 1988,
pp. 237-253.
[25] J. VAPAH, Stability restrictions on second-order, three level finite difference schemes for
parabolic equations, SIAM J. Numer. Anal., 17 (1980), pp. 300-309.
[26] R. VAR(A, Matrix Iterative Analysis, Prentice-Hall, Englewood Cliffs, N J, 1962.

Partial: Implicit-Explicit Methods For Time-Dependent Differential Equations Ascher, Wetton

Uploaded by

Copyright:

Available Formats

Partial: Implicit-Explicit Methods For Time-Dependent Differential Equations Ascher, Wetton

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Partial: Implicit-Explicit Methods For Time-Dependent Differential Equations Ascher, Wetton

Uploaded by

Copyright:

Available Formats

SIAM J. NUMER. ANAL.

1995 Society for Industrial and Applied Mathematics

IMPLICIT-EXPLICIT METHODS FOR TIME-DEPENDENT

PARTIAL DIFFERENTIAL EQUATIONS*

AMS subject classifications. 65J15, 65M20

1. Introduction. Various methods have been proposed to integrate dynamical

be because the Jacobian of f(u) is nonsymmetric and nondefinite and an iterative

(3) Ut aUx + ’Uxx

In 4 we perform a variety of numerical experiments using finite difference and

(p 1)! (p- 1)! + (p- 1)!

j’-I j--1 j--1

It is not difficult to prove the following theorem [22].

aDIU + D2U, 1 <_ <_ M.

2i i3jxi + aixi, j 1,...,

(8) (oj,/3j) ( 2[cs(2rjh)-l]’a sin(2rjh)

For notational convenience, we write

and pj is constant for -j simple and a polynomial in n otherwise. Clearly, stability

and we restrict 0 <_ /_< 1 to prevent large truncation error.

(3) un+ un kf(un) + l]]gg(u n+l).

X n+l (0, )X n, where (c, )

The stability region is thus {(c,/): I(c,)l _< 1}.

I( 0,)1 11+ ik] v/l+ k2 2 > 1.

[9(un+l) + 9(un)+ g(n-l

Finally, setting (7, c) (1, 0) yields

(7+c ) +(1-7- c)r +

(20) c= for 7 > O,

high frequency modes can be weak, tending to 1 as a

FIG. 6. Time step restriction for various where c /.

FIG 7. Time step restriction for various y where c

FIG. 8. Time step restriction for various "y where c O.

[23]), providing no damping of high frequency error components anywhere. To obtain

consider higher-order schemes.

Each term is divided by 2

(30) 0o 6(7o2 7o) 0,

We begin by examining if it is possible to arrive at a stable scheme for any fixed

minimum of T,0,c if 0.4661. Using these parameters, the scheme simplifies to

small time steps are needed for moderate or large

3.4. Fourth-order IMEX schemes and a comparison. From 2, the general

4f(un) 6f(u n-l) + 4f(un-2) f(un-3) + Pg(un+l).

(z) (25) (4 + + (3 + (4)

100 CNLF // O( CNLF CNLF

FIc. 15. Time step restrictions for various IMEX schemes.

apply to more complicated problems. In addition, strongly damping schemes are

FIG. 16. Large viscosity behaviour for various IMEX methods.

FIG. 1T. Small viscosity example for various IMEX methods.

As in the previous two sections, we use second-order centered differences to ap-

u(x, n) sin(2rx), n 0, 1, 2,...,

(36) u + UUx UUxx

FIG. 18. Burgers equation for various for v 0.01.

4.2.1. Fourier spectral collocation. Since the problem of Example 3 is pe-

FIG. 19. Fourier spectral collocation for Burgers equation.

xi cos[(2i- 1)u/2N], i-1,...,N

20. Chebyshev spectral collocation for Burgers equation.

4--!-1’ ’l ":92 .57’5 "l .0056 .0072’ .0045

The value of the approximate solution at time t k is obtained using first-order

(37) ut + (u. V)u

FI(. 21. Aliased and non-aliased computations for CNAB.

u(x, y, O) sin[2-(x + y)],

HF(z, y) 0.005 cos[2r(64x + 63y)]