On Krylov Subspace Approximations to the Matrix Exponential Operator.
Krylov subspace methods for approximating  the action of matrix exponentials are analyzed in this paper. We derive error bounds via a functional calculus of Arnoldi and Lanczos methods that reduces the study of Krylov subspace approximations of functions of matrices to that of linear systems of equations. As a side result, we obtain error bounds for Galerkin-type Krylov methods for linear equations, namely, the biconjugate gradient method and the full orthogonalization method. For Krylov approximations to matrix exponentials, we show superlinear error decay from  relatively small iteration numbers onwards, depending on the geometry of the numerical range, the spectrum, or the pseudospectrum. The convergence to exp$(\tau A)v$ is faster than that of corresponding Krylov methods for the solution of linear equations $(I-\tau A)x=v$, which usually arise in the numerical solution of stiff ordinary differential equations (ODEs). We therefore propose a new class of time integration  methods for large systems of nonlinear differential equations which use Krylov approximations to the exponential function of the Jacobian instead of solving linear or nonlinear systems of equations in every time step.
Introduction
. In this article we study Krylov subspace methods for the approximation
of exp(-A)v when A is a matrix of large dimension, v is a given vector,
scaling factor which may be associated with the step size in a time integration
method. Such Krylov approximations were apparently first used in Chemical
Physics [20, 22, 17] and were more recently studied by Gallopoulos and Saad [10, 24];
see also their account of related previous work. They present Krylov schemes for exponential
propagation, discuss the implementation, report excellent numerical results,
and give some theoretical error bounds. As they also mention, these bounds are however
too pessimistic to explain the numerically observed error reductions. Moreover,
their error bounds do not make evident that - or when and why - Krylov methods
perform far better than standard explicit time stepping methods in stiff problems. A
further open question concerns the relationship between the convergence properties of
Krylov subspace methods for exponential operators and those for the linear systems
of equations arising in implicit time integration methods. In the present paper we
intend to clear up the error behavior.
When we wrote this paper, we were unaware of the important previous work by
Druskin and Knizhnerman [3, 4, 14, 15] who use a different approach to the analysis.
We will comment on the relationship of some of their results to ours in a note at the
end of this paper.
Our error analysis is based on a functional calculus of Arnoldi and Lanczos methods
which reduces the study of approximations of exp(-A)v to that of the corresponding
iterative methods for linear systems of equations. Somewhat oversimplified, it
may be said that the error of the mth iterate for exp(-A)v behaves like the minimum,
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
taken over all ff ? 0, of e ff multiplied with the error of the mth iterate for the solution
of by the same Krylov subspace method. This minimum is usually
attained far from especially for large iteration numbers. Unless a good preconditioner
for I \Gamma -A is available, the iteration for exp(-A)v converges therefore faster
than that for v. We do not know, however, of a way to "precondition"
the iteration for exp(-A)v .
Gallopoulos and Saad showed that the error of the mth iterate in the approximation
of exp(-A)v has a bound proportional to k-Ak m =m!, which gives superlinear
convergence for m AE k-Ak. In many cases, however, superlinear error decay begins
for much smaller iteration numbers. For example, we will show that for symmetric
negative definite matrices A this occurs already for m -
k-Ak, whereas for skew-Hermitian
matrices with uniformly distributed eigenvalues substantial error reduction
begins in general only for m near k-Ak. We will obtain rapid error decay for m- k-Ak
also for a class of sectorial, non-normal matrices. Convergence within the required tolerance
ensures that the methods become superior to standard explicit
time stepping methods for large systems. For m - k-Ak, our error bounds improve
upon those of [10] and [24] typically by a factor 2 \Gammam e \Gammack- with a c ? 0. The analysis
explains how the error depends on the geometry of critical sets in the complex
plane, namely the numerical range of A for Arnoldi-based approximations, and the
location of the spectra or pseudospectra of A and the Krylov-Galerkin matrix Hm for
both Lanczos- and Arnoldi-based approximations. In our framework, it is also easily
seen that clustering of eigenvalues has similar beneficial effects in the Krylov subspace
approximation of exp(-A)v as in the iterative solution of linear systems of equations.
As mentioned above, exp(-A)v can often be computed faster than
Krylov subspace methods. This fact has implications in the time integration of very
large systems of ordinary differential equations arising, e.g., in many-particle simulations
and from spatial discretizations of time-dependent partial differential equations.
It justifies renewed interest in ODE methods that use the exponential or related functions
of the Jacobian instead of solving linear or nonlinear systems of equations in
every time step. Methods of this type in the literature include the exponential Runge-Kutta
methods of Lawson [18] and Friedli [8], the adaptive Runge-Kutta methods of
Strehmel and Weiner [27] in their non-approximated form, exponentially fitted methods
of [5], and the exponential multistep methods of [9].
In the last section of this paper, we propose a promising new class of "exponen-
tial" integration methods. With Krylov approximations, substantial savings can be
expected for large, moderately stiff systems of ordinary differential equations which
are routinely solved by explicit time-stepping methods despite stability restrictions of
the step size, or when implicit methods require prohibitively expensive Jacobians and
linear algebra.
The paper is organized as follows: In Section 2, we describe the general framework
and derive a basic error bound for the Arnoldi method. In Section 3, this leads to
specific error bounds for the approximation of exp(-A)v for various classes of matrices
A. Lanczos methods are studied in Section 4, which contains also error bounds for
BiCG and FOM. In Section 5, we introduce a class of time-stepping methods for large
systems of ODEs which replace the solution of linear systems of equations by multiplication
with '(-A), where whose Krylov subspace approximations
converge as fast as those for exp(-A)v.
Throughout the paper, k \Delta k is the Euclidean norm or its induced matrix norm.
2. Arnoldi-based approximation of functions of matrices. In the sequel,
let A be a complex square matrix of (large) dimension N , and v 2 C N a given
vector of unit length, 1. The Arnoldi process generates an orthonormal basis
of the Krylov space
matrix Hm of dimension m (which is the upper left part of its successor Hm+1 ) such
that
where e i is the ith unit vector in R m . By induction this clearly implies, as noted in
[3, Theorem 2] and [24, Lemma 3.1],
A standard use of the Arnoldi process is in the solution of linear equations [23], where
one approximates
when - is not an eigenvalue of A, and hopefully not of Hm . The latter condition is
always satisfied when - is outside the numerical range
since (2.1) implies
m AVm and therefore
We now turn to the approximation of functions of A. Let f be analytic in a neighborhood
of F(A). Then
Z
where \Gamma is a contour that surrounds F(A). In view of (2.3), we are led to replace this
Z
so that we approximate
Such an approximation was proposed previously [22, 33, 3, 10], with different derivations

In practice, we are then left with the task of computing the lower-dimensional
expression f(Hm )e 1 , which for m - N is usually much easier to compute than f(A)v,
e.g., by diagonalization of Hm . The above derivation of (2.7) also indicates how to
obtain error bounds: Study the error in the Arnoldi approximation (2.3) of linear
systems and integrate their error bounds, multiplied with jf(-)j, for - varying along
a suitable contour \Gamma. This will actually be done in the present paper.
Our error bounds are based on Lemma 1 below. To prepare its setting, let E be a
convex, closed bounded set in the complex plane. Let OE be the conformal mapping
that carries the exterior of E onto the exterior of the unit circle fjwj ? 1g, with
We note that ae is the logarithmic capacity
of E. Finally, let \Gamma be the boundary curve of a piecewise smooth, bounded region G
that contains E, and assume that f is analytic in G and continuous on the closure of
G.
Lemma 1. Under the above assumptions, and if the numerical range of A is contained
in E, we have for every polynomial q m\Gamma1 of degree at most
Z
with is the length of the boundary curve @E of
E, and where d(S) is the minimal distance between F(A) and a subset S of the complex
plane. If E is a straight line segment or a disk, then (2:8) holds with
Remark. It will be useful to choose the integration contour dependent on m, in order
to balance the decay of OE \Gammam away from E against the growth of f outside E. For
functions f , such as the exponential function studied in detail below, this will
ultimately yield superlinear convergence. On the other hand, the liberty in choosing
the polynomial q m\Gamma1 will not materialize in the study of the exponential function.
Proof. (a) We begin by studying the error of (2.3) and consider a fixed
the moment. Our argumentation in this part of the proof is inspired by [25]. Using
m v, we rewrite the error as
with
. By (2.1) and the orthogonality of
Hence we have for arbitrary y
We note that
is a polynomial of degree - m with Conversely, for every such
(A)v is of the above form. To bound \Delta m , we recall kVm use
the estimates k(-I
which follow from [26, Thm.4.1] and (2.4). We thus obtain
for every polynomial p m of degree at most m with
(b) It remains to bound p m (A). Since
Z
we have
For the special case when E is a line segment, we have A of the form A
with a Hermitian B and complex coefficients ff; fi, so that then
When E is a disk jz \Gamma -j - ae, then p m (z) will be chosen as a multiple of (z
and an inequality of Berger (see [1, p. 3]) then tells us that
In all these cases we thus have
with M as stated in Lemma 1.
(c) To proceed in the proof, we use near-optimality properties of Faber polynomials.
These have been employed previously in analyses of iterative methods by Eiermann
[6] and Nevanlinna [21]. Let OE m (z) denote the Faber polynomial of degree m associated
with the region E. This is defined as the polynomial part of OE(z) m , i.e.,
1. We now choose the polynomial p m (z) with the
normalization
cf. [21, p.76]. A theorem of K-ovari and Pommerenke [16, Thm.2] provides us with the
bound
This implies max z2E jOE m
(d) Combining inequalities (2.10),(2.11), and (2.14) gives us
The proof is now completed by inserting this bound into the difference of formulas
(2.5) and (2.6) and taking account of (2.2).
Remark. Part (c) of the above proof, combined with Cauchy's integral formula, shows
that there exists a polynomial \Pi m\Gamma1 (z) of degree at most m \Gamma 1 such that
Z
where ffi is the minimal distance between \Gamma and E. This holds for \Pi
R
of (2.12). Polynomial approximation
bounds of this type are closely tied to Bernstein's theorem [19, Thm.III.3.19].
3. Approximation of the matrix exponential operator. In this section we
give error bounds for the Arnoldi approximation of e -A v for various classes of matrices
A. We may restrict our attention to cases where the numerical range of A is contained
in the left half-plane, so that ke -A k and, in view of (2.4), also are bounded by
unity. This assumption entails no loss of generality, since a shift from A to A
changes both e -A v and its approximation by a factor e - ff .
The same bounds as in Theorems 2 to 6 below (even slightly more favorable
are valid also for Krylov subspace approximations of '(-A)v, with
Theorem 2. Let A be a Hermitian negative semi-definite matrix with eigenvalues
in the interval [\Gamma4ae; 0]. Then the error in the Arnoldi approximation of e -A v, i.e.,
is bounded in the following ways:
Remark. It is instructive to compare the above error bounds with that of the conjugate
gradient method applied to the linear system which is given by
This bound becomes small for m AE p ae- but only with a linear decay.
Proof. We use Lemma 1 with Then the conformal mapping is
We start by applying the linear transformation to the
interval [\Gamma1; 1]. As contour \Gamma we choose the parabola with right-most point fl that is
mapped to the parabola \Pi given by the parametrization
This parabola osculates to the ellipse E with foci \Sigma1 and major
us the error bound
Z
Z
where \Phi(-+
1. The absolute value of \Phi(-) is constant along every ellipse
with foci \Sigma1. Since the parabola \Pi is located outside the ellipse E , we have along \Pi
Hence we obtain from (3.3)
Z 1e \Gammaae- ' 2
d'
s
Fig. 3.1. Errors and error bounds for the symmetric example
Moreover, we have r - e ff
2ffl with ff ? 0:96 for ffl - 1=2 (and ff ? 0:98 for ffl - 1=4).
Minimizing e 2ae- ffl\Gammam
2ffl with respect to ffl yields
Inserting this ffl in (3.4) results in the bound (with
which together with " m - 2 is a sharper version of (3.1). The condition ffl - 1=2 is
equivalent to m - 2ae- .
To obtain the bound (3.2), we note that 1 insert in (3.4)
which is close to the minimum for m AE ae- . This yields for m - 2ae-
which is a sharper version of (3.2).
Finally we remark, in view of the proof of Theorem 3 below, that the bounds (3.1)
and (3.2) are also obtained when \Gamma is chosen as a composition of the part of the above
parabola contained in the right half-plane and two rays on the imaginary axis.
To give an illustration of our error bounds we consider the diagonal matrix A with
equidistantly spaced eigenvalues in the interval [\Gamma40; 0] and a random unit vector v of
dimension 1001. Fig. 3.1 shows the errors of the approximation to exp(A)v and those
of the cg approximation to nearly a straight line. Moreover, the
dashed line shows the error bounds (3.5) and (3.6), while the dotted line corresponds
to 2k 1
which is the error bound of [24, Corollary 4.6] for symmetric, negative
semi-definite matrices A.
It is well known that Krylov subspace methods for the solution of linear systems of
equations benefit from a clustering of the eigenvalues. The same is true also for the
Krylov subspace approximation of e -A v. This is actually not surprising in view of the
Cauchy integral representations (2.5), (2.6). As an example of such a result, we state
the following theorem. This might be generalized in various directions for different
types of clusterings and different types of matrices, but we will not pursue this further.
Theorem 3. Let A be a Hermitian negative semi-definite matrix with eigenvalues
contained in f- 1 \Gamma4ae. Then the (m 1)st error " m+1 in the
Arnoldi approximation of e -A v is bounded by the right-hand sides of (3:1) and (3:2).
Proof. The result is proved by using the polynomial
instead of (2.12). The absolute value of the first factor is bounded by unity for z 2
[\Gamma4ae; 0] and Re- 0. Hence we obtain the same error bounds as in Theorem 2 with
replaced by m \Gamma 1.
For skew-Hermitian matrices A (with uniformly distributed eigenvalues) we cannot
show superlinear error decay for m ! ae- . The reason is, basically, that here the
conformal mapping OE maps the vertical line contour with jOE(-)j -
in the symmetric negative definite case we have jOE(-)j
ffl. This
behavior affects equally the convergence of Krylov subspace methods for the solution
of linear systems v. For skew-Hermitian A, there is, in general, no
substantial error reduction for m ! ae- , and convergence is linear with a rate like
Theorem 4. Let A be a skew-Hermitian matrix with eigenvalues in an interval on
the imaginary axis of length 4ae. Then the error in the Arnoldi approximation of e -A v
is bounded by
Proof. We use Lemma 1 with the conformal mapping is
After applying the linear transformation we choose the integration
contour as an ellipse with foci \Sigma1 and minor semiaxis
semiaxis is then
and the length of the contour is bounded by 2-a. In
addition, we have 1). The absolute value
constant along the ellipse. With Lemma 1, we get for the error
Fig. 3.2. Errors and error bounds for the skew-Hermitian example
with Inserting gives the stated error bound. A sharper
bound for ae- ? 1
2 , which is obtained by integrating over the parabola that osculates
to the above ellipse at the right-most point, reads
For a numerical illustration we choose the diagonal matrix A with 1001 equidistant
eigenvalues in [\Gamma20i; 20i], and a random vector v of unit length. Fig. 3.2 shows the
errors of the approximation to exp(A)v and those of the BiCG approximation to
nearly a straight line. The dashed line shows the error
bound (3.7), the dotted line corresponds to which is the bound given in
[10, Corollary 2.2].
Theorem 5. Let A be a matrix with numerical range contained in the disk jz
Then the error in the Arnoldi approximation of e -A v is bounded by
Proof. We use Lemma 1 with a circle with radius rae centered
at \Gammaae. Lemma 1 gives the bound
Setting gives the stated result.
The following is a worst-case example which shows nearly no error reduction for m -
ae- .
Example. Let A be the bidiagonal matrix of dimension N that has \Gamma1 on the diagonal
and +1 on the subdiagonal. The numerical range of A is then contained in the disk
The Arnoldi process gives as the m-dimensional version of A,
so that
The error vector thus contains the entries e \Gamma- k =k! for k ? m. The largest of these is
close to (2-) \Gamma1=2 by Stirling's
Similar to Theorem 2, the onset of superlinear convergence begins already for m- ae-
when F(A) is contained in a wedge-shaped set. In particular, consider the conformal
mapping
which maps the exterior of the unit disk onto the exterior of the
bounded sectorial set in the left half-plane
S ' has a corner at 0 with opening angle '- and is symmetric with respect to the real
axis.
Theorem 6. For some ae ? 0 and let the numerical range of A be contained
in ae \Delta S ' . Then the error in the Arnoldi approximation of e -A v is bounded by
ae-
with
1\Gamma' . The constants C and c ? 0 depend only on '.
Proof. In the course of this proof, C denotes a generic constant which takes on
different values on different occurrences. After the transformation -=ae, we use
choose ffi such that r=(1
note that then ae . The
right-most point of the integration contour ae Hence we have from
For m AE (ae- ) ff , the right-hand side is minimized near
Inserting this ffl gives (3.8) with
ff ff
For any r - 2 and ffi ? 0, Lemma 1 gives
e /(r)ae-
which becomes (3.9) upon choosing
4. Lanczos-based approximation of functions of matrices. The Arnoldi
method unfortunately requires long recurrences for the construction of the Krylov
basis. The Lanczos method overcomes this difficulty by computing an auxiliary basis
which spans the Krylov subspace with respect to A   and w 1 . The
Lanczos vectors v j and w j are constructed such that they satisfy a biorthogonality
condition, or block biorthogonality in case of the look-ahead version [7, 28], i.e., Dm :=
m Vm is block diagonal. The look-ahead process ensures that Dm is well conditioned
when the index m terminates a block, which will be assumed of m in the sequel.
The Lanczos vectors can be constructed by short (mostly three-term) recurrences.
This results again in a matrix representation (2.1), but now with a block tridiagonal
m AVm . However, unlike the Arnoldi case, neither Vm nor Wm
are orthogonal matrices. It is usual to scale the Lanczos vectors to have unit norm,
in which case the norms of Vm and Wm are bounded by p
m. Since Hm is now an
oblique projection of A, the numerical range of Hm is in general not contained in
F(A). Variants of Lemma 1, which apply in this situation, are given in the following
two lemmas. For the exponential function, Lemmas 7 and 8 lead to essentially the
same error bounds as given for the Arnoldi method in Theorems 5 and 6, except for
different constants. In Theorems 2, 3, and 4, Arnoldi and Lanczos approximations
coincide.
The first lemma works with the ffl-pseudospectrum of A [32], defined by
Otherwise, the setting is again the one described before Lemma 1.
Lemma 7. If  then the error of the Lanczos approximation
of f(A)v is bounded by (2:8) with
Proof. The proof modifies the proof of Lemma 1. For the Lanczos process we have
I , and therefore
Noting
m v, we thus obtain
for every polynomial p m of degree - m with 1. By assumption we have that
the norms of both (-I are bounded by fl \Gamma1 for - 2 \Gamma. Using
further kp m (A)k - '(@E)=(2-ffl) \Delta max z2E jp m (z)j leads to
(4.
which in turn yields the estimate stated in the lemma.
For a diagonalizable matrix A we let
where X is the matrix that contains the eigenvectors of A in its columns.
The following lemma involves only the spectrum  (A) of A and uses once more
the setting of Lemma 1.
Lemma 8. Let A be diagonalizable and assume that  (A) ae E,  (H m
\Gamma. Then the Lanczos approximation of f(A)v
satisfies (2:8) with
, where ffi is the minimal
distance between  (A) and \Gamma.
Proof. The result follows from (4.1) along the lines of parts (c) and (d) of the proof
of Lemma 1.
Remarks. (a) It is known that in generic situations, extreme eigenvalues of A are well
approximated by those of Hm for sufficiently large m [34]. For a contour \Gamma that is
bounded away from  (A), one can thus expect that usually k(-I uniformly
bounded along \Gamma.
(b) Lemmas 7 and 8 apply also to the Arnoldi method, where
and kVm
(c) The convexity assumption about E can be removed at the price of a larger
factor M . For E a continuum containing more than one point, one can use instead of
inequality (2.13) the estimate in the lemma on pp. 107f. in Volume III of [19].
The proofs of Lemmas 1, 7, and 8 provide error bounds for iterative methods for the
solution of linear systems of equations whose iterates are defined by a Galerkin condition
(2.3). This gives new error bounds for the biconjugate gradient method, where the
Krylov basis is constructed via the Lanczos process, and for the full orthogonalization
method, which is based on the Arnoldi process. The proofs can be extended to give
similar error bounds also for the (quasi-) minimization methods QMR and GMRES,
see [13].
5. A class of integration methods for large systems of ODEs. In the
numerical integration of very large stiff systems of ordinary differential equations y
f(y), Krylov subspace methods have been used successfully for the solution of the
linear systems of equations arising in fully or linearly implicit integration schemes [11,
2, 25]. These linear systems are of the form A is the Jacobian
of f evaluated near the current integration point, h is the step size, and fl is a method
parameter. The attraction with Krylov subspace methods lies in the fact that they
require only the computation of matrix-vector products Aw. When it is convenient,
these can be approximated as directional derivatives Aw :
that the Jacobian A need never be formed explicitly. Our theoretical results as well
as computational experiments indicate that Krylov subspace approximations of e flhA v
or '(flhA)v, with
converge faster than the corresponding iterations for v, at least unless
a good preconditioner is at hand. This suggests the use of the following class of
integration schemes, in which the linear systems arising in a linearly implicit method
of Rosenbrock-Wanner type are replaced by multiplication with '(flhA). Starting
from y 0 - y(t 0 ), the scheme computes an approximation y 1 of
s
Here are the coefficients that determine the method. The
internal stages are computed one after the other, with one multiplication
by '(flhA) and a function evaluation at each stage. The simplest method of this type
is the well-known exponentially fitted Euler method
which is of order 2 and exact for linear differential equations y
A and b. It appears well suited as a basis for Richardson extrapolation. Here is another
example of such a method:
Theorem 9. The two-stage methods with coefficients
eter),
are of order 3. For arbitrary step
sizes, they provide the exact solution for every linear system of differential equations
with constant matrix A and constant inhomogeneity b.
Proof. Taylor expansion in h of the exact and the numerical solutions shows that
the order conditions up to order 3, which correspond to the elementary differentials
are given by
Here all sums extend from 1 to s, and we have set ff Cf. with the
order conditions for Rosenbrock methods in [12], p.116, which differ from the present
order conditions only in the right-hand side polynomials in fl.
For the right-hand side of the last order condition vanishes, and hence
this condition is automatically satisfied for every two-stage method with
With ff remaining three equations yield the stated
method coefficients. Direct calculation shows that the method applied to y
which is the claimed property.
Remarks. (a) With 3=4, the method satisfies the order condition fi 2 ff 3
which corresponds to the fourth-order elementary differential f 000 (f; f; f ). The order
conditions corresponding to f are satisfied independently of ff,
so that the order condition corresponding to f 00 (f; f 0 f) is then the only fourth-order
condition that remains violated.
(b) For non-autonomous problems y it is useful to rewrite the equation
in autonomous form by adding the trivial equation t taking the Jacobian
e
In particular, the method is then exact for every linear equation of the form y
tc, since this is rewritten as
y
c A
y
which is again a linear system with constant inhomogeneity.
An efficient implementation and higher-order methods are currently under investigation

Note added in the revised version. After finishing this paper we learned that
Druskin and Knizhnerman [3, 4] previously obtained an estimate similar to (3.5) for
the symmetric case, using a different proof. They give the asymptotic estimate
a
'-
a
a 3
with which they prove using the Chebyshev series expansion of the exponential
function. In an extension of this technique to the non-Hermitian case, Knizhner-
man [14] derived error bounds in terms of Faber series for the Arnoldi method (2.7).
He showed
k=m
are the Faber series coefficients of f and the exponent ff depends on the
numerical range of A. As one referee emphasizes, the Faber series approach could be
put to similar use as our Lemma 1. In fact, Leonid Knizhnerman showed to us in a
personal communication how it would become possible to derive a result of the type
of our Theorem 6 using (5.9). Our approach via Lemma 1 makes it more obvious to
see how the geometry of the numerical range comes into play. An example similar to
that after Theorem 5 is given in [15, x3]. We thank Anne Greenbaum and two referees
for pointing out these references and Leonid Knizhnerman for providing a commented
version of the Russian paper [14]. Error bounds via Chebyshev and Faber series, for
the related problem of approximating matrix functions by methods that generalize
semi-iterative methods for linear systems, were given by Tal-Ezer [29, 30, 31].

Acknowledgement

. We are grateful to Peter Leinen and Harry Yserentant for providing
the initial motivation for this work.



--R

Numerical Ranges of Operators on Normed Spaces and of Elements of Normed Algebras

Two polynomial methods of calculating functions of symmetric matrices
Krylov subspace approximations of eigenpairs and matrix functions in exact and computer arithmetic
Uniform Numerical Methods for Problems with Initial and Boundary Layers
On semiiterative methods generated by Faber polynomials
An implementation of the look-ahead Lanczos algorithm for non-Hermitian matrices
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
Iterative solution of linear equations in ODE codes
Solving Ordinary Differential Equations II
analysis of Krylov methods in a nutshell
Computation of functions of unsymmetric matrices by means of Arnoldi's method
bounds in Arnoldi's method: The case of a normal matrix

Propagation methods for quantum molecular dynamics
Generalized Runge-Kutta processes for stable systems with large Lipschitz con- stants
Theory of Functions of a Complex Variable
New approach to many-state quantum dynamics: The recursive- residue-generation method
Convergence of Iterations for Linear Equations
Unitary quantum time evolution by iterative Lanczos reduction
Krylov subspace methods for solving large unsymmetric linear systems
Analysis of some Krylov subspace approximations to the matrix exponential operator



Analysis of the look-ahead Lanczos algorithm
Spectral methods in time for hyperbolic problems
Spectral methods in time for parabolic problems
Polynomial approximation of functions of matrices and applications
Pseudospectra of matrices
An iterative solution method for solving f(A)x
A convergence analysis for nonsymmetric Lanczos algorithms
--TR

--CTR
P. Novati, An explicit one-step method for stiff problems, Computing, v.71 n.2, p.133-151, October
Vladimir Druskin, Krylov Subspaces and Electromagnetic Oil Exploration, IEEE Computational Science & Engineering, v.5 n.1, p.10-12, January 1998
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
C. Gonzlez , A. Ostermann , M. Thalhammer, A second-order Magnus-type integrator for nonautonomous parabolic problems, Journal of Computational and Applied Mathematics, v.189 n.1, p.142-156, 1 May 2006
Christian Lubich, A variational splitting integrator for quantum molecular dynamics, Applied Numerical Mathematics, v.48 n.3-4, p.355-368, March 2004
Marlis Hochbruck , Alexander Ostermann, Exponential Runge-Kutta methods for parabolic problems, Applied Numerical Mathematics, v.53 n.2, p.323-339, May 2005
S. Krogstad, Generalized integrating factor methods for stiff PDEs, Journal of Computational Physics, v.203 n.1, p.72-88, 10 February 2005
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
Serhiy Kosinov , Stephane Marchand-Maillet , Igor Kozintsev , Carole Dulong , Thierry Pun, Dual diffusion model of spreading activation for content-based image retrieval, Proceedings of the 8th ACM international workshop on Multimedia information retrieval, October 26-27, 2006, Santa Barbara, California, USA
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
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
Paolo Novati, A polynomial method based on Fejr points for the computation of functions of unsymmetric matrices, Applied Numerical Mathematics, v.44 n.1-2, p.201-224, January
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
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
James V. Lambers, Practical Implementation of Krylov Subspace Spectral Methods, Journal of Scientific Computing, v.32 n.3, p.449-476, September 2007
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
matrix exponential function;matrix-free time integration methods;arnoldi method;superlinear convergence;lanczos method;conjugate gradient-type methods;krylov subspace methods