Splitting methods with complex coefficients

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 sixth-order integrators and analyze their main features on a pair of numerical examples, in particular how the errors are propagated along the evolution.


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.
It turns out that ψ h can also be written in terms of χ h and χ * h ψ h = ϕ [2] α 2s h • ϕ as long as Equivalently, [15]. A relevant consequence of this property is that, starting with the coefficients a i , b i of a given splitting method, we can get the coefficients α i for the composition (7), which can be then applied in a more general setting with the maps χ h and χ * h 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, 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, time-symmetry 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 r ≥ 3 [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 with k = s (respect. k = 2s). Obviously, this sum vanishes only if at least one of the α i is negative. In consequence, the flows ϕ h , j = 1, . . . , m − 1 in (4) and ϕ [j] h , j = 2, . . . , m − 1 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 k = 2s leads to In a similar way, using the same condition with α 0 = α 2s+1 = 0, one has and then at least one a i as well as one b i 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 u t = ∆u on the unit interval with homogeneous Dirichlet boundary conditions. Then the corresponding generated semigroup is well defined only for t ≥ 0 [11].
More generally, consider the nonlinear heat equation with functions v i real and positive, and D i ≡ ∂/∂x i , on a certain domain Ω ∈ R d . 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 r = 2 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 ( = 1): A technique used in practice to obtain the eigenvalues and eigenfunctions for a given potential V consists in numerically integrating the equation (after spatial discretization) along pure imaginary times (τ = −it). Equivalently, the equation to be analyzed is 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 time-dependent 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 exp(z∆) (z ∈ C) associated with the Laplacian is well defined (in a reasonable distributional sense) if and only if Re(z) ≥ 0 [7]. More generally, it is possible to extend the semigroup related with parabolic PDEs into a sector in the right half plane of C [11].
The aim of this paper is to review some of the splitting methods with complex coefficients published in the literature, propose new sixth-order 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.

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 S [2] . Thus, one gets a third-order method as 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 s = 2 [8,11]. A fourth-order integrator can be obtained with the symmetric composition αh .
Methods (15), (16) and (17) can be used to generate recursively higher order composition schemes as S βh .
Here the coefficients have to verify the conditions α + β = 1, α n+1 + β n+1 = 0, whence if n is odd, and β = 1 − α. The choice l = 0 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: with 2α + β = 1, 2α n+1 + β n+1 = 0.
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 α j 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 Re(α j ) > 0. Observe that any method of the form (10) with coefficients having positive real part can be expressed in terms of the elementary flows ϕ λh with Re(λ) > 0 when S [2] h is taken as the leapfrog (5). This is also true, of course, in the more general case when f is split in m parts and S [2] h is taken as the symmetric second order basic method (4). For instance, suppose that f in (1) is separable in two parts, so that S [2] h is given by (5). Then it is straightforward to check that the third order scheme (15) can be written as 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 Re(α j ) > 0 involving the minimum number s of compositions for a prescribed order. For instance, the minimal number of compositions for achieving order 6 is s = 7. The corresponding order conditions can be written as [6,10,18] where for each j = 1, . . . , s, This system of algebraic equations has several solutions with Re(α j ) > 0. Among them, we have chosen the two sets of coefficients collected in Table 1. The first one corresponds to a symmetric method, α s+1−i = α i , as scheme (16), and was already found by Chambers [8]. The second method is apparently new, and possesses the special symmetry α s+1−i = α * i , as scheme (15) (or (20) when expressed as (6)).

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, Table 1: Coefficients of two 7-stage sixth-order methods of type (10): S 7 6 is a symmetric method and S * 7 6 is conjugate to a symmetric method (symmetric in the real part of the coefficients and skew-symmetric in the imaginary part).
we take the Hamiltonian function H(q, p) = 1 2 (p 2 + q 2 ), with q, p ∈ R. The corresponding equations of motion are linear and can be written as so that the numerical solution at time t = h furnished by method (6) is given by hB e ashA e bshB · · · e b 2 hB e a 1 hA e b 1 hB x 0 .
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 (q, p) = (1, 1) and integrate the system (23) for t ∈ [0, 20000π] 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 q out = Re(q), p out = Re(p). Figure 1 shows the results obtained with the following methods: (i) S 2 3, the 2-stage third-order non-symmetric method (15), (ii) S 3 4, the 3-stage fourth-order symmetric method (16), (iii) S * 7 6, the 7-stage sixth-order non-symmetric method, (iv) S 7 6, the 7-stage sixth-order symmetric method. The coefficients of these two 6th-order methods are collected in Table 1. The time step is chosen such that all methods require 27-28 evaluations per period. Notice the significant difference in the qualitative behavior of the numerical solution. Whereas the error grows exponentially for the symmetric methods S 3 4 and S 7 6, this is not the case for S 2 3 and S * 7 6, 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 K(h) in (24) is given explicitly by
In this way, one gets where p(h), d(h) (respectively, q(h), e(h)) are even (resp. odd) polynomial functions having in general complex coefficients and det K(h) = p(h) 2 If the splitting method (6) is such that (as happens, in particular, when it comes from a composition of the form (10) with α s−j+1 = α * j ), then K(h) −1 = K(−h) * . More specifically, This implies that p(y), q(h), and e(h) are real polynomials, whereas the coefficients of d(y) are purely imaginary. Notice that this is precisely the case of methods S 2 3 and S * 7 6. If, on the other hand, the splitting method is symmetric, i.e., it is of the form (6) satisfying a s−j+1 = a j , b s−j+2 = b j (as happens, in particular, when it comes from a composition of the form (10) with α s−j+1 = α j ), then K(h) −1 = K(−h). This clearly implies that d(h) ≡ 0, but in general the polynomials p(h), q(h), and e(h) have complex coefficients. For instance, methods (16) (S 3 4) and S 7 6 are such that p(h) has non-real coefficients.
When a splitting method with matrix K(h) is used to integrate the harmonic oscillator, it is essential that p(h) ∈ R. Otherwise K(h) n grows exponentially with the number n of steps. As a matter of fact, the eigenvalues of K(h) are λ 1 = e iφ(h) and λ 2 = e −iφ(h) , where φ(h) = arccos(p(h)), and thus max(|λ 1 |, |λ 2 |) > 1 if p(h) ∈ R (and also if p(h) ∈ R and |p(h)| > 1). That is precisely the situation with methods S 3 4 and S 7 6, 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 −1 ≤ p(h) ≤ 1. In fact, the numerical solution can still be (weakly) unstable when p(h) 2 = 1 with q(h)d(h)e(h) = 0 [5]. Furthermore, it is shown in [5] that, for stable numerical solutions (that is, with a suitable 2 × 2 matrix Q(h) (typically close to the identity matrix). In consequence, the numerical solution x n = (q n , p n ) is such thatx n := Q(h)x n corresponds to the exact solution at t n = nh of a harmonic oscillator with frequencyω = 1/hφ(h) ≈ 1. This feature explains why schemes S 2 3 and S * 7 6, when applied to the harmonic oscillator (23) with h = π/7 and h = π/2 respectively, exhibit a linear error growth in positions and a bounded error in energy, since for such methods, p(h) = 1 − h 2 /2 + · · · is real and satisfies p(h) ∈ (−1, 1) for the values of h considered in the numerical experiments.

Example 2: The Volterra-Lotka problem
Consider now the Volterra-Lotka probleṁ 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 (u * , v * ) = (1, 2) 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 0 < u, v determined by the first integral I(u, v) = ln(uv 2 ) − (u + v).
The vector field f (u, v) = (u(v−2), v(1−u)) can be separated in two solvable parts and this can be done in different ways. We consider the following split: , 0) and f B = (0, v(1 − u)) (although the linear and nonlinear separation can also be considered).
We take as initial conditions (u 0 , v 0 ) = (2, 4), integrate up to t = 20000 × 2π and measure the relative error in the first integral, |I − I 0 |/|I 0 |. As in the previous example, we integrate using complex arithmetic and take the real parts of u and v only for representing the output. Figures 2-(a) and (b) show the results obtained for time steps h = 4mπ 210 and four times smaller h = mπ 210 , with m 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 so-called pseudosymplectic methods (integrators of order n which preserve symplecticity up to order p > n), where the dominant errors behave as E I = Ch n + tDh p for some constants C and D. If p > n 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 u and v 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.

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 non-real 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 sixth-order 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.