Exponential Integrators for Large Systems of Differential Equations.
We study the numerical integration of large stiff systems of differential equations by methods that use matrix--vector products with the exponential or a related function of the Jacobian. For large problems, these can be approximated by Krylov subspace methods, which typically converge faster than those for the solution of the linear systems arising in standard stiff integrators. The exponential methods  also offer favorable properties in the integration of differential equations whose Jacobian has large imaginary eigenvalues. We derive methods up to order 4 which are exact for linear constant-coefficient equations. The implementation of the methods is discussed. Numerical experiments with reaction-diffusion problems and a time-dependent Schrdinger equation are included.
Introduction
. The idea to use the exponential function of the Jacobian in a
numerical integrator is by no means new, but it has mostly been regarded as rather
impractical. Since the mid-eighties, Krylov subspace approximations to the action of
the matrix exponential operator have, however, been found to be useful in Chemical
Physics [16, 20, 22] and subsequently also in other fields [6, 8, 9, 21, 24, 29]. On the
numerical analysis side, the convergence of such Krylov approximations was studied
in [4, 5, 13, 26]. It was shown in [13], and previously in [4] for the symmetric case,
that Krylov approximations to exp(-A)v converge substantially faster than those for
the solution of linear systems v, at least unless a good preconditioner
is available. Such linear systems arise in the numerical integration of stiff differential
equations by standard integrators. For large problems, their solution often dominates
the computational work.
For nonlinear differential equations, the exponential of the Jacobian combined
with Krylov approximations has previously been used in generalizations of Adams-
type multistep methods in [8]. On the other hand, the use of matrix exponentials has
for a long time been prominent in the exponential fitting literature, see e.g. [7, 17],
and [1, 2] as recent examples.
In this paper, we study new numerical methods for the integration of large stiff
systems of nonlinear initial value problems
The methods proposed here use matrix-vector multiplications '(-A)v, where A is the
Jacobian of f , - is related to the step size, and This choice allows
us to obtain methods that are exact for constant-coefficient linear problems
Mathematisches Institut, Universit?t T-ubingen, Auf der Morgenstelle 10, D-72076
T-ubingen, Germany. E-mail: marlis@na.uni-tuebingen.de, lubich@na.uni-tuebingen.de,
hubert@na.uni-tuebingen.de
We remark that Krylov subspace approximations to '(-A)v converge about as fast as
those for exp(-A)v, see [13].
Potential advantages for exponential integrators can thus originate from two different
sources: Computing '(-A)v can be less expensive than solving
and the exponential integration method itself may behave more favorably than standard
integrators. The latter case occurs in particular for mildly nonlinear differential
equations whose Jacobian has large imaginary eigenvalues, e.g., wave equations,
Schr-odinger equations, flexible mechanical systems, and oscillatory electric circuits.
Standard stiff integrators either damp high frequencies or map them to one and the
same frequency (or nearly so) in the discretization, neither of which may be desirable.
In Section 2 we give some simple methods of order 2 that are exact for (1.2) or
for linear second-order differential equations. They include new symmetric methods,
which appear useful for long-time integration of conservative, time-reversible problems.
In Section 3 we consider a class of methods that would reduce to explicit Runge-Kutta
methods if were replaced by '(z) j 1, and to Rosenbrock-
Wanner methods for We give order conditions, both for exact and
inexact Jacobian, and derive sufficient and necessary conditions to ensure that (1.2)
is solved exactly.
In Section 4 we extend the methods to differential-algebraic systems. We derive
order conditions up to order 3 for such problems and for singularly perturbed systems.
In Section 5 we construct methods of classical order 4 which are exact for (1.2)
and have further favorable properties when applied to stiff problems. In particular,
we use a reformulation that reduces the computational work for the Krylov processes.
Section 6 deals with implementation issues. Important topics are how to take into
account the computational work and storage requirements of the Krylov process in
the step size control, and when to stop the Krylov process.
Based on the considerations of Sections 5 and 6, we have written a code exp4,
which can be obtained via anonymous ftp from na.uni-tuebingen.de in the directory
pub/codes/exp4.
In Section 7 we describe numerical experiments with this code for reaction-diffusion
problems and for a Schr-odinger equation with time-dependent Hamiltonian, which
show both the scope and the limitations of using Krylov approximations in exponential
integrators.
In a final section we discuss conclusions and perspectives for the methods proposed
in this article.
We will describe the methods only for autonomous differential equations (1.1). For
non-autonomous problems y the methods should be applied to the extended
formally autonomous system obtained by adding the trivial differential equation t
The methods are then exact for linear differential equations whose inhomogeneity is
linear in t.
2. Simple methods of order 2.
2.1. The exponentially fitted Euler method. The prototype exponential
method, which seems to have appeared repeatedly in the literature under various
disguises, is
(2.
where h is the step size,
The method is of order 2, and exact for linear differential equations (1.2).
2.2. A symmetric exponential method. For long-time integration of conservative
time reversibility is an important property. A symmetric method of
order 2 is given by the two-step formula
with given by (2.2). The method is exact for linear problems (1.2)
provided that the starting values y 0 and y 1 are exact. The method can be viewed
as a generalization of the explicit midpoint rule, to which it reduces for
characteristic roots of the method applied to y are e h- and \Gamma1 which shows that
the method is A-stable. The oscillatory error component (\Gamma1) n can be eliminated by
taking the average of two successive values, (y n as an approximation to
2.3. A cosine method for second-order differential equations. We now
consider
For the linear problem
the exact solution satisfies
with the entire function
This motivates the scheme
with which is a symmetric method of order 2 for (2.4). Because of (2.6),
the scheme is exact for linear problems (2.5).
Derivative approximations that are exact for (2.5) are obtained via
where
\Gammaz
\Gammaz
3. Higher-order exponential one-step methods: Order conditions and
stability. In this section we study a general class of exponential integration methods
introduced in [13]. Starting with y 0 as an approximation to y(t 0 ), an approximation
to computed via
s
Here are the coefficients that
determine the method. The internal stages can be computed one after the
other, with one multiplication by '(flhA) and a function evaluation at each stage. The
scheme would become an explicit Runge-Kutta method for '(z) j 1 (and
a Rosenbrock-Wanner method for the choice As in the exponential
Euler method (2.1), we choose instead the function (2.2).
3.1. Order conditions when using the exact Jacobian. Our aim now is to
construct higher-order methods. The order conditions for the exponential methods
can be derived similarly to Rosenbrock-Wanner methods, see, e.g., [12, Section IV.7].
Therefore, we only state the conditions here. For abbreviation we define
Theorem 3.1. An exponential method (3:1)-(3:3) with is of order p
iff
s
for elementary differentials - up to order p. Here, \Phi j (-) and the polynomials P - (fl)
are listed in Table 3:1 for p - 5.
The only difference to the order conditions for Rosenbrock-Wanner methods is in
the polynomials P - (fl).
Theorem 3.2. The method (3:1)-(3:3) is exact for linear differential equations
(1:2), iff for all
These conditions can be fulfilled if fl is the reciprocal of an integer. Then only a finite
number of these conditions are needed. The others are satisfied automatically because
for sufficiently large n, both sides of (3.5) then vanish.
Proof. For the linear problem (1.2), both the exact and the numerical solution
depend analytically on h. Since only the elementary differentials f , f
are nonvanishing for (1.2), it thus suffices to show that their order conditions are
given by (3.5). Like for Rosenbrock methods, one obtains that they are of the form
(3.
Elementary
differential -
k;l ff jk ff jl 1=3
f 000 (f; f; f)
k;l;m ff jk ff jl ff jm 1=4
f (4) (f; f; f; f)
f 000 (f 0 f; f; f)

Table
Order conditions for exponential methods up to order 5
where P n\Gamma1 is a polynomial of degree at most n \Gamma 1 which depends on the choice of
' but not on the method coefficients. It remains to show that P n\Gamma1 (fl) is given by
the right-hand side of (3.5). If then the method applied to (1.2) is just a
Runge-Kutta method with coefficients fi jk and weights b i . From the order conditions
for Runge-Kutta methods, we thus have
The exponential Euler method (2.1) is a one-stage method (3:1)-(3:3) with b
Obviously,
method. Since we already know that the exponential Euler method
is exact for (1.2) we conclude from (3.6) that
Similarly, two consecutive steps of the exponential Euler method with step size h=2
can be viewed as one step of a two-stage method (3:1)-(3:3) with 1=2. For such
a method, (3.7) is valid for n ? 2. As before, we conclude from (3.6) that
Continuing this argument for 3; steps of the exponential Euler method with step
sizes
Elementary
differential -
Af
k;l ff jk ff jl 1=3
k;l ff jk ff kl 1=6
k;l ff jk fl kl \Gammafl =4
k;l
AAf
k;l

Table
Order conditions for exponential W-methods up to order 3
It follows that P n\Gamma1 (fl) is given by the right-hand side of (3.5).
3.2. Order conditions for inexact Jacobians. One may also want to use the
method with an approximate Jacobian A. This requires further restrictions on the
method parameters. For order 3 the conditions are given in Table 3.2. They are the
same as for W-methods, see [12, p. 124], except for different polynomials in fl.
If the first five conditions of Table 3.2 are satisfied, then the method is of order 3
when for the analogous situation in W-methods.
3.3. Stability. When the method is exact for linear differential equations, it is
trivially A-stable. Much more can then in fact be shown about stability, including
the practical situation where '(flhA)v is computed only approximately. Consider a
perturbed method (3.1)-(3.3) applied to the linear problem (1.2):
Here, ffi i is a perturbation at the ith stage, and e y 0 is a perturbed starting value.
Subtracting from the unperturbed scheme yields for the error "
It is easy to see that
s
is a polynomial of degree k with p k whose coefficients are products
of In particular, when the numerical range of A is contained in the left half-
plane, then we have the stable error recurrence
s
The stability analysis could be extended to nonlinear problems y
Lipschitz-bounded g, to singularly perturbed problems, and to nonlinear parabolic
problems in a similar way to what has been done for Rosenbrock methods, cf. [10, 18,
30].
4. Exponential methods for differential-algebraic and singularly perturbed
problems. As with Rosenbrock-Wanner methods [12, Section VI.3], the
method (3.1)-(3.3) can be extended to differential-algebraic equations
z invertible
by applying it formally to the singularly perturbed differential equation
and letting ffl ! 0 in the scheme. This will give us the following method:
l i
I 0
\Gammag
z g y I
I \Gammaf z g \Gamma1
z
I
y g z
l i
z g y , where the Jacobians are evaluated at (y
and v i are defined by
z 0
l j
Finally we set
z 0
s
l i
The derivation uses the auxiliary assumption that the eigenvalues of g z have negative
real part. The Jacobian of (4.2) is block diagonalized by
I +O(ffl) fflf z g \Gamma1
z
\Gammag
z g y +O(ffl) I +O(ffl)
viz.,
ffl g z +O(1)
Since
'( flh
z
the method (3.1)-(3.3) applied to (4.2) reads
l i
z
l i
We note that
I 0
I 0
ffl I
I \Gammaf z g \Gamma1
z
I
lead to the above method (4.3).
Remark. The matrix B need not be computed when one uses Krylov methods
to approximate '(flhB)u. Matrix vector multiplications with B are cheap when the
action of g \Gamma1
z is inexpensive to compute. For example, this is the case in constrained
mechanical systems, cf. [12, p. 542],
a
Here, q and v are position and velocity variables, respectively, a is acceleration and -
represents the Lagrange multipliers. In this system, g z corresponds to
In suitable multibody formulations, linear equations with this matrix can be solved in
an amount of work proportional to the dimension.
When the exponential method is exact for linear differential equations with constant
inhomogeneity, then method (4.3)-(4.5) is exact for linear differential-algebraic
equations
with constant matrices F y ; F z ; G y ; G z (G z invertible) and constant vectors b; c. Apart
from a direct calculation, this may be seen as follows: When the eigenvalues of G z
have negative real part, the exactness is again obtained by letting ffl ! 0 in the singularly
perturbed problem, which is solved exactly by the method. From this case, the
exactness in the general situation of invertible G z follows by analytical continuation.
In general, the application of this scheme to differential-algebraic equations results
in an order reduction to order 2, unless the method coefficients satisfy additional
conditions.
Theorem 4.1. The method (4:3)-(4:5) is convergent of order 3 for the differential-algebraic
equation (4:1) if it satisfies the order conditions of Table 3:1 up to order 3,
(3:5), and in addition
where [! jk
The additional order condition is the same as for Rosenbrock methods applied to
(4.1) [12, p. 446]. Instead of giving a cumbersome formal proof of the theorem, we make
the reappearance of condition (4.6) for exponential methods plausible as follows: Like
the order conditions of Section 3, also the differential-algebraic order conditions are of
the same form as for Rosenbrock methods, but possibly with different right-hand sides
involving fl. We know that the theorem is valid for z). The appearance
of the ! jk is related only to the term (flhg z ) \Gamma1 in (4.3), which is independent of '.
The terms ff jk are also unrelated to '(flhB). Therefore, the condition remains the
same as for Rosenbrock methods.
The differential-algebraic order condition (4.6) is important not only for differential-algebraic
systems but also for stiff differential equations. For example, the third-order
error bound of Rosenbrock methods for singularly perturbed problems (4.2) in Theorem
of [10] can be shown to be valid also for exponential methods.
5. Construction of fourth-order methods.
5.1. Reduced methods. We recall that one step of the exponential method
evaluated in the form (3.1)-(3.3) contains s multiplications of '(flhA) with a vector.
Since this vector is different in each of these s steps, the approximation with a Krylov
subspace method requires the construction of bases of s Krylov spaces with respect
to the same matrix A but with different vectors. This turns out to be prohibitively
expensive. One may think of exploiting techniques for solving linear systems with
multiple right-hand sides [25, 27], but in our experiments the savings achieved were
minor. Therefore, we will present an alternative formulation of the method.
A key point for the construction of efficient methods is that one can compute
recursively from '(z):
Once we have computed '(flhA), we can thus compute '(jflhA)v for any integer j ? 1
with the expense of matrix vector multiplications.
The recurrence (5.1) is equally useful for the more interesting case where '(jflhA)v
is approximated by Krylov methods. The Krylov subspace approximation is of the
is the matrix containing the Arnoldi (or Lanczos) basis of the
mth Krylov subspace with respect to A and v, and Hm is the orthogonal (oblique)
projection of A to the mth Krylov subspace, which is an m \Theta m upper Hessenberg
(block tridiagonal, respectively) matrix. Further, e 1 is the first m-dimensional unit
vector.
The iteration number m is typically very small compared to the dimension of
the matrix A, so that the matrix '(flhH m ) can be computed quite cheaply (see x6 for
details). Then the recurrence (5.1) can be used to compute '(jflhH m )e 1 by performing
matrix vector multiplications with the small matrices Hm and '(flhH m ). If we denote
the identity matrix of dimension m by I m , then
We can exploit the recurrences (5.1) by reformulating the method. For this we
introduce auxiliary vectors
s
Note that for corresponds to a first-degree Taylor expansion of f
around y 0 . Hence the vectors d i are usually small in norm and would vanish for linear
f . With (3.4) and (5.3) we have
s
Because of (5.1) we can choose fi kl such that for
All the coefficients fi kl are uniquely determined by (5.4). In order to apply the recurrence
formulas (5.1) in (5.4) we further choose
which gives
This reduces the number of f-evaluations and of evaluations of '(flhA) by a factor of
n compared to the general scheme (3.1)-(3.3). This is particularly important when
this reduced method is combined with a Krylov process for approximating '(flhA)v
since in this case we need to compute a basis of a new Krylov space only at every
nth intermediate step. Moreover, since the vectors d i are usually small in norm, the
Krylov approximation of '(iflhA)d nj+1 (j - 1) typically takes only few iterations to
achieve the required accuracy. The cost for building up the Krylov space of A with
respect to the vector f(y 0 ) thus dominates the computational cost.
We note finally that we can reorganize the computations in (5.4) as
~
and we can use the values ~ k l in (3.2) and (3.3), with
appropriately modified weights:
~
and
5.2. Methods of order 4. Next we show that the reduced scheme proposed
above still allows the construction of higher-order methods. Here, we concentrate on
start with a 3-stage method for uses 2 function
evaluations per step. The parameters fi kl satisfying (5.4) are given by
To fulfill the conditions for order 4, there remain two free parameters ff 3;1 ; ff 3;2 , and the
weights 4. The order conditions from Table 3.1 have a unique solution
This yields the scheme
On k 3 we have omitted the tilde corresponding to (5.5). This method is of order 4,
and exact for linear differential equations (1.2). However, it is only of first order when
used with inexact Jacobian and of second order when applied to DAEs. Moreover,
it is impossible to construct an embedded method of order 3, which makes it hard
to perform a reliable estimation of local errors for step size control. The only cheap
variant is to use the exponential Euler method (2.1), which is only of order 2 and thus
tends to overestimate the local error.
The method (5.7) with embedded (2.1) is however of interest as a very economical
method in situations where the time step is not restricted by accuracy, but only by
the convergence of the Krylov process for computing '(hA)f(y 0 ). We note that k 3
is usually well approximated in a low dimensional Krylov space, because d 3 is often
much smaller in norm than f(y 0 ).
A more sophisticated method can be constructed with using
3 function evaluations per step. The parameters for (5.4) are given by
With these parameters fi kl , all the order conditions (3.5) for linear problems are satisfied
automatically for
For our method we choose to evaluate the function f at both end points and at
the middle of the time interval, i.e.,X
The solution is obtained by first solving the order condition up to order 4 from Table
3.1. The equations for f 00 (f; f) and f 0 f immediately yield b
the conditions for f f) and (4.6) result in a linear system with
four equations for the unknowns b j , 7. This system has the unique solution
which also satisfies the second order W-condition.
From the equation for f we obtain b It remains to fulfill the equation for
further we satisfy the third-order W-condition for f 0 f 0 f in order to
obtain order 3 when the approximation to the Jacobian is O(h) close to the true Ja-
cobian. This yields ff 2. We
still have some freedom so that we can solve the fifth-order conditions for f 000 (f 0 f; f; f)
and f 00 (f 0 f 0 f; f ). This gives ff
other fifth-order conditions can be satisfied, we now minimizeX
7;j
which yields ff
This construction gives us the following method:
Again, we have omitted the tilde on k used in (5.5). We summarize the
properties of this method in a theorem.
Theorem 5.1. The scheme (5:8) is of order 4 for differential equations (1:1),
and exact for linear differential equations (1:2). It converges of order 3 for differential-algebraic
equations (4:1) and to smooth solutions of singularly perturbed problems (4:2)
uniformly for ffl - h 2 . For differential equations (1:1), it is of second order when used
with inexact Jacobian, and of order 3 when the approximation to the Jacobian is O(h)
close to the true Jacobian.
The method satisfies three of the order-5 conditions. The residuals of the other
conditions appear to be rather small, the largest one being 0:1.
Although the scheme (5.8) is a 7-stage method, it requires only three function
evaluations. When using Krylov approximations, the computational cost is dominated
by computing k 1 . As discussed before, the reason is that k 2 , k 3 , k 5 , and k 6 can be
computed recursively from (5.1) or the more stable recurrence (6.2) below, and that
k 4 to k 7 are typically well approximated in very low dimensional Krylov subspaces,
because d 4 and d 7 are usually much smaller in norm than f(y 0 ). For these reasons, and
because of its superior theoretical properties, we prefer (5.8) to a "standard" 3-stage
fourth-order scheme of type (3.1)-(3.3).
5.3. Embedded methods. We have constructed two embedded methods with
different properties for the scheme (5.8). The first one is of order 3 for differential
equations (1.1) and differential-algebraic equations (4.1), and exact for linear equations
(1.2). Solving the third-order conditions of Table 3.1 and condition (4.6), and choosing
gives the embedded scheme
This method does not satisfy the fourth-order conditions, except that for f 0 f 0 f 0 f . It
is however only of order 1 as a W-method, i.e., when used with inexact Jacobian.
The second embedded method is of order two as a W-method. It is not exact for
linear differential equations (1.2), and it does not satisfy the third-order conditions of

Table

3.1. It reads
e
5.4. Dense output. Like for Runge-Kutta and Rosenbrock methods, a continuous
numerical solution defined via
s
with polynomials b i satisfying b i This approximation is of order
iff
s
for all elementary differentials - of order ae - p, see [12, p. 452].
For the 3-stage method (5.7) a continuous numerical solution of order 3 is given
by
For the 7-stage method (5.8) a continuous numerical solution of order 3, which is also
of order 3 for differential-algebraic equations and of order 2 when used with inexact
Jacobian, is given by
The actual computation uses
s
are defined as in (5.6) and ~
k i are the k i from (5.8).
6. Implementation issues.
6.1. Step size control. The step size control for the scheme (5.8) uses the two
embedded methods proposed in Section 5.3. As an estimate for the local error, we
choose the minimum of the local error estimates of these two methods. A step size
selection strategy due to Gustafsson, see [12, p. 31-35] and the Radau5 code [12,
p. 550ff], then yields a new step size proposal h err .
However, if Krylov subspace methods are used to approximate the matrix exponential
operator, then in addition to the local error estimate it is necessary to take the
work and storage requirements of the Krylov process into account. We propose the
following strategy: First choose a suitable "window" [-; M ] for the number of Krylov
steps m required in the approximation of k (recall that in reduced methods,
the overall work of the Krylov processes is dominated by this first Krylov process).
In this window we choose a desirable number of Krylov steps m opt . We preserve the
actual step size h of the integration method whenever m 2 [-; M ]. If m ? M the new
Krylov step size is reduced until the required accuracy is achieved with an m 2 [-; M ].
in two consecutive steps, we set
where we have found 1=3 as a reasonable value in our numerical experiments. It
also turned out that a more drastical enlargement of the step size is possible if m is
very small for more than two consecutive steps. For example we used
in the last j time steps.
Finally we choose the new step size as
6.2. Savings from previous steps. The scheme may reuse the Jacobian of a
previous time step as an approximation to the actual Jacobian. This is done if the
local error of the embedded method (5.10) is acceptable and in addition h kry - h err ,
i.e., the step size is determined by the Krylov process.
Further savings can be achieved if the Jacobian A and the step size h are the same
as in the previous time step. We then write
If f(y n ) is close to f(y n\Gamma1 ), then the initial vector for the Krylov process is small in
norm and thus the Krylov process becomes less expensive.
6.3. Stopping criterion for the Krylov method. We need to decide when
the Krylov approximation (5.2) is to be considered sufficiently accurate. Since exact
errors are inaccessible, the stopping criterion in the iterative solution of linear systems
usually based on the residual
instead of the error of the mth iterate
For Galerkin-type methods like FOM and BiCG, the residual vectors can be computed
from
where hm+1;m is the (m+1;m) entry of Hm+1 , and [ \Delta ] m;1 denotes the (m; 1)-entry of a
matrix. Using Cauchy's integral formula, the error of the mth Krylov approximation
to '(-A)v can be written as
Z
where \Gamma is a contour enclosing the eigenvalues of -A and -Hm , cf. [13]. Thus, the
can be interpreted as a linear combination of errors e m (-) of linear systems.
Replacing e m (-) by r m (-) in this formula, we get a generalized residual
Z
which can be computed at no additional cost. This suggests to use ae m instead of the
unknown error ffl m in the stopping criterion. The use of ae m was also proposed by Saad
[26], who used a different derivation that is plausible only for small k-Ak.
In the scheme (5.8), the Krylov approximations to k j are multiplied by the step
size h. It is therefore reasonable to stop the iteration if
tol is the weighted norm used in the integrator:
with are the given absolute
and relative error tolerances.
In our numerical experiments we found that (6.1) is on the safe side, but sometimes
rather pessimistic. Then it may pay off to apply an idea attributed to Shampine in [12,
p. 134], which consists in using a smoothed residual instead of the true
residual. Since solving a linear system with coefficient matrix prohibitively
expensive when A is large, one can perform a smoothing in the m-dimensional subspace
and use
instead of (6.1) for m - 5, say. For smaller m, this criterion may be overly optimistic
when -A has large norm.
6.4. Computation of '(-Hm ). To reduce the computational costs, we evaluate
figures in an index sequence, e.g., m 2 f1; 2; 3; 4; 6; 8; 11; 15; 20;
48g. This sequence is chosen such that the computation of '(-Hm ) is about as
expensive as the total of the previously computed '(-H j ), since the computation of
If A is Hermitian, then Hm is Hermitian tridiagonal. In this case, one can simply
diagonalize Hm .
In the non-Hermitian case, we suggest to use Pad'e approximation similarly to the
third method described in [19] to compute the matrix exponential. Here, the matrix
is first scaled by a factor of 2 \Gammak such that k2 \Gammak -Hm k ! 1=2. Then we evaluate the
approximation to '(z) for the scaled matrix:
26 z
Next, '(-Hm ) is computed recursively from '(2 \Gammak -Hm ) by applying the following
coupled recurrences:
This recurrence is stable for all z in the left half-plane, whereas (5.1) becomes unstable
for large jzj because of the multiplication with z.
Alternatively, in the non-Hermitian case, one can use a formula due to Saad [26,
Section 2.3]:
exp
exp(-Hm
This appears favorable when the dimension m is not too small.
7. Numerical experiments. We have implemented the method (5.8) with (and
approximations in a Matlab code exp4. The program is written in
the format used in the Matlab ODE suite [28], which is available via anonymous ftp
on ftp.mathworks.com in the pub/mathworks/toolbox/matlab/funfun directory. The
code exp4 can be obtained from na.uni-tuebingen.de in the pub/codes/exp4 directory.
A C version of exp4 is also available from this ftp site.
7.1. A reaction-diffusion equation with nonstiff chemical reaction: the
Brusselator. To illustrate the behavior of the exponential integrator with Krylov
approximations to '(flhA)v in the transition from a nonstiff to a stiff problem, we
have chosen the two-dimensional Brusselator [11, pp. 248ff]:
@t
together with Neumann boundary conditions
@n
@n
and initial conditions
The Laplacian is discretized on a uniform 100 \Theta 100 grid by central differences, so
that the dimension of the resulting ODE problem is 20; 000. The eigenvalues of the
discretized Laplacian lie between \Gamma80; 000 and zero. We present numerical experiments
with three different values of the diffusion coefficient
which mark the transition from a nonstiff to a stiff problem. The solution of
the problem for shown in the movie on pp. 250ff in [11].
In Figs. 7.1-7.3 we show work-precision diagrams for our exponential integrator
exp4, and for the explicit Runge-Kutta integrator ode45 from the Matlab ODE
suite [28], which is based on a fifth-order method of Dormand and Prince [3]. The
vertical axis shows the error at the end point the horizontal axis gives the
exp4
error
flops
Fig. 7.1. Brusselator for
flops
error
exp4
Fig. 7.2. Brusselator for
exp4
error
flops
Fig. 7.3. Brusselator for
number of flops. The markers \Theta for exp4 and ffi for ode45 correspond to the error
While the computational work of the
nonstiff integrator increases drastically with growing ff, the performance of exp4 is
considerably less affected.
7.2. A reaction-diffusion equation with stiff chemistry: the Robertson
example. The following example shows the behavior of the exponential integrator
for a very stiff problem. We consider the Robertson reaction [12, pp. 3f] with one-dimensional
diffusion:
together with Neumann boundary conditions u
The diffusion coefficient is chosen as . The second spatial derivative
is discretized on a uniform grid with points. In this problem, the stiffness
originates from the reaction terms. We have chosen such a small problem because
we intend to illustrate the influence of the Krylov approximation procedure to the
performance of the integrator. In Fig. 7.4 we show the step sizes as a function of time
in a double logarithmic scale with and without Krylov approximation of '(flhA)v. As
this example has only dimension 90, '(flhA) can here be computed by diagonalization.
For comparison, we also show the step size of the explicit integrator ode45 and the
stiff integrator ode15s from the Matlab ODE suite [28], which uses a variant of a
method. All the methods have been run with the same tolerances atol=rtol=
\Gamma6 . It is seen that in this example the step size is always limited by the Krylov
process. The step size restriction does not appear very severe on the considered
time interval. Similar step size sequences are obtained for the Krylov-approximated
exponential method for higher-dimensional versions of the problem. However, the
limits of the Krylov approach show up when the integration is continued to very long
times. There, the step size remains essentially on the level seen at the right-most part
of Fig. 7.4. It has been observed that this behavior is largely due to roundoff error
effects.
7.3. A Schr-odinger equation with time-dependent potential. As an example
of a problem whose Jacobian has large imaginary eigenvalues we consider, following
[23], the one-dimensional Schr-odinger equation for
with the Hamiltonian
This equation models an atom/molecule interacting with a high intensity CW laser.
The parameter values used were 100. The initial value was /(x;
exp4 without Krylov
exp4 with Krylov
step size
time
Fig. 7.4. Step sizes versus time for the Robertson example.
which corresponds to the eigenstate of the unforced harmonic oscillator
to the lowest energy level. Semi-discretization in space is done by a pseudospectral
method with Fourier modes on the space interval x 2 [\Gammaa; a] for
with periodic boundary conditions. This leads to the non-autonomous linear system
of differential equations for
diag
Here, y j (t) is an approximation to /(x
N , FN is the discrete
Fourier-transform operator, and
a diag(0;
The Jacobian is full but matrix-vector multiplications are obtained with O(N log N)
flops using FFT. In Fig. 7.5 we show the work-precision diagram at
the standard nonstiff and stiff solvers ode45 and ode15s from the Matlab ODE suite,
and for a Matlab implementation of Hairer and Wanner's [12] Radau5 implicit Runge-Kutta
code. The codes were used with tolerances atol=rtol=
The surprisingly good behavior of the stiff integrators ode15s and radau5 is due
to the following matrix-free implementation: In the simplified Newton iterations the
Jacobian was approximated by i=2F
so that the linear systems could be
solved in O(N log N) operations using FFT. Therefore, the computational cost per
time step was essentially the same as for an explicit method. Using the full Jacobian
would make the implicit methods completely inefficient for this problem. We note,
however, that the performance of the versions with the simplified Jacobian deteriorates
when the parameters - and - increase.
The exponential code exp4 is clearly superior to the explicit integrator ode45.
For accuracy requirements more stringent than 10 \Gamma4 , Fig. 7.5 shows an advantage for
exp4
error
flops
Fig. 7.5. Work-precision diagram for the Schr-odinger equation.
exp4 also with respect to the implicit methods in their optimized versions discussed
above. This stems from the fact that exp4 is able to take much larger time steps than
the other integrators.
In computations with Schr-odinger equations with time-independent Hamiltonian,
the use of Chebyshev approximations to the matrix exponential operator is very popular
[16]. We therefore also implemented a version of exp4 where the Arnoldi process is
replaced by a Chebyshev approximation. In our numerical experiments the Chebyshev
version needed about twice as many flops as the Arnoldi-based implementation.
[We thank S. Gray for pointing out references [16] and [23].]
8. Conclusions and perspectives. In this paper we have introduced new integration
methods which use matrix-vector multiplications with the exponential of the
Jacobian. In particular, we have studied Rosenbrock-like exponential methods. Since
a straightforward implementation of these methods is computationally expensive, we
have identified a subclass of "reduced" methods which are reformulated such that they
allow for an efficient implementation. Two promising fourth-order methods, which are
exact for linear constant-coefficient problems, have been given in formulas (5.7) and
(5.8). The method (5.8), which offers superior properties at slightly higher cost per
time step, has been implemented with Krylov subspace approximations to the matrix
exponential operator in a code exp4. This implementation requires only function
evaluations and matrix-vector multiplications with the Jacobian.
Numerical experiments and theoretical considerations indicate that exponential
integrators are highly competitive for the following problem classes of large systems
of initial-value problems:
Mildly stiff problems (e.g., reaction-convection-diffusion problems with nonstiff
reaction terms): The most efficient traditional methods are explicit integrators which
are used despite stability restrictions of the step size. For the special case where the
eigenvalues of the Jacobian are on the negative real axis, the high-stage Runge-Kutta-
Chebyshev methods of van der Houwen and Sommeijer [14, 31] are known to be very
successful. Here, the theory in [5, 13] and [31] tells us that the number of necessary
matrix-vector multiplications with the Jacobian in Krylov iterations for exponential
methods and the number of function evaluations in Runge-Kutta-Chebyshev methods
needed to attain stability are both of the magnitude of
hkAk. However, the Krylov
methods take advantage of clustered eigenvalues and of vectors with small components
in some eigendirections. There is no restriction to problems with eigenvalues near the
real axis for the exponential methods with Krylov approximations, and much larger
time steps than with standard explicit Runge-Kutta methods (such as the Dormand-
Prince methods) can be taken.
Stiff problems (e.g., reaction-diffusion problems with stiff reaction terms): For
high-dimensional systems, the standard approach is to use implicit methods (such as
BDF or Radau) where the linear systems are solved iteratively with the help of a
hopefully good and cheap preconditioner. If - and only if - an efficient preconditioner
is available, those methods are clearly favorable over the exponential methods proposed
here, since it is not known how to precondition the iterative computation of
the matrix exponential operator. Due to the superlinear error reduction of the Krylov
approximations to the matrix exponential, exponential methods are often competitive
even without a preconditioner. We hope that future developments will allow to effectively
use ideas of preconditioning in the computation of the exponentials and hence
further enlarge the range of stiff problems on which exponential methods are efficiently
applicable.
Highly oscillatory problems (e.g., wave equations, Schr-odinger equations, elasto-
dynamics, and oscillatory electric circuits): Here, the proposed exponential methods
are able to resolve high frequencies to the required error tolerance without the severe
time step restrictions of standard schemes. Time step restrictions of an often milder
type still occur because of nonlinear effects and because of limitations of the iteration
number in the Krylov process. The latter are less severe when the eigenvalues of
the Jacobian are clustered. The good resolution of high frequencies with exponential
methods is in contrast to usual implicit integrators used with large time steps, which
either damp high frequencies or map them to one and the same frequency (or nearly
so) in the discretization.
It will be interesting to see how the new methods perform in real-life scientific
problems.



--R

Sufficient conditions for uniformly second-order convergent schemes for stiff initial-value problems
A new efficient numerical integration scheme for highly oscillatory electric circuits
A family of embedded Runge-Kutta
bounds in the simple Lanczos procedure for computing functions of symmetric matrices and eigenvalues
Krylov subspace approximations of eigenpairs and matrix functions in exact and computer arithmetic
Krylov methods for the incompressible Navier-Stokes equations
Verallgemeinerte Runge-Kutta Verfahren zur L-osung steifer Differentialgleichungen
A method of exponential propagation of large systems of stiff nonlinear differential equations
Efficient solution of parabolic equations by Krylov approximation methods
of Rosenbrock methods for stiff problems studied via differential algebraic equations
Solving Ordinary Differential Equations I
Solving Ordinary Differential Equations II
On Krylov subspace approximations to the matrix exponential operator
On the internal stability of explicit m-stage Runge-Kutta methods for large values of m
Rosenbrock methods using few LU-decompositions
Propagation methods for quantum molecular dynamics
Generalized Runge-Kutta processes for stable systems with large Lipschitz con- stants
Linearly implicit time discretization of non-linear parabolic equations
Nineteen dubious ways to compute the exponential of a matrix
New approach to many-state quantum dynamics: The recursive- residue-generation method
Applications of the Lanczos algorithm
Unitary quantum time evolution by iterative Lanczos reduction
The solution of the time dependent Schr-odinger equation by the (t

On the Lanczos method for solving symmetric systems with several right hand sides
Analysis of some Krylov subspace approximations to the matrix exponential operator

The Matlab ODE suite
Expokit: Software package for computing matrix exponentials

Explicit Runge-Kutta methods for parabolic partial differential equations
--TR

--CTR
Xin-Yuan Wu , Jian-Lin Xia, New vector forms of elemental functions with Taylor series, Applied Mathematics and Computation, v.141 n.2-3, p.307-312, 5 September
S. Koikari, An error analysis of the modified scaling and squaring method, Computers & Mathematics with Applications, v.53 n.8, p.1293-1305, April, 2007
R. Weiner , B. A. Schmitt , H. Podhaisky, Parallel 'Peer' two-step W-methods and their application to MOL-systems, Applied Numerical Mathematics, v.48 n.3-4, p.425-439, March 2004
Philip W. Livermore, An implementation of the exponential time differencing scheme to the magnetohydrodynamic equations in a spherical shell, Journal of Computational Physics, v.220 n.2, p.824-838, January, 2007
A. Schmitt , R. Weiner, Design, analysis and testing of some parallel two-step W-methods for stiff systems, Applied Numerical Mathematics, v.42 n.1, p.381-395, August 2002
X.-Y. Wu , J.-L. Xia , F. Yang, A dynamic method for weighted linear least squares problems, Computing, v.68 n.4, p.375-386, September 2002
Ya Yan Lu, Computing a matrix function for exponential integrators, Journal of Computational and Applied Mathematics, v.161 n.1, p.203-216, 1 December
Jrg Wensch, Krylov-ROW methods for DAEs of index 1 with applications to viscoelasticity, Applied Numerical Mathematics, v.53 n.2, p.527-541, May 2005
Marlis Hochbruck , Alexander Ostermann, Exponential Runge-Kutta methods for parabolic problems, Applied Numerical Mathematics, v.53 n.2, p.323-339, May 2005
Z. Jackiewicz , H. Podhaisky , R. Weiner, Construction of highly stable two-step W-methods for ordinary differential equations, Journal of Computational and Applied Mathematics, v.167 n.2, p.389-403, 1 June 2004
Jrg Wensch , Fernando Casas, Extrapolation in lie groups with approximated BCH-formula, Applied Numerical Mathematics, v.42 n.1, p.465-472, August 2002
Tobin A. Driscoll, A composite Runge-Kutta method for the spectral solution of semilinear PDEs, Journal of Computational Physics, v.182 n.2, p.357-367, November 2002
M. Lintner, The eigenvalue problem for the 2D laplacian in H-matrix arithmetic and application to the heat and wave equation, Computing, v.72 n.3-4, p.293-323, May 2004
S. Krogstad, Generalized integrating factor methods for stiff PDEs, Journal of Computational Physics, v.203 n.1, p.72-88, 10 February 2005
M. Caliari , M. Vianello , L. Bergamaschi, Interpolating discrete advection-diffusion propagators at Leja sequences, Journal of Computational and Applied Mathematics, v.172 n.1, p.79-99, 1 November 2004
F. Carbonell , J. C. Jimenez , R. Biscay, A numerical method for the computation of the Lyapunov exponents of nonlinear ordinary differential equations, Applied Mathematics and Computation, v.131 n.1, p.21-37, 10 September 2002
Roger B. Sidje, Expokit: a software package for computing matrix exponentials, ACM Transactions on Mathematical Software (TOMS), v.24 n.1, p.130-156, March 1998
Hvard Berland , Brd Skaflestad , Will M. Wright, EXPINT---A MATLAB package for exponential integrators, ACM Transactions on Mathematical Software (TOMS), v.33 n.1, p.4-es, March 2007
M. A. Botchev , D. Harutyunyan , J. J. W. van der Vegt, The Gautschi time stepping scheme for edge finite element discretizations of the Maxwell equations, Journal of Computational Physics, v.216 August 2006
Elena Celledoni , Arieh Iserles , Syvert P. Nrsett , Bojan Orel, Complexity theory for lie-group solvers, Journal of Complexity, v.18 n.1, p.242-286, March 2002
M. Tokman, Efficient integration of large stiff systems of ODEs with exponential propagation iterative (EPI) methods, Journal of Computational Physics, v.213
high-dimensional differential equations;numerical integrator;matrix exponential;krylov subspace methods