1 Introduction
Given a differential equation
(1) 
composition methods constitute a powerful technique to raise the order of a given integrator applied to (1) with timestep , as high as might be required, by considering expressions of the form
(2) 
where the coefficients are appropriately chosen so as to satisfy some universal algebraic conditions [HLW06, MSS99, CM09]. It is known in particular that if is of order , i.e. satisfies
where denotes the exact flow of (1), then will be at least of order (i.e., local error ) if the following two conditions are satisfied
(3) 
Given that these two equations have no real solution for odd
and arbitrary , a series of authors (e.g. [Suz90, Yos90]) suggested to start from a secondorder method and to consider symmetric compositions only, i.e., schemes with coefficients satisfying the additional conditionThis has led to socalled triplejump compositions (, ) obtained by iterating the process described above to construct a sequence of symmetric methods with even orders (see, e.g., [HLW06] pp. 44).
In spite of its simplicity, the triplejump rationale leads to inefficiencies for high orders as compared to methods obtained by solving directly the order conditions [HLW06]
. On top of this, it also suffers from the occurrence of negative timesteps, although this fact is not specific to triplejump methods and concerns all composition or splitting methods of orders higher than two. This, of course, is a severe limiting factor for equations where the vector field (usually an operator) is not reversible, the prototypical example of which being the heat equation. To circumvent this difficulty, several authors have suggested to use complex timesteps (or complex coefficients) in the context of parabolic equations
[CCDV09, HO09]. One indeed easily sees that, already for , solutions of equations exist in .Generally speaking, suppose that is an integrator of order , denoted in the sequel for clarity, and consider the composition (2) with ,
(4) 
Then, if the coefficients verify conditions , that is to say if
(5) 
then (4) results in a method of order , which can subsequently be used to generate recursively higher order composition schemes by applying the same procedure. The choice ,
(6) 
gives the solutions with the smallest phase. If we start with a symmetric method of order 2, , and apply composition (4) with corresponding coefficients (6), we can construct the following sequence of methods:
all of which have coefficients with positive real part [HO09]. The final method of order involves evaluations of the basic scheme . By contrast, there are composition methods of order (both with real and complex coefficients) involving just evaluations of [BCCM13, Yos90]. It is thus apparent that this direct approach does not lead to costefficient highorder schemes.
One should remark that the composition (4) does not provide a timesymmetric method, i.e., is not the identity map, even if happens to be symmetric. As we have mentioned before, symmetry allows to raise the order by two at each iteration by considering the triplejump composition
(7) 
starting from a symmetric method. Apart from the real solution, the complex one with the smallest phase is
(8) 
and symmetric methods up to order with coefficients having positive real part are possible if one starts with a symmetric secondorder scheme^{5}^{5}5It is actually possible to reach order 14 if, in the construction, one uses formula (7) alternatively with coefficients and coefficients [CCDV09].. These order barriers has been rigorously proved in [BCCM13].
The simple thirdorder scheme (4) corresponding to has been in fact rediscovered several times in the literature [BS91, CCDV09, Cha03, HO09, Suz92]. In particular, it was shown in [Cha03] that the method, when applied to the twobody Kepler problem, behaves indeed as a fourthorder integrator, the reason being attributed to the fact that the main error term in the asymptotic expansion is purely imaginary. In this note we elaborate further the analysis and provide a comprehensive study of the general composition (4), paying special attention to the qualitative properties the method shares with the continuous system (1). In addition, we show how it is possible to combine compositions and a trivial linear combination to raise the order, while still preserving the qualitative properties of the basic integrator up to an order higher than of the method itself.
2 Composition and pseudosymmetry or pseudosymplecticity
In what follows, we will assume for convenience that all values of in (1) lie in a compact set where the function is smooth. Before starting the analysis, it is worth recalling the notions of adjoint method and symplectic flow.
The adjoint method of a given method is the inverse map of the original integrator with reversed time step , i.e., . A symmetric method satisfies [Cha15, HLW06].
The vector field in (1) is Hamiltonian if there exists a function such that , where and is the basic canonical matrix. Then, the exact flow of (1) is a symplectic transformation, for [BC16, SSC94].
Definition 1
Let be a smooth and consistent integrator:

it is pseudosymmetric of pseudosymmetry order if for all sufficiently small , the following relation holds true:
(9) where the constant in the term depends on bounds of derivatives of on .

it is pseudosymplectic of pseudosymplecticity order if for all sufficiently small , the following relation holds true when is applied to a Hamiltonian system:
(10) where the constant in the term depends on bounds of derivatives of on .
Remark 1
A symmetric method is pseudosymmetric of any order , whereas a method of order is pseudosymmetric of order . A similar statement holds for symplectic methods.
As a first illustration of Definition 1, let us consider again a symmetric 2ndorder method and form the composition
with . Then, if the vector field under consideration is realvalued, its real part
is a method of order 4 and pseudosymmetric of pseudosymmetry order 7. This result is a consequence of the fact that
and the following general statement, which lies at the core of the construction procedure described in this paper.
Proposition 1
Let be any consistent smooth method for equation (1) and consider the new method
Assume also that is pseudosymmetric of order . Then is of pseudosymmetry order . If is furthermore of pseudosymplecticity order , then is of pseudosymplecticity order .
Proof:
By assumption, there exists a smooth function , defined for all in a compact set and for all sufficiently small real , such that
(11) 
so that
Composing the third relation of (11) from the left by , we obtain
(12) 
where the term depends on bounds of the derivatives of and on . Similarly, composing the second relation of (11) from the right by , we get
(13) 
As a consequence, we have
We are then in position to write
which proves the first statement. Now, if is in addition of pseudosymplecticity order , then its adjoint is also of pseudosymplecticity order , so that relation (11) leads to
which implies that
As an immediate consequence, we have that
which proves the second statement.
This result can be rendered more specific as follows:
Proposition 2
Let be a smooth method of order and pseudosymmetry order . Let us consider the composition method
(14) 
where the coefficients and satisfy both relations and . Then the method
(15) 
is of order
(16) 
when the vector field in (1) is real, and of pseudosymmetry order
(17) 
If in addition, is a (real) Hamiltonian vector field and is of pseudosymplecticity order , then is of pseudosymplecticity order
(18) 
Remark 2
Note that in Proposition 2, one has necessarily . This can be seen straightforwardly by a direct computation of with .
Proof:
Noticing that and are complex conjugate and (1) is real, and taking into account that is of pseudosymmetry order , we have
(19) 
Moreover, by construction, is at least of order , so that
(20) 
and altogether
(21) 
Now, since the pseudosymmetry order of is at least , the method
is, according to Proposition 1, of pseudosymmetry order and of pseudosymplecticity order . The first (16), second (17) and third (18) statements on orders then follow from (21).
In the Appendix we provide an alternative proof of Proposition 2 based on the Lie formalism, which allow us, in addition, to generalise the previous result on pseudosymplecticity to other geometric properties the continuous system may possess (such as in volume preserving flows, isospectral flows, differential equations evolving on Lie groups, etc.).
Notice that, according with Proposition 2, if we start from , that is to say from a basic symmetric () and symplectic () method of order , we get a method of order that is pseudosymmetric and pseudosymplectic of order just by considering the simple composition (14) and taking the real part of the output at each time step. If this technique is applied to a symmetric and symplectic method of order , i.e. with , then is of order and pseudosymmetric and pseudosymplectic of order .
Let us consider, in particular, the order symmetric scheme (7) with as basic scheme. Then, the resulting order integrator only requires the evaluation of secondorder methods , whereas the corresponding order scheme obtained by the triplejump technique involves evaluations. This number is reduced to by considering general compositions of [BCCM13]. If we take this order composition of schemes as the basic method , the resulting integrator of order , , involves the evaluation of , whereas evaluations are required by pure composition methods. Notice that is pseudosymmetric and pseudosymplectic of order 15, so that for values of sufficiently small, it preserves effectively the symmetry up to roundoff error while the drift in energy for Hamiltonian systems is hardly noticeable.
3 Families of pseudosymplectic methods
There is another possibility to increase the order, though, and it consists in applying the technique of Proposition 2 recursively. Thus, if denote by the method of eq. (15), we propose to apply the following recurrence:
(22)  
where is given by (6). Then, according with Proposition 2, it is possible to raise the order up to the pseudosymmetry order of the underlying basic method . Thus, in particular, the maximum order one can achieve by applying this technique to the basic symmetric method is 7, whereas if we start with a basic symmetric method of order 4, , the maximum order is . It is from a symmetric method of order and so on and so forth.
To give an assessment of the computational cost of the methods obtained by applying this type of composition, we notice that the computation of and required to form by (22) at the intermediate stages can be done in parallel, whereas at the final stage it only requires taking the real part. Thus, the method of order 6 constructed recursively from only requires the effective computation of 4 basic methods .
Starting from a symmetric secondorder method , say Strang splitting for instance, it is important to monitor the sign of the real part of all coefficients involved in the previous iteration. It is immediate to see that in the recursive construction
envisaged in the recurrence (22), the basic method is used with the following coefficients
Given the expression of (see (6)), these coefficients have arguments of the form
so that their maximum argument is
For all the coefficients to have positive real parts, a necessary and sufficient condition is thus that
It clearly holds for methods , and , of respective orders , and , since . Similarly, starting form a symmetric method of order having real or complex coefficients with maximum argument , the condition becomes
For instance, suppose that in (1) can be split as , so that the exact flows and corresponding to and , respectively, can be computed exactly. Then, the following composition
(23) 
with
provides a 4thorder symmetric scheme (see [CCDV09]). Taking (23) as basic method we get so that
and thus all methods
of respective orders , , , and obtained by the procedure (22) have all their coefficients with positive real parts. As far as the part is concerned, the maximum argument is less than .
4 Numerical experiments
In this section we illustrate the previous results on several numerical examples, comprising Hamiltonian systems and partial differential equations of evolution previously discretised in space.
4.1 Harmonic oscillator
We first consider the simple harmonic oscillator, with Hamiltonian
If we denote by the exact matrix evolution associated with the Hamiltonians , and , i.e., , then
respectively. We take as basic symmetric (and symplectic) scheme the leapfrog/Strang splitting:
(24) 
and compute the first three iterations in (22). In Table 1 we collect the main term in the truncation error for the resulting integrators , . We also check their timesymmetry and the preservation of the symplectic character of the approximate solution matrix by computing its determinant (a matrix is symplectic iff
). One can observe that these results are in agreement with the previous estimates.
Next we take initial conditions , integrate until the final time with , , and and compute the relative error in energy along the evolution. The result is depicted in Figure 1. We see that for and the error in energy is almost constant for a certain period of time, and then there is a secular growth proportional to .
4.2 Kepler problem
Next, we consider the twodimensional Kepler problem with Hamiltonian
(25) 
Here , , is the gravitational constant and is the sum of the masses of the two bodies. Taking and initial conditions
(26) 
if , then the solution is periodic with period , and the trajectory is an ellipse of eccentricity . Note that the gradient function must here be implemented carefully so as to be analytic for complex values of . Here, we define it using the following determination of the complex logarithm (analytic on the complex plane outside the negative real axis):
As a consequence, the analytic continuation of the function writes
where and .
Here, as with the harmonic oscillator, we take as basic method the 2ndorder Strang splitting
(27) 
where (respectively, ) corresponds to the exact solution obtained by integrating the kinetic energy (resp., potential energy ) in (25).
We take , integrate until the final time with Strang and the schemes obtained by the recursion (22) with for several time steps and compute the relative error in energy at the final time. Figure 2 show this error as a function of the inverse of the step size to illustrate the order of convergence: order 2 for Strang, order 4 for and order 6 for . For , and contrary to what happens to the harmonic oscillator, the observed numerical order is higher than expected, varying between 7 and 8. We do not have at present a theoretical explanation for this phenomenon. Figure 2 (right) depicts the time evolution of this error when the final time is .
4.3 The semilinear reactiondiffusion equation of Fisher
Our third testproblem is the scalar equation in onedimension
(28) 
with periodic boundary conditions on the interval . Here is a nonlinear reaction term. For the purpose of testing our methods, we take Fisher’s potential [Sar15]
as considered for example in [BCCM13].
The splitting corresponds here to solving, on the one hand, the linear equation with the Laplacian (as
), and on the other hand, the nonlinear ordinary differential equation
with initial condition , whose analytical solution is given by the welldefined (for small enough complex time ) formula
Here we aim to solve Eq. (28) with periodic boundary conditions on the interval , and initial condition . Numerically, the interval is discretised on a uniform grid, i.e., , and is approximated by Fourier pseudospectral methods. In this way we construct a vector with components , . If we denote by the whole numerical solution computed by a certain integrator with step size from until the final time, and by the corresponding numerical solution computed by the same integrator with step size , then the quantity is a good indicator of the convergence order.
Numerical simulations were carried out in quadruple precision (with Intel Fortran) such that roundoff errors are suppressed. Figure 3 shows the successive errors , at final time , of the methods obtained with the sequence (22) with the Strang splitting as the basic method (left) and the fourth order scheme given by (23) (right) with different time steps . One can clearly observe that the convergence order matches the previous analysis with a slightly better performance for the highest order, analogously to the Kepler problem. Figure 4 shows the successive errors versus the number of basic integrators in each case.
4.4 The semilinear complex Ginzburg–Landau equation
Our final test problem is the complex Ginzburg–Landau equation on the domain ,
(29) 
with , and initial condition . Here, , and denote real coefficients. In physics, the GinzburgLandau appears in the mathematical theory used to model superconductivity. For a broad introduction to the rich dynamics of this equation, we refer to [vS95]. Here, we will use the values , and , for which plane wave solutions establish themselves quickly after a transient phase (see [WMC05]). In addition, we set
so that the solution can be represented in Figure 5.
To apply the composition methods presented in previous sections, it seems natural to split equation (29) as
(30) 
whose solution is
Comments
There are no comments yet.