B -SERIES AND ORDER CONDITIONS FOR EXPONENTIAL INTEGRATORS

. We introduce a general format of numerical ODE-solvers which include many of the recently proposed exponential integrators. We derive a general order theory for these schemes in terms of B -series and bicolored rooted trees. To ease the construction of speciﬁc schemes we generalize an idea of Zennaro [ Math. Comp., 46 (1986), pp. 119–133] and deﬁne natural continuous extensions in the context of exponential integrators. This leads to a relatively easy derivation of some of the most popular recently proposed schemes. The general format of schemes considered here makes use of coeﬃcient functions which will usually be selected from some ﬁnite dimensional function spaces. We will derive lower bounds for the dimension of these spaces in terms of the order of the resulting schemes. Finally, we illustrate the presented ideas by giving examples of new exponential integrators of orders 4 and 5.


Introduction.
Numerical integration schemes which use the matrix exponential go back all the way to Certaine [4], but there are also early papers by Lawson [15], Nørsett [20], Ehle and Lawson [6], and Friedli [7] to mention just a few.Recently there has been a revived interest in these schemes, in particular for the solution of nonlinear partial differential equations; see for instance [11,17,5,3,14,13].For a thorough review of the history of exponential integrators; see [16] and the references therein.The integrators found in these papers are derived in rather different ways, and they are formulated for different types of systems of differential equations.On this note, we consider the autonomous nonlinear system of ordinary differential equations Here L is a matrix and N (u) a nonlinear mapping.The order theory we consider is valid for a large class of exponential integrators, including the Runge-Kutta-Munthe-Kaas (RKMK) schemes [17], the commutator-free Lie group integrators [3], and those schemes of Cox and Matthews [5], as well as Krogstad [14] which reduce to classical Runge-Kutta schemes when L = 0.
We present the general format for integrators of (1.1) as a j r (hL) N j , r = 1, . . ., s (1.2) Here we assume that the functions a j r (z) and b r (z) are at least p times continuously differentiable at z = 0 for integration schemes of order p.
(b) Commutator-free, order 4 0 Table 1 gives the coefficient functions a j r (z) and b r (z) for the fourth order RKMK scheme introduced in [18] in this general format when applied to the problem (1.1) with an affine Lie group action, and the commutator-free scheme of order 4 from [3]; in both tables φ 0 (z) = (e z − 1)/z.
For deriving order conditions, we expand the coefficient functions in powers of z, where the sum may terminate with a remainder term.For the schemes we consider here, these functions are in fact all entire.If N (u) = 0 in (1.1), then any scheme in the above class will reproduce the exact solution in every step.Whereas if L = 0, the scheme (1.2)-(1.3)reduces to a classical Runge-Kutta method with coefficients a j r = α j,0 r and b r = β r,0 .This scheme is henceforth called the underlying Runge-Kutta scheme.We will always assume that c r = j α j,0 r , 1 ≤ r ≤ s.The schemes proposed by Friedli [7] closely resemble the format (1.2)-(1.3), the difference being that the coefficient functions a rj (resp.b r ) are evaluated in c r hL rather than in hL, thus a nontrivial discrepancy may occur whenever c r = 0.And even though Friedli explicitly requires that the functions a rj (z) and b r (z) be of the form his analysis holds also for the case of more general coefficient functions, so that the order conditions he obtains for p ≤ 4 are almost identical to those derived in section 2 here.However, the order theory presented here is general.We will discuss conditions on the coefficients α j,k r and β r,k under which the scheme (1.2)-(1.3)has order of consistency p for problems of the type (1.1).We will use the well known approach involving rooted trees; see, for instance, [9,2].The conditions we find will depend only on the first α j,k r for k ≤ p − 2 and on β j,k for k ≤ p − 1.On this note we will not address issues related to the behavior of the coefficient functions a j r (z) and b r (z) for large values of z.
In the recent paper [12], an order theory for explicit exponential integrators is presented and its application to semilinear parabolic problems is discussed.While classical or nonstiff order conditions are usually derived by assuming that a Lipschitz constant exists, one needs to account for the unboundedness of the operator L whenever PDEs are considered.It is found that a set of additional order conditions must be satisfied to guarantee convergence order p under suitable assumptions; one requires the linear operator L to be the infinitesimal generator of an analytic semigroup, and that the nonlinear function satisfies a Lipschitz condition.The authors are also able to give an example where order reduction is seen numerically for schemes not satisfying the additional conditions.But the conditions are rather restrictive, and in [13] exponential integrators of (nonstiff) order four are tried out numerically for a number of well-known semilinear PDEs, and no order reduction is seen, despite the fact that these integrators do not satisfy all the required conditions for order four as given in [12].This shows that the issue of determining the order behavior of exponential integrators for PDEs is indeed a subtle one, and remains today in an unsatisfactory state of nonresolution.
2. B-series and order conditions.Repeated differentiation of (1.1) with respect to time yields etc.The exact solution of (1.1) has a formal expansion where each term in the qth derivative corresponds in an obvious way to a rooted bicolored tree.Let for instance • ∼ F (•) = N (u) and • ∼ F (•) = Lu be the two trees with one node.Next, define B + as the operation which takes a finite set of trees {τ 1 , . . ., τ μ } and connects their roots to a new common black root.Similarly, τ = W + (τ ) connects the root of τ to a new white root resulting in the tree τ associated to F (τ ) = L • F (τ ).It suffices here to allow W + to act on a single tree and not on a set of trees.To each tree τ with q nodes formed this way, there exists precisely one term, F (τ ) called an elementary differential, in the qth derivative of the solution of (1.1).For q > 1 it is defined recursively as We may denote by T the set of all bicolored trees such that each white node has at most one child, and set T = T b ∪ T w the union of trees with black and white roots, respectively.Introducing the empty set ∅, and using the convention B + (∅) = •, W + (∅) = •, we may write The same bicolored trees used here also appear in the linearly implicit W -methods; see Steihaug and Wolfbrandt [21] as well as the text [10] by Hairer and Wanner.Following, for instance, the text by Hairer, Lubich, and Wanner [8], we may work with formal B-series.For an arbitrary map c : T ∪ ∅ → R, we let the formal series where the m i s count the number of equal trees among τ 1 , . . ., τ μ .
The further derivation of order conditions is based on the assumption that both the exact and numerical solution possess B-series of the form (2.4), say B(e, u 0 ) and B(u 1 , u 0 ), respectively.We refer to [1] for details, and present only the final result.
Theorem 2.1.Let T ⊂ T be the set of bicolored rooted trees such that every white node has precisely one child.An exponential integrator defined by (1.2)-(1.3)has order of consistency p if where The trees in T with at most four nodes are listed in Table 2.Note that even though all trees in the set T feature in the B-series for the exact and numerical solutions, it suffices to consider a subset T consisting of all trees in T except those with a terminal white node.There is an interesting connection between the set of trees T and the trees used to develop the order theory for composition methods in [19].White nodes appear as connected strings of nodes which, except for the root, have exactly one parent and one child, and always terminate in a black node.Therefore one can remove all white nodes and assign to the terminating black node the number of removed nodes plus one.Black nodes not connected to a white node are assigned the number one.These multilabelled trees are precisely those appearing in [19], they can also be identified as the set of rooted trees of nonempty sets.The generating function for these trees is well known, The number of order conditions for each order 1 to 9 is 1, 2, 5, 13, 37, 108, 332, 1042, and 3360.

Construction of exponential integrators.
The schemes of Lawson [15] are exponential integrators derived simply by introducing a change of variable, w(t) = e −tL u(t) in (1.1), and by applying a standard Runge-Kutta scheme to the resulting ODE.This approach results in a formula for w 1 in terms of w 0 .By setting u n = e tL w n one gets a scheme of the form (1.2)-(1.3) in which as noted by Lawson in [15].This scheme has order p if the underlying scheme determined by α j,0 r and β r,0 is of order p.This gives us a very useful tool for constructing exponential integrators with given underlying Runge-Kutta schemes.We express this in the following proposition.
Proposition 3.1.Suppose that the coefficients α j,0 r and β r,0 , 1 ≤ r, j ≤ s define a Runge-Kutta scheme of order p.Then, any exponential integrator of the form is of order p.In the above expression we use 0 0 := 1.
Proof.Order conditions for exponential integrators of order p involve α j,m r for 0 ≤ m ≤ p − 2 and β r,m for 0 ≤ m ≤ p − 1.On the other hand, the Lawson schemes must satisfy the order conditions for exponential integrators, and their values for these coefficients are precisely those specified in the proposition.
It is convenient to introduce finite dimensional function spaces V a and V b to which the respective coefficient functions a j r (z) and b r (z) will belong.For the purpose of calculations, it is also useful to work with basis functions ψ k for these spaces, There is a technical assumption that we will adopt to the end of this note.
Assumption 3.2.Any finite dimensional function space V of dimension K used for coefficient functions a j r (z) or b r (z) has the property that the map from is injective.Equivalently, any function in V is uniquely determined by its first K Taylor coefficients.

|τ |
Tree 3.1.Deriving schemes with natural continuous extensions.The approach of Krogstad in [14] is to approximate the nonlinear function N (u(t 0 + θh)), 0 < θ < 1 with a polynomial in θ.Assuming that the functions a j r (z) for the internal stages are given, one lets N (u(t n + θh)) be approximated by where N r = N (U r ) are the stage derivatives and w r (θ) are polynomials of degree d, with w(0) = 0, such that N (t 0 + θh) approximates N (u(t 0 + θh)) uniformly for 0 < θ < 1 to a given order.Replacing the exact problem with the approximate one, v = Lv + N (t), v(t 0 ) = u 0 one finds we then define the functions Thus, here the function space [5] presented a fourth order scheme using these basis functions with K b = 3. Krogstad [14] also derived a variant of their method by using a continuous extension as just explained.In [22] Zennaro developed a theory which generalizes the collocation polynomial idea to arbitrary Runge-Kutta schemes.The approach was called natural continuous extensions (NCE).By making a slight modification to the approach of Zennaro, one can find a useful way of deriving exponential integrators as well as providing them with a continuous extension.
Definition 3.3.We call N (t) of (3.4) a natural continuous n-extension (NCNE) of degree d of the exponential integrator where u(t) is the exact solution of (1.1) satisfying u(t 0 ) = u 0 ; 3.

7) for every smooth matrix-valued function G(t).
It is important to note that the polynomial N (t) only depends on the stages N r and the weights b r (0) = β r,0 corresponding to the underlying Runge-Kutta scheme.We also observe that since the w r (θ) does not depend on L, an NCN E as defined above is also an NCE in the sense of Zennaro for the system u = N (u).Before discussing the existence of NCN Es, we motivate their usefulness in designing exponential integrators.Suppose an underlying Runge-Kutta method has been chosen, and that an NCN E has been found.Then we can determine the functions b r (z) in order to obtain an exponential Runge-Kutta method of the same order as the underlying scheme.Proof.The exponential integrator we consider is obtained by replacing (1.1) by over the interval [t 0 , t 1 ] and by solving (3.8) exactly.We subtract (3.8) from (1.1) to obtain We may solve this equation to obtain the last equality is thanks to (3.7).
A reinterpretation of a result by Zennaro [22] combined with Proposition 3.1 leads to the following theorem.
Theorem 3.5.Suppose that an underlying Runge-Kutta scheme with coefficients α j,0 r and β r,0 of order p is given.Then it is possible to find a set of coefficient functions a j r (z) with a j r (0) = α j,0 r such that an NCNE of degree d = p+1 where ν * is the number of distinct elements among c 1 , . . ., c s .Corollary 3.6.For every underlying Runge-Kutta scheme, there exists an exponential integrator whose coefficient functions b r (z) are in the linear span of the functions {φ 0 (z), . . ., φ d−1 (z)}, where d = p+1 2 .Note, in particular, that one can derive fourth order exponential integrators using linear combinations of just φ 0 (z) and φ 1 (z) for b r (z), which is one less than what Cox and Matthews used; we present a specific example in section 4.

3.2.
Lower bounds for K a and K b .We start establishing lower bounds for the number of necessary basis functions ψ k by proving an ancillary result.
Lemma 3.7.Let q ≥ 0 be an integer.The matrix T q ∈ R d×d with elements Proof.Let w = (w 1 , . . ., w d ) T ∈ R d be arbitrary, and consider the polynomial We compute So T q w = 0 is equivalent to p (j) (1) = 0 for 0 ≤ j ≤ d − 1.Since p(x) is of the form x q+d r(x) where r(x) is a polynomial of degree at most d − 1, it follows that p(x) ≡ 0 so that w = 0.As φ ! for φ k defined by (3.5), we get as an immediate consequence of this lemma that the function spaces V = span(φ q , . . ., φ q+K−1 ), q ≥ 0 satisfy Assumption 3.2.
Theorem 3.8.For an exponential integrator of order p, the dimension of the function spaces V a and V b are bounded from below as follows: Proof.We will show that using smaller values of K a or K b than dictated by (3.9) is incompatible with the order conditions for a scheme of order p.Let V a and V b be arbitrary function spaces, satisfying Assumption 3.2, let V denote either of them, and let d = dim V .If f ∈ V , then there are numbers w 0 , . . ., w d−1 such that Consider the bicolored trees τ m,k q defined by which consist of a string of q ≥ 0 black nodes followed by a string of m > 0 white nodes with a bushy tree of k + 1 black nodes grafted onto the topmost leaf of the white nodes.We shall use these trees with q = 0 for proving the bound on K b and with q = 1 for K a .The density of τ m,k q is given by The trees corresponding to order conditions for a scheme of order p have at most p nodes, |τ m,k q The conditions corresponding to τ da,k 1 can be expressed as The conditions for τ d b ,k In both cases (d = d a or d b ), we end up with a (d + 1) × d linear system of equations for determining w m , m = 0, . . ., d − 1.This system is of the form for q ∈ {0, 1} and is solvable only if the matrix with elements is singular.However, Lemma 3.7 implies that the matrix T q is invertible so the linear system is inconsistent.It is hence not possible to choose Some remarks regarding the implications of Theorem 3.8 are in order.First, note that the bounds in the theorem are not proved to be sharp; however, Theorem 3.5 ensures that the lower bound is attainable for the dimension of V b if a basis is given by the functions φ k of (3.5).However, this result does not apply to the space V a of the functions a j r (z).For instance, in the case p = 5, one can prove that it is indeed possible to take K a = 2, but V a cannot be the span of φ 0 and φ 1 .But an example of a feasible two-dimensional space is that with basis ψ 0 (z) = φ 1 (z) and ψ 1 (z) = φ 1 ( 3 5 z).A particular scheme is given in Table 4, though the usefulness of the bounds are questionable in this particular example.Using say V b = span{φ 0 , φ 1 , φ 2 } combined with the above choice of V a requires the computation with a total of 4 basis functions, whereas only 3 are necessary if one instead chooses V a = V b .
Furthermore, we note that the minimum attainable value of the parameters K a and K b depend only on the order p of the underlying Runge-Kutta scheme and the choice of the basis functions ψ k .Specifically, the coefficients of the underlying Runge-Kutta scheme do not influence the minimum values of K a and K b .

Examples of exponential integrators.
In this section we will present examples of exponential integrators.For fourth order methods, one will notice that some well-known schemes are obtained for particular choices of the free parameters, suggesting that a search on the entire space of parameters may result in schemes which in some sense may have better properties than the known methods.The scheme of order 5 presented at the end is only included as an illustration of the proposed procedure for solving the order conditions.It remains a subject of future research to establish to which extent higher order exponential integrators are useful for practical purposes.
The procedure we have used in constructing schemes may be summarized as follows: 1. Choose an underlying Runge-Kutta scheme.This determines α j,0 r and β r,0 .
Table 3 Coefficient function for a fourth order ETD scheme with classical RK4 as underlying scheme.Basis functions given by (3.5).
2. Choose basis functions ψ k (z) for the coefficient functions and determine K a and K b .3. Use the order conditions for the trees of the form W m + (τ C ), where τ C is a tree with only black nodes, and determine β r,m , for 1 ≤ m ≤ K b − 1; see also (3.2). 4. Identify order conditions which are linear in c r = s j=1 α j,1 r and which otherwise depend only on β j,m r , 0 ≤ m ≤ K b − 1 and α j,0 r , and solve for c r . 5. Identify remaining conditions which depend linearly on α j,1 r .Solve for α j,1 r together with c r = s j=1 α j,1 r .Repeat this procedure to solve for α j,m r , 2 ≤ m ≤ K a − 1. 6. β r,m are now uniquely determined for m ≥ K b and α j,m r for m ≥ K a by (3.10).Verify all remaining order conditions for β r,m , K b ≤ m ≤ p − 1 and for α j,m r , K a ≤ m ≤ p−2.If inconsistencies appear, the basis functions are not feasible.7. Verify all remaining order conditions.
In most cases we have considered, once α j,0 r and β r,0 have been chosen, one can find the remaining α j,m r independently of the β r,m .Most of the exponential integrators we find in the literature are based on the classical fourth order scheme of Kutta, and it is typical that one can combine a j r (z) from one scheme with b r (z) from another scheme and still get overall order four.
In the class of ETD schemes, proposed by Cox and Matthews in [5] and Krogstad in [14], the space V b is spanned by the three functions φ 0 , φ 1 , and φ 2 of (3.5).However, in the former reference, dim V a = 2 with a basis {φ 0 (z/2), zφ 0 (z/2) 2 }.This V a coincides with the one used in [3] given in Table 1(a).
Another choice is to use φ k (z) of (3.5) both for V a and V b .In Table 3 we characterise all resulting schemes with K a = 2 and K b = 3.It is interesting to note that Theorem 3.8 predicts K a ≥ 2 and K b ≥ 2, and indeed, by choosing γ 1 = 1 3 and γ 2 = − 1  3 , we see that φ 2 disappears from the b r (z)-functions.Choosing γ 1 = γ 2 = 0, we recover the b r (z)-functions obtained in [5].
Letting V b be spanned by ψ 0 (z) = φ 0 (z) and ψ 1 (z) = φ 0 (z/2), one obtains the   These weights coincide with the ones derived in the fourth order scheme in [3] given in Table 1.Yet another choice is to let V b consist of functions of the form p(z)φ 0 (z), where p(z) is a polynomial of degree 1, and we recover b r (z) as in Table 1(b).Finally, we give an example of a fifth order exponential integrator based on a scheme of Fehlberg.As indicated in section 3.2, we take dim V a = 2 with basis ψ 0 (z) = φ 1 (z) and ψ 1 (z) = φ 1 ( 3 5 z).For V b we use the basis ψ k (z) = φ k (z) for k = 0, 1, 2. The resulting coefficient functions are given in Table 4.
In summary, this paper presents a complete order theory for exponential integrators of the form (1.2)-(1.3).From deriving order conditions by means of bicolored trees to proving bounds for the lowest possible number of basis functions, the results presented herein provide a general framework for constructing schemes of this type.A number of issues are, however, not addressed in the present paper.These include systematically choosing basis functions ψ k , and how to construct schemes with low error constants.
Exponential integrators are interesting from the point of view of handling unbounded or stiff operators, yet the order theory does not say anything about what happens for large eigenmodes of L in (1.1).Determining conditions for favorable behavior in light of such operators should be an arena for future work.

Table 1
Examples of schemes in general format for exponential integrators.

Table 4
Coefficient functions for a fifth order exponential integrator with Fehlberg's fifth order RK as the underlying scheme.Here a j i