Splitting methods with complex coefficients
Abstract
Splitting methods for the numerical integration of differential equations of order greater than two involve necessarily negative coefficients. This order barrier can be overcome by considering complex coefficients with positive real part. In this work we review the composition technique used to construct methods of this class, propose new sixthorder integrators and analyze their main features on a pair of numerical examples, in particular how the errors are propagated along the evolution.

Instituto de Matemática Multidisciplinar, Universidad Politécnica de Valencia, E46022 Valencia, Spain.

Institut de Matemàtiques i Aplicacions de Castelló and Departament de Matemàtiques, Universitat Jaume I, E12071 Castellón, Spain.

Konputazio Zientziak eta A.A. saila, Informatika Fakultatea, EHU/UPV, Donostia/San Sebastián, Spain.
1 Introduction
Splitting methods for the numerical integration of differential equations constitute an appropriate choice when the associated vector field can be decomposed into several pieces and each of them is explicitly integrable.
Given the initial value problem
(1) 
with and solution , let us suppose that can be expressed as for certain functions , in such a way that the equations
(2) 
can be integrated exactly, with solutions at , the time step. Splitting methods intend to approximate the exact flow by a composition of flows . For instance,
(3) 
both provide firstorder approximations to the exact solution, since and similarly for (which is called the adjoint of and verifies ).
It is possible to get higher order approximations by introducing more maps with additional real coefficients, , in (3). Perhaps the most popular splitting method is the second order symmetric composition
(4) 
When in (1) is separable in two parts the above particularizes to
(5) 
and is known as the Strang splitting [22], the leapfrog or the Störmer–Verlet method [26], depending on the context where it is used. More generally, one may choose the coefficients , to achieve order with the composition
(6) 
It turns out that can also be written in terms of and
(7)  
as long as
(8) 
Equivalently,
(9) 
with . This relation remains valid if [15]. A relevant consequence of this property is that, starting with the coefficients of a given splitting method, we can get the coefficients for the composition (7), which can be then applied in a more general setting with the maps and of (3). A particular case widely used in practice to achieve high order approximations consists in considering compositions using the Strang splitting (4) as basic method,
(10) 
Splitting methods are, in general, explicit, easy to implement and preserve structural properties of the exact solution, thus conferring to the numerical scheme a qualitative superiority with respect to other standard integrators, especially when long time intervals are considered (see [6] for a review). Examples of these structural features are symplecticity, volume preservation, timesymmetry and conservation of first integrals. In this sense, splitting methods constitute an important class of geometric numerical integrators [10, 12, 14, 16, 17, 20].
It has been shown that some of the coefficients in splitting schemes (6) are negative when the order [9, 21, 24]. In other words, the methods always involve stepping backwards in time. An elementary proof of this feature can be worked out as follows. It is quite straightforward to check that one of the necessary condition for the composition (7) (respectively, (10)) to have order is
(11) 
with (respect. ). Obviously, this sum vanishes only if at least one of the is negative. In consequence, the flows in (4) and in (3) evolve with at least one negative fractional time step. On the other hand, by taking into account the link (8) among coefficients of (6) and (7), condition (11) with leads to
In a similar way, using the same condition with , one has
and then at least one as well as one are negative. It must be stressed that condition (11) still persists when the processing technique is used, so that the same conclusion also follows in this case [4].
In summary, the presence of negative coefficients in splitting methods of order higher than two is unavoidable if one restricts oneself to real coefficients. Of course, this does not suppose any special impediment when the flow of the ODE evolves in a group (such as in the Hamiltonian case), but may be unacceptable when the differential equation is defined in a semigroup [16], as occurs, for instance, with the simple heat equation on the unit interval with homogeneous Dirichlet boundary conditions. Then the corresponding generated semigroup is well defined only for [11].
More generally, consider the nonlinear heat equation
(12) 
with functions real and positive, and , on a certain domain . If a space discretization is carried out (either by finite differences or by a pseudospectral method), a large system of ODEs results which has to be numerically integrated in time. To this end, we can split the resulting equation into linear and nonlinear parts, but schemes of the form (7) or (10) of at most order can only be applied, since the resulting discrete Laplacian with negative fractional time steps is not well conditioned.
A closely related problem is the linear Schrödinger equation ():
(13) 
A technique used in practice to obtain the eigenvalues and eigenfunctions for a given potential consists in numerically integrating the equation (after spatial discretization) along pure imaginary times (). Equivalently, the equation to be analyzed is
(14) 
which can be considered as a linear heat equation. The system evolves to the ground state whose norm decreases exponentially in proportion to the value of its energy (eigenvalue). By orthogonalization, one can make the system to evolve to any other eigenfunction [1, 13]. In any case, whereas there is no special difficulties with numerically integrating equation (13) using a splitting methods with negative fractional time steps, this is not the case for (12) and (14) due to the presence of the Laplacian.
It has been noticed, however, that higher order splitting methods with complex coefficients having positive real part do exist [3, 16, 23, 24, 25]. These schemes were reported mainly for theoretical purposes but received very little attention as practical numerical tools. Perhaps the main reason was that working with complex arithmetic makes the schemes more involved and, in many cases, also considerably more costly from a computational point of view (usually, four times more expensive).
It is only within recent years that a systematic search for new methods with complex coefficients has been carried out and the resulting schemes have been tested in different settings: Hamiltonian systems in celestial mechanics [8], the timedependent Schrödinger equation in quantum mechanics [2, 19] and also in the more abstract setting of evolution equations with unbounded operators generating analytic semigroups [7, 11]. In this sense, we recall that the propagator () associated with the Laplacian is well defined (in a reasonable distributional sense) if and only if [7]. More generally, it is possible to extend the semigroup related with parabolic PDEs into a sector in the right half plane of [11].
The aim of this paper is to review some of the splitting methods with complex coefficients published in the literature, propose new sixthorder schemes in the class (10) and analyze them on a pair of simple numerical examples, to get a glance of the performance and main features of this kind of integrators and some of the difficulties involved.
2 Integrators with complex coefficients
Most of the existing splitting methods with complex coefficients have been constructed by applying the composition technique to the symmetric secondorder leapfrog scheme . Thus, one gets a thirdorder method as
(15) 
where the coefficients have to satisfy (11) together with the consistency condition
Due to its simplicity, this scheme has been rediscovered several times, either as the composition (15) [3, 23, 7] or by solving the order conditions required by (6) with [8, 11].
A fourthorder integrator can be obtained with the symmetric composition
(16) 
Although the necessary order conditions are the same, the timesymmetry of the composition rises the order by one (all the error terms at odd orders vanish identically):
with . Notice that for it is true that .
Another fourthorder method can be obtained by symmetrizing the thirdorder scheme (15), i.e.,
(17) 
Methods (15), (16) and (17) can be used to generate recursively higher order composition schemes as
(18) 
Here the coefficients have to verify the conditions , , whence
and . The choice gives the solutions with the smallest phase and allows one to build methods up to order six with coefficients having positive real part. This feature was stated in [25] and rediscovered in [11].
In a similar way, one may use recursively a symmetric three term composition, which allows to increase the order by two at each iteration:
(19) 
with
The solutions providing coefficients with the smallest phase are
and methods up to order eight with coefficients having positive real part are possible. Moreover, methods up to order fourteen of the more general form (10) with coefficients with positive real part are attainable [7, 11]. An interesting (and open) question is to determine whether arbitrarily high orders can be attained or wether, as for the previous compositions, there is an order barrier for methods of the form (10) with . Observe that any method of the form (10) with coefficients having positive real part can be expressed in terms of the elementary flows with when is taken as the leapfrog (5). This is also true, of course, in the more general case when is split in parts and is taken as the symmetric second order basic method (4).
For instance, suppose that in (1) is separable in two parts, so that is given by (5). Then it is straightforward to check that the third order scheme (15) can be written as
(20) 
with . This particular symmetry of the coefficients results in a method whose leading error terms at order 4 are all strictly imaginary [8].
Another question of practical nature is the construction of methods of the form (10) with involving the minimum number of compositions for a prescribed order. For instance, the minimal number of compositions for achieving order 6 is . The corresponding order conditions can be written as [6, 10, 18]
(21)  
(22) 
where for each ,
This system of algebraic equations has several solutions with . Among them, we have chosen the two sets of coefficients collected in Table 1. The first one corresponds to a symmetric method, , as scheme (16), and was already found by Chambers [8]. The second method is apparently new, and possesses the special symmetry , as scheme (15) (or (20) when expressed as (6)).

3 Numerical examples
3.1 Example 1: the harmonic oscillator
We consider the simple harmonic oscillator to illustrate some qualitative properties of the previous composition methods with complex coefficients. That is, we take the Hamiltonian function , with . The corresponding equations of motion are linear and can be written as
(23) 
so that the numerical solution at time furnished by method (6) is given by
(24) 
As is well known, for splitting methods with real coefficients the average error in energy remain constant for exponentially long times under suitable general conditions on the Hamiltonian. For the particular case of the harmonic oscillator and with a sufficiently small time step, this is true for all times, and the average error in positions grows only linearly.
We propose here to check whether this also holds for methods with complex coefficients. To do that, we take as initial conditions and integrate the system (23) for using a constant time step. We measure the error in position and energy of the output obtained by propagating the solution with the splitting method and then computing the real parts of the results , . Figure 1 shows the results obtained with the following methods: (i) , the 2stage thirdorder nonsymmetric method (15), (ii) , the 3stage fourthorder symmetric method (16), (iii) , the 7stage sixthorder nonsymmetric method, (iv) , the 7stage sixthorder symmetric method. The coefficients of these two 6thorder methods are collected in Table 1. The time step is chosen such that all methods require 2728 evaluations per period. Notice the significant difference in the qualitative behavior of the numerical solution. Whereas the error grows exponentially for the symmetric methods and , this is not the case for and , which show a performance analogous to standard splitting methods with real coefficients: bounded energy error and linear growth of error in positions. Of course, such a behavior deserves a theoretical explanation, which we pursue next.
The matrix in (24) is given explicitly by
In this way, one gets
where , (respectively, , ) are even (resp. odd) polynomial functions having in general complex coefficients and .
If the splitting method (6) is such that
(25) 
(as happens, in particular, when it comes from a composition of the form (10) with ), then . More specifically,
This implies that , , and are real polynomials, whereas the coefficients of are purely imaginary. Notice that this is precisely the case of methods and .
If, on the other hand, the splitting method is symmetric, i.e., it is of the form (6) satisfying
(as happens, in particular, when it comes from a composition of the form (10) with ), then . This clearly implies that , but in general the polynomials , , and have complex coefficients. For instance, methods (16) () and are such that has nonreal coefficients.
When a splitting method with matrix is used to integrate the harmonic oscillator, it is essential that . Otherwise grows exponentially with the number of steps. As a matter of fact, the eigenvalues of are and , where
and thus if (and also if and ). That is precisely the situation with methods and , and thus they are useless when integrating harmonic oscillators or systems that can be considered as close perturbations of harmonic oscillators with the partition (23).
From the previous comments, it is clear that instability will take place when integrating the harmonic oscillator unless . In fact, the numerical solution can still be (weakly) unstable when with [5]. Furthermore, it is shown in [5] that, for stable numerical solutions (that is, either or with ), one has
with a suitable matrix (typically close to the identity matrix). In consequence, the numerical solution is such that corresponds to the exact solution at of a harmonic oscillator with frequency . This feature explains why schemes and , when applied to the harmonic oscillator (23) with and respectively, exhibit a linear error growth in positions and a bounded error in energy, since for such methods, is real and satisfies for the values of considered in the numerical experiments.
3.2 Example 2: The Volterra–Lotka problem
Consider now the Volterra–Lotka problem
(26) 
This is a very simple nonlinear system which allows us to make a preliminary study about the behavior and performance of methods with complex coefficients in the transition process from a linear to a nonlinear problem. In a neighborhood of the steady state at the system can be considered close to a harmonic oscillator. The nonlinear contributions are manifest as we move away from it. The system evolves along periodic trajectories around the equilibrium point in the region determined by the first integral .
The vector field can be separated in two solvable parts and this can be done in different ways. We consider the following split: and (although the linear and nonlinear separation can also be considered).
We take as initial conditions , integrate up to and measure the relative error in the first integral, . As in the previous example, we integrate using complex arithmetic and take the real parts of and only for representing the output. Figures 2(a) and (b) show the results obtained for time steps and four times smaller , with the number of stages of each method. In this way, all methods require the same number of evaluations. Contrarily to the pure harmonic oscillator, we observe a secular error growth in the determination of the first integral for all methods which diminishes considerably when the time step is reduced. The observed behavior resembles what takes place with the socalled pseudosymplectic methods (integrators of order which preserve symplecticity up to order ), where the dominant errors behave as for some constants and . If the secular part of the error does not manifest for relatively long times when the time step is reduced.
We have repeated the same experiment but, after each time step, we discard the imaginary part of and and initiate the next step only with their real part. In other words, we project each component on the real axis at the end of each integration step. The results obtained are shown in Figures 2(c) and (d). Obviously, this way of proceeding does not preserve symplecticity any more but the results obtained suggest that a significant improvement in accuracy can be achieved.
4 Conclusions and outlook
We have presented a short review of the splitting and composition technique to build methods of order greater than two with complex coefficients with positive real part. This procedure allows to overcome the order barrier where splitting methods of order greater than two involve necessarily negative coefficients in the real space. In general, splitting methods with complex coefficients are considerably more expensive than the corresponding methods with real coefficients (about four times more expensive), and this make them hardly competitive in practice. For this reason, one can think that the main application of the new methods could be on parabolic PDEs, where higher order methods with real coefficients (which necessarily have some negative coefficientes) can not be used. However, there is a number of problems which evolve in the complex space where using methods which complex coefficients does not necessarily mean increasing the cost of the algorithm. This can be the case, for instance, of the Schödinger equation (13).
As for the practical implementation of splitting methods with complex coefficients, in [8] it is claimed that one has to carry the numerical integration in complex variables, and (for problems with real solutions) one should take either the real part of the variables or their modulus only for the output. However, we have observed that removing the imaginary part at each step, i.e. projecting on the real space at each step, the error grow can be considerably diminished in some cases. In the numerical examples considered in previous section, the linear error grow in the first integrals originate from different sources depending on wether the projection onto the real domain is performed after each step or not. In the first case, the projection after each step destroys symplecticity but only at a higher order, and the schemes can be considered as pseudosymplectic. In the second case, the method is actually symplectic and can thus be (formally) considered as an exact solution of a Hamiltonian system in the complex domain, which have qualitatively different properties to trajectories in the real domain. We have also noticed that the higher order methods present a considerably reduced error grow. Then, it seems appropriate to look for efficient higher order methods with complex coefficients. In general, symmetric splitting methods are desirable. However, we have shown that for the harmonic oscillator symmetric methods (with nonreal stability polynomial) present an exponential error grow, which is not the case for methods with the special symmetry (25). In a preliminary search of methods, we have presented a new sixthorder method with that special symmetry. This is an interesting subject to be further explored since many problems in different applications can be considered as perturbations to the harmonic oscillator.
Acknowledgements
This work has been supported by Ministerio de Ciencia e Innovación (Spain) under project MTM200761572 (cofinanced by the ERDF of the European Union). SB also acknowledges financial support from Generalitat Valenciana through project GV/2009/032.
References
 J. Auer, E. Krotscheck, and S.A. Chin. A fourthorder realspace algorithm for solving local Schrödinger equations. J. Chem Phys., 115:6841–6846, 2001.
 A.D. Bandrauk, E. Dehghanian, and H. Lu. Complex integration steps in decomposition of quantum exponential evolution operators. Chem. Phys. Lett., 419:346–350, 2006.
 A.D. Bandrauk and Hai Shen, Improved exponential split operator method for solving the timedependent Schrödinger equation. Chem. Phys. Lett., 176:428–432, 1991.
 S. Blanes and F. Casas. On the necessity of negative coefficients for operator splitting schemes of order higher than two. Appl. Numer. Math., 54:23–37, 2005.
 S. Blanes, F. Casas, and A. Murua. On the linear stability of splitting methods. Found. Comp. Math., 8:357–393, 2008.
 S. Blanes, F. Casas, and A. Murua. Splitting and composition methods in the numerical integration of differential equations. Bol. Soc. Esp. Math. Apl., 45:87–143, 2008.
 F. Castella, P. Chartier, S. Decombes, and G. Vilmart. Splitting methods with complex times for parabolic equations. BIT, 49:487–508, 2009.
 J.E. Chambers. Symplectic integrators with complex time steps. Astron. J., 126:1119–1126, 2003.
 D. Goldman and T.J. Kaper. thorder operator splitting schemes and nonreversible systems. SIAM J. Numer. Anal., 33:349–367, 1996.
 E. Hairer, Ch. Lubich, and G. Wanner. Geometric Numerical Integration. StructurePreserving Algorithms for Ordinary Differential Equations. SpringerVerlag, Second edition, 2006.
 E. Hansen and A. Ostermann. High order splitting methods for analytic semigroups exist. BIT, 49:527–542, 2009.
 A. Iserles, H. Z. MuntheKaas, S. P. Nørsett, and A. Zanna. Liegroup methods. Acta Numerica, 9:215–365, 2000.
 L. Lehtovaara, J. Toivanen, and J. Eloranta. Solution of timeindependent Schrödinger equation by the imaginary time propagation method. J. Comput. Phys., 221:148–157, 2007.
 B. Leimkuhler and S. Reich. Simulating Hamiltonian Dynamics. Cambridge University Press, 2004.
 R. I. McLachlan. On the numerical integration of ordinary differential equations by symmetric composition methods. SIAM J. Numer. Anal., 16:151–168, 1995.
 R.I. McLachlan and R. Quispel. Splitting methods. Acta Numerica, 11:341–434, 2002.
 R.I. McLachlan and R. Quispel. Geometric integrators for ODEs. J. Phys. A: Math. Gen., 39:5251–5285, 2006.
 A. Murua and J.M. SanzSerna. Order conditions for numerical integrators obtained by composing simpler integrators. Philos. Trans. Royal Soc. London, ser. A, 357:10791100, 1999.
 T. Prosen and I. Pizorn. High order nonunitary splitstep decomposition of unitary operators. J. Phys. A: Math. Gen., 39:5957–5964, 2006.
 J. M. SanzSerna and M. P. Calvo. Numerical Hamiltonian Problems. AMMC 7. Chapman & Hall, 1994.
 Q. Sheng. Solving linear partial differential equations by exponential splitting. IMA J. Numer. Anal., 9:199–212, 1989.
 G. Strang. On the construction and comparison of difference schemes. SIAM J. Numer. Anal., 5:506–517, 1968.
 M. Suzuki, Fractal decomposition of exponential operators with applications to manybody theories and Monte Carlo simulations. Phys. Lett. A, 146:319323, 1990.
 M. Suzuki. General theory of fractal path integrals with applications to manybody theories and statistical physics. J. Math. Phys., 32:400–407, 1991.
 M. Suzuki. Hybrid exponential product formulas for unbounded operators with possible applications to Monte Carlo simulations. Phys. Lett. A, 201:425–428, 1995.
 L. Verlet. Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard–Jones molecules. Phys. Rev., 159:98–103, 1967.