High order integrators obtained by linear combinations of symmetric-conjugate compositions

A new family of methods involving complex coefficients for the numerical integration of differential equations is presented and analyzed. They are constructed as linear combinations of symmetric-conjugate compositions obtained from a basic time-symmetric integrator of order 2n (n $\ge$ 1). The new integrators are of order 2(n + k), k = 1, 2, ..., and preserve time-symmetry up to order 4n + 3 when applied to differential equations with real vector fields. If in addition the system is Hamiltonian and the basic scheme is symplectic, then they also preserve symplecticity up to order 4n + 3. We show that these integrators are well suited for a parallel implementation, thus improving their efficiency. Methods up to order 10 based on a 4th-order integrator are built and tested in comparison with other standard procedures to increase the order of a basic scheme.


Introduction
Composition methods constitute a standard tool to construct high-order numerical integrators for the initial value problemẋ = f (x), in particular when the vector field f possesses some qualitative property whose preservation by numerical approximations is deemed relevant [6,16]. Let S [2n] h denote a 2n-th order method, so that S [2n] h (x 0 ) = ϕ h (x 0 ) + O(h 2n+1 ), where x(h) = ϕ h (x 0 ) is the exact solution of Eq. (1.1) for a time step h. Then, if the coefficients α 1 , α 2 , . . . , α s satisfy some algebraic conditions, the composition of the basic scheme with step sizes α 1 h, α 2 h, . . . , α s h, i.e., is a new method of higher order 2n + m [9]. If in particular f is Hamiltonian and S [2n] h is symplectic, then the composition method (1.2) is also symplectic [16]. In general, any geometric property the basic method has in common with the exact solution is still shared by the higher-order scheme (1.2) if this property is preserved by composition [18]. Moreover, suppose S [2n] h is time-symmetric, namely, it satisfies where id is the identity map, for any h. Then, method (1.2) is also time-symmetric if the composition is left-right palindromic, i.e., α s+1−j = α j , j = 1, 2, . . ..
A well known class of composition methods is obtained by applying the triple-jump procedure [21,24]: with is a new method of order 2n + 2. The same technique can be applied again to S [2n+2] h , so that one can construct recursively time-symmetric methods of any order 2n + 2k, k = 1, 2, . . ..
When constructing high-order composition methods, real coefficients α 1 , . . . , α s are not the only option, however. In fact, the unavoidable existence of negative α j in (1.2) when the order is higher than two [5,15,20,22] typically imposes stability restrictions on the step size. This occurs in particular when Eq. (1.1) is the outcome of a parabolic differential equation discretized in space. In that case, considering complex coefficients with positive real part is also a valid alternative [12,17]. Even for problems where the presence of some α j < 0 is not particularly troublesome, composition methods with complex coefficients have also been proposed and analyzed from the preservation of properties viewpoint [10,7,13].
In the particular case of the triple-jump composition (1.3), in addition to the real solution (1.4), the complex one with the smallest phase is 5) and the resulting method has in fact smaller truncation errors than its real counterpart (1.4). If the basic scheme is time-symmetric and of order 2, then time-symmetric methods up to order 14 with coefficients having positive real part are possible by applying this technique [8].
The order can be raised by one instead with the simplest composition [3,22] ψ The choice = 0 gives the solution with the smallest phase, which we denote by γ [2n] : When the vector field f in (1.1) is real, then is complex, and so it is quite natural to project x 1 on the real axis and proceed to the next step only with (x 1 ). This is equivalent of course to integrating with the scheme R is. Nevertheless, it has been shown in [11] that R (1) h is pseudo-symmetric of order 4n + 3, in the sense that h is also pseudo-symplectic of order 4n + 3. In other words, projecting ψ [2n+1] h at each integration step leads to a numerical method that preserves geometric properties of the exact solution up to an order that is much higher than the order of the method itself. Pseudo-symplectic integrators have been previously considered in the literature, both in the context of Runge-Kutta [2] and polynomial extrapolation methods [4,14].
Moreover, as shown in [11], R h can be taken as the basis of the recursion producing methods of order 2(n + k), also pseudo-symmetric of order 4n + 3. Here the coefficients γ [2k] are given by Eq. (1.7). For future reference, we call (1.9) R-methods. Scheme (1.6) is a particular example of a symmetric-conjugate composition. These are composition methods of the form i.e., compositions (1.2) with α j ∈ C and Methods of this class, as shown in [7], possess remarkable preservation properties when considering its real part, In particular, if one takes a time-symmetric 2nd-order scheme as the basic method and the coefficients α 1 , α 2 , . . . are chosen in such a way that ψ h is of order 2n − 1, then (ψ h ) is of order 2n and pseudosymmetric of order 4n − 1 when the vector field f in (1.1) is real. If in addition f is a (real) Hamiltonian vector field and S [2] h is a symplectic integrator, then (ψ h ) is pseudo-symplectic of order 4n − 1. Since taking the real part of a symmetric-conjugate method is just a very special linear combination, it is quite natural to ask what happens when one considers a more general linear combination of symmetricconjugate compositions and their complex-conjugate, ψ h : is it possible to construct new methods of higher order whereas still preserving time-symmetry (and symplecticity) up to the order prescribed by the composition ψ (j) h ? If yes, how the new methods are built? Addressing these questions is precisely the subject of the present paper. In doing so, we present a new family of schemes of increasingly higher order well adapted for implementation in a parallel environment, requiring less computational effort than the Rmethods (1.9) but with the same qualitative properties.

Construction of the family of T -methods
In this section we construct the new family of integrators T (k) h and show explicitly that they are of order 2(n + k) and pseudo-symmetric of order 4n + 3 for k = 1, 2, 3. The same procedure can be formally extended to any k > 3. The analysis is based on the Lie formalism applied to the series of differential operators associated to the integrators.

Series of differential operators
As is well known, given a time-symmetric integrator S [2n] h of order 2n ≥ 2 one can associate a series of linear operators exp(Y (h)) so that for all functions g [9], with Here Y k are certain operators depending on the particular method and, for consistency, Y 1 = F , where F is the Lie derivative associated with f : The composition (1.2) then has the associated series which can be formally written as Ψ(h) = exp(V (h)) by repeated application of the Baker-Campbell-Hausdorff formula, with Here V 2n+1 , V 2n+2 , . . . are linear combinations of Lie brackets involving the operators Y 1 , Y 2n+1 , Y 2n+3 , . . . [18]. In the particular case of a symmetric-conjugate composition (1.10), terms V 2k in V (h) of even powers in h are pure imaginary, whereas terms V 2k+1 are real [7]. For a consistent symmetric-conjugate composition (1.10), i.e., verifying we get explicitly where µ n,k , σ n,k are homogeneous real polynomials of degree n in the coefficients α l , l = 1, . . . , s, and E n,k are elements Y j and independent Lie brackets involving these operators. In particular

Linear combinations of symmetric-conjugate compositions
Let us now consider the linear combination h is a consistent symmetric-conjugate composition of the form (1.10) with different coefficients α as the associated series of operators, where each V j (h) is of the form (2.4). Now, by following the same approach as in [11], we express Φ(h) as (2.7) Here This is done by applying the symmetric Baker-Campbell-Hausdorff formula to each product , a straightforward calculation shows that Therefore, and Φ(h) can also be written as In consequence, each term in φ h is time-symmetric up to terms h 4n+3 , with independence of the polynomials µ On the other hand, one has so that it is also true that

Order conditions
It is thus possible to obtain the order conditions for the method φ h in (2.5) by analyzing just the exponent of the central term in (2.8). From (2.7) it follows that In consequence, for consistent compositions ψ (j) h , j = 1, . . . , k, the conditions to be satisfied so that φ h is a method of order r are the following: • r = 2n + 2: c 2n+1,1 = 0

New schemes
Once identified the relevant order conditions, our next goal is to solve these equations with the minimum number of basic schemes in the compositions ψ Order r = 2n + 2. One needs to solve two equations to get a method φ h of order 2n + 2: consistency and c 2n+1,1 = 0. These can be satisfied by taking k = 1 and the simplest composition In other words, we recover the composition (1.6) and the R-method (1.8). Our first T -method (1.11) is thus Order r = 2n+4. Now we have to solve 3 order conditions in addition to consistency for the compositions ψ (j) h involved. As before, one could take in principle k = 1. In that case, the minimum number of basic maps in ψ (1) h is 4, just to have enough parameters to satisfy the order conditions. It turns out, however, that there are no solutions with the required symmetry α 4 =ᾱ 1 , α 3 =ᾱ 2 . In fact, if we take 2n+3,1 = 0, but µ 2n+3,2 = 0. On the other hand, if we take with the same values of α 1 , α 2 as before, then µ 2n+3,2 = −µ (1) 2n+3,2 , whereas still verifying that µ 2n+1,1 = µ (2) 2n+3,1 = 0. In consequence, by combining both compositions, one gets a method of order 2n + 4 and pseudo-symmetric of order 4n + 3. This corresponds to our second T -method, which reads explicitly Order r = 2n + 6. A total of 7 equations (including consistency) have to be solved in this case, so that we take a symmetric-conjugate composition involving s = 8 basic maps, With the choice it turns out that conditions c 2n+1 = c 2n+3,1 = c 2n+5,1 = 0 are automatically satisfied. By following the same approach as before, we permute the position of the coefficients and take the composition Then, one has µ 2n+3,2 = −µ (1) 2n+3,2 , so that ψ h leads to a method of order 2n + 4. More composition have to be incorporated, however, in order to verify conditions c 2n+5,2 = 0 and c 2n+5,3 = 0. The former is accomplished by both sums ψ but the later is satisfied only by adding up the four compositions. In summary, the linear combination leads to a method of order 2n + 6, denoted as T h . More explicitly, (2.11) The same procedure can be carried out in general, although more order conditions (and consequently more compositions involving more basic maps) have to be dealt with. This class of methods can be represented in a convenient way as follows. If we introduce the matrix of coefficients in the sense that each file of Γ 2n+2 corresponds to a particular symmetric-conjugate composition entering into the formulation of T (2) h . We can write analogously and moreover T h ; Γ 2n+4 ⊗ (Γ 2n+2 ⊗ Γ 2n ). In general, the coefficients in the T -method of order r = 2n + 2k are distributed according with the pattern

Numerical examples
We illustrate next the behavior of some of the previously constructed T -methods on a pair of numerical examples. The first one (the 2-dimensional Kepler problem) allows one to check preservation properties, whereas the second (a simple diffusion equation) is used as a test of their relative performance. In all cases we take as basic scheme S  previously considered in [8]. This integrator is intended for Eq. (1.1) when f can be decomposed as f (x) = f a (x) + f b (x) in such a way that each sub-probleṁ The implementation of all the integrators has been done in Python 3.7 running on Debian GNU/Linux 10 and the operations with complex arithmetics have been coded using the complex class of the numpy library.
Kepler problem. The Hamiltonian function for the planar two-body problem reads Here q = (q 1 , q 2 ), p = (p 1 , p 2 ), r = q , µ = GM , G is the gravitational constant and M is the sum of the masses of the two bodies. The corresponding equations of motion are theṅ Taking µ = 1 and initial conditions the resulting trajectory is an ellipse of eccentricity 0 ≤ e < 1. In this case ϕ h (respectively, ϕ h ) corresponds to the exact solution obtained by integrating the kinetic energy T (p) (resp., potential energy V (q)) in (3.3).
We take e = 0.6, integrate until the final time t f = 20π with the basic splitting method S [4] h given by (3.1) and schemes T (k) h , with k = 1, 2, 3 for several time steps and then we compute the average error in energy along the integration interval. Figure 1 (left) shows this error as a function of the number of evaluations of the basic scheme S [4] h . The diagram clearly exhibits the order of convergence of each method: order 4 for S [4] h , and orders 6, 8 and 10 for T  In the right panel we show the long-time behavior of the error in energy for each method when the step size is chosen so that all of them involve the same computational cost. We see that the error in energy is almost constant for t ≤ 2000π, as is the case for symplectic integrators. In other words, the lack of symplecticity at order h 12 has no effect in this integration interval. In addition, the scheme T h (red). The blue line corresponds to S [4] h .
We take N = 128 and integrate until t f = 1, where we compute the relative error U −U ex / U ex with each method T (k) h , k = 1, 2, 3, in addition to the basic scheme (3.1). The 'exact' solution U ex is taken as the output of the 8th-order composition method P8S15 of [8]. The corresponding efficiency diagram is shown in Figure 2, where the same notation is used for the curves depicted. Here also the higher degree integrators provide the best efficiency.

T -methods and R-methods
Methods T   h involves the sum of 2 2 k −2 compositions of 2 k appropriately weighted basic schemes: h . These numbers are collected in Table 1, when schemes R  schemes with respect to the former ones. In this respect, one should take into account that both explicit formulations (1.11) and (4.1) are directly amenable to parallelization, whereas this is less obvious for the recursion (1.9).
If one has a computer with, say, 2 threads, it is easy to estimate the effective number of evaluations of S • if > k − 1 then the number of evaluations is 2 k · 2 k−1− .
In Table 2 we collect these numbers for the first values of k in the particular case of 2 2 = 4 and 2 5 = 32 threads. We see that the implementation of the explicit expression of the R-methods is more advantageous than the recursive procedure already with a relatively small number of threads, and that, in any case, Tmethods require less computational effort.  To better illustrate this issue, we next compare the efficiency of the different methods when implemented on a computer able to execute 4 threads without loss of performance. The corresponding results are displayed h is the same, i.e., 2 k in both cases, the latter turn out to be more efficient. This is clearly visible in Figure 4, obtained again by applying the previous schemes to the Kepler problem (left) and the linear parabolic equation (right).
Finally, it is also illustrative to compare the efficiency of the new T -methods with the standard triplejump procedure, Eqs. (1.3)-(1.4), both applied to the same basic scheme (3.1). Thus, in Figure 5 we depict the results achieved by projecting S [6] h , S [8] h , and S [10] h at each step, together with T (k) h , k = 1, 2, 3 for the Kepler problem with the same parameters and final time t f = 20π. Here the effective number of evaluations of the basic scheme has been taken as 2 k for T -methods and 3 k for triple-jump. Not surprisingly, the new schemes turn out to be much more efficient.

Concluding remarks
The standard triple-jump procedure is a popular technique that allows one to construct numerical integrators for differential equations of arbitrarily high order by composition of a basic integrator of low order. It has nevertheless certain limitations: the number of basic maps grows rapidly with the order, and the main error terms are quite large in comparison with other specially built integrators. Moreover, they involve some negative coefficients when the order r ≥ 3, so that the resulting schemes cannot be used in particular when the initial value problem (1.1) results from the space discretization of a parabolic partial differential equation involving the Laplace operator. In this context it is quite natural to explore whether it is still possible using the triple-jump technique (1.3), but with the complex coefficients furnished by (1.5) as long as their real part is positive. It has been established that this is indeed the case, although once again they require an exceedingly large number of basic methods. For this reason, other alternatives for constructing high-order composition methods have also been proposed [8,12,17]. Among them, the class of schemes (1.8) possess some special features: starting from a time-symmetric basic scheme S [2n] h of order 2n, it is possible to construct recursively methods of order 2n + 2k, k = 1, 2, . . . that are still time-symmetric up to order 4n + 3. Moreover, if the differential equation in (1.1) has some qualitative properties (such as symplecticity or volume preservation) then these properties are still shared by the numerical solution up to order 4n + 3 [11].
Methods (1.9) are based on the simple symmetric-conjugate composition (1.6). As shown in [7], symmetric-conjugate composition methods still possess remarkable preservation properties when projected on the real axis at each integration step, and so it makes sense to consider more general linear combinations of methods within this class. The corresponding analysis has been carried out here, and as a result we have built a new class of schemes that essentially have the same preservation properties as methods (1.9), but requiring a much reduced computational cost. In addition, these methods are particularly well suited for their parallel implementation. The examples included show a significant improvement in efficiency with respect to schemes (1.9) and those obtained by applying the triple-jump procedure. h for T -methods in comparison with schemes obtained by triple-jump for the Kepler problem.