Composition Methods for Dynamical Systems Separable into Three Parts

New families of fourth-order composition methods for the numerical integration of initial value problems defined by ordinary differential equations are proposed. They are designed when the problem can be separated into three parts in such a way that each part is explicitly solvable. The methods are obtained by applying different optimization criteria and preserve geometric properties of the continuous problem by construction. Different numerical examples exhibit their improved performance with respect to previous splitting methods in the literature.


Introduction
Splitting methods are particularly useful for the numerical integration of ordinary differential equations (ODEs)̇ ≡ = ( ), when the vector field can be written as ( ) = ∑ =1 ( ), so that each subprobleṁ = ( ), ( 0 ) = 0 , = 1, … , is explicitly solvable, with solution ( ) = [ ] ( 0 ). Then, by composing the different flows with appropriate chosen weights it is possible to construct a numerical approximation to the exact solution (ℎ) for a time-step ℎ of arbitrary order [1]. Although splitting methods have a long history in numerical mathematics and have been applied, sometimes with different names, in many different contexts (partial differential equations, quantum statistical mechanics, chemical physics, molecular dynamics, etc. [2]), it is in the realm of Geometric Numerical Integration (GNI) where they play a key role, and in fact some of the most efficient geometric integrators are based on the related ideas of splitting and composition [3].
In GNI the goal is to construct numerical integrators in such a way that the approximations they furnish share one or several qualitative (often, geometric) properties with the exact solution of the differential equation [4]. In doing so, the integrator has not only an improved qualitative behavior, but also allows for a significantly more accurate long-time integration than it is the case with generalpurpose methods. In this sense, symplectic integration algorithms for Hamiltonian systems constitute a paradigmatic example of geometric integrators [5,6]. Splitting and composition methods are widely
Requirement (6) leads to a set of polynomial equations (the so-called order conditions), whose number and complexity grows enormously with the order. In particular, if = 1 (i.e., for a consistency method) one has The specific number of order conditions is determined in fact by the dimension of the homogeneous subspace of grade , 1 ≤ ≤ , of the free Lie algebra ( , , ) generated by the Lie derivatives , , corresponding to [ ] , [ ] , [ ] , respectively [1]. These dimensions are collected Table 1 for 1 ≤ ≤ 8.
where id is the identity map. Condition (7) is verified by left-right palindromic compositions, i.e., if in (5). Then all the conditions at even order are automatically satisfied. Thus, a symmetric method of order 4 requires solving 11 order conditions (instead of 32). Still, within this approach, one has to solve 56 polynomial equations to construct a method of order 6. Methods of this class have been systematically analyzed in [11]. In particular, it has been shown that if one aims to get schemes (5) of order 2 with the minimum number of maps, then the Strang splitting (4) is recovered. With respect to order 4, the following scheme was presented: In fact, 13 is the minimum number of maps required. More efficient schemes involving 17 and 25 maps can also be found in [11]. For simplicity, we denote method (8) as ( 1 1 1 2 2 3 2 3 2 2 1 1 1 ). More recently, in [12] a method involving 21 maps of the form has also been proposed and tested on several numerical examples.

Second Approach: Composition Methods
As it is clear from the previous considerations, constructing high order splitting methods for systems separable into three parts requires solving a large number of polynomial equations involving the coefficients, and this is a very challenging task in general. For this reason, we turn our attention to another strategy based on compositions of the first order with appropriately chosen weights. In other words, we look for integrators of the form verifying in addition the time-symmetry condition 2 +1− = for all .

Remark 1 Methods of the form
and +1− = (commonly referred in the literature as symmetric compositions of symmetric methods [1]) verify a much reduced number of order conditions and allows one to construct very efficient high-order schemes [3]. Notice that, since the Strang splitting (4) verifies  [2] ℎ = ℎ∕2 • * ℎ∕2 , then it is clear that when analyzing methods (10) we also recover schemes of the form (11).

Analysis in Terms of Exponentials of Operators
The analysis of the composition methods considered here can be conveniently done by considering the Lie operators associated with the vector fields involved and the graded free Lie algebra they generate.
It is worth remarking that the order conditions (15) are general for any composition method of the form (10), with independence of the particular basic first-order scheme ℎ considered, as long as ℎ and its adjoint * ℎ are included in the sequence. Thus, for instance, one might take the explicit Euler method as ℎ and the implicit Euler method as * ℎ , and also a symplectic semi-implicit method and its adjoint, leading to the symplectic partitioned Runge-Kutta schemes considered in [16].

Composition Methods of Order 4
Although one already gets a method of order 4 with only three stages, it is well known that the scheme (17) has large high-order error terms. A standard practice to construct more efficient integrators consists in adding more stages in the composition and determine the extra free parameters thus introduced according with some optimization criteria. Although assessing the quality of a given integration method applied to all initial value problem is by no means obvious (the dominant error terms are not necessarily the same for different problems), several strategies have been proposed along the years to fix these free parameters in the composition method (10). Thus, in particular, one looks for solutions such that the absolute value of the coefficients, i.e., is as small as possible, the logic being that higher order terms in the expansion (14) involve powers of these coefficients. In fact, methods with small values of 1 ( ) usually have large stability domains and small error terms [1]. In addition, for a number of problems, the dominant error term is precisely the coefficient 5 multiplying 5 in the expansion (14), so that it makes sense to minimize for a given composition to take also into account the computational effort measured as the number 2 of basic schemes considered. Here, as in [17], we construct symmetric methods with small values of 1 which, in addition, have also small values of 2 . For future reference, the corresponding values of the objective functions for the triple-jump (17) are 1 = 4.40483 and 2 = 4.55004, respectively.
Next we collect the most efficient schemes we have obtained with = 4, 5, 6 by applying this strategy.

= 4 stages. The composition is
and involves 17 maps when the basic scheme ℎ is given by (3). Now we have a free parameter, which we take as 1 . The minima of both 1 and 2 are achieved at approximately 1 = 0.358, and the resulting coefficients are collected in Table 2 as method XA 4 . In that case, 1 = 2.9084 and 2 = 3.1527.
= 5 stages. The resulting composition involves 21 maps when applied to a system separable into three parts. Minimum values for 1 and 2 are achieved when In consequence, the method can be written as ℎ with = 2 1 , = 2 5 . Then 1 = 2.3159 and 2 = 2.6111. This method, denoted XA 5 , was first proposed in [18] and analyzed in detail in [19].
= 6 stages. Analogously we have considered a composition involving three free parameters (and 25 maps when ℎ is given by (3)): The proposed solution is collected in Table 2 as method XA 6 leading to 1 = 2.0513, 2 = 2.4078. Notice how, by increasing the number of stages, it is possible to reduce the value of 1 and 2 as a measure of the efficiency of the schemes. This integrator has been tested in the numerical integration of the so-called reduced 1 + 1∕2 Vlasov-Maxwell system [20].
We could of course increase the number of stages. It turns out, however, that with = 7 one has the sufficient number of parameters to satisfy all the order conditions up to order 6, resulting in a method of the form (11) [15] involving 29 maps. More efficient 6th-order schemes can be obtained indeed by increasing the number of stages. Thus, in particular, with = 9 and = 11 one has the methods designed in [21] (37 maps) and [22] (45 maps), respectively, when the basic scheme is given by (3).

Third Approach: Splitting via Composition
We have already seen that there exists a close relationship between composition methods of the form (10) and splitting methods. This connection can be established more precisely as follows [24]. Let us assume that in the ODE (1) Since  (22) can be rewritten as the splitting scheme if 1 = 1 and (with 2 +1 = 0). Conversely, any integrator of the form (23) can be expressed in the form (10) with 0 = 0 for consistency. In consequence, any splitting method in principle designed for systems of the forṁ = ( ) + ( ) with no further restrictions on or can be formulated as a composition (10) which, in turn, can also be applied when is split into three (or more) pieces, . The performance will be in general different, since different optimization criteria are typically used. Notice that the situation is different, however, if splitting methods of Runge-Kutta-Nyström type are considered.
A particularly efficient 4th-order splitting scheme designed for problems separated into two parts has been presented in ( [23] Table 2) (method S 6 ) and will be used in our numerical tests. It is a time-symmetric partitioned Runge-Kutta method of the form (23), since the role played by and are interchangeable. When formulated as a composition method, it has six stages, i.e., it is of the form (21), with coefficients listed in Table 2. For comparison, the corresponding values of 1 and 2 are 1 = 2.4668 and 2 = 3.1648.

An Optimization Criterion Based on the Error in Energy
Very often, the class of problems to integrate are derived from a Hamiltonian function. In that case, Equation (1) is formulated aṡ In other words, { , } is the Poisson bracket of and . In this context, then, the Lie bracket of operators can be replaced by the real-valued Poisson bracket of functions [13].
It is well known that the flow corresponding to (25) is symplectic and in addition preserves the total energy of the system. If can be split as = + , then [ ] = , [ ] = and the splitting method (23) is also symplectic. Important as it is that the method shares this feature with the exact flow, one would like in addition that the energy be preserved as accurately as possible (since a numerical scheme cannot preserve both the symplectic form and the energy). A possible optimization criterion would be then to select the free parameters in such a way that the error in the energy (or more in general, in the conserved quantities of the continuous system) is as small as possible.
This criterion can be made more specific as follows [25]. First, we expand the modified Hamilto-niañ ℎ in the limit ℎ → 0 for a 4th-order splitting method (23). A straightforward calculation shows that̃ ℎ = + ℎ 4 can be taken as a measure of the energy error, and consequently, constitutes a possible objective function to minimize. The previous analysis can be also carried out for a composition method (10), resulting in Δ = 2 5 + 2 14 + 2 113 + 2 1112 + 2 23 + 2 212 .
The -stage methods XB whose coefficients are collected in Table 3 have been obtained by minimizing 3 with (29) and in addition provide small values for (27) when applied with ℎ = [ ] ℎ • [ ] ℎ . We should emphasize again that, although methods XB have been obtained by minimizing (29), and thus the local error in the energy, their applicability is by no means limited to Hamiltonian systems. As a matter of fact, both classes of schemes XA and XB can be used with any first-order basic method and its adjoint. Their efficiency may depend, of course, of the type of problem one is approximating and the particular basic scheme taken to form the composition. Moreover, due to the close relationship between symplectic and composition methods, these schemes can also be seen as symplectic partitioned Runge-Kutta methods that, in contrast to splitting schemes, do not require the knowledge of the solution of the elementary flows.

Numerical Examples
Although optimization criteria based on the objective functions 1 , 2 and 3 allow one in principle to construct efficient composition schemes, it is clear that their overall performance depends very much on the particular problem considered, the initial conditions, etc. It is, then, worth considering some illustrative numerical examples to test the methods proposed here with respect to other integrators previously available in the literature. In particular, we take as representatives the splitting method (9) designed in [12] for problems separated into three parts (referred to as ABC 21 in the sequel) and the splitting scheme of [23] considered as a composition (10) (referred as S 6 in Table 2). When a specific composition method (10) is applied to a particular problem of the forṁ = + + and the first-order method is , the implementation is in fact very similar as for a splitting method of the form (9). Thus, in particular, for the integrator (21) one has to apply the following procedure for the time step ⟼ +1 , where one has to take into account the symmetry of the coefficients: 12 = 1 , etc. and = 6: end +1 = It is worth remarking that the examples considered here have been chosen because they admit an straightforward separation into three parts that are explicitly solvable and thus may be used as a kind of testing bench to illustrate the main features of the proposed algorithms. Of course, many other systems could also be considered, including non linear oscillators and the time integration of Vlasov-Maxwell equations [8,20]. In addition, the general technique proposed in [26] for obtaining explicit symplectic approximations of non-separable Hamiltonians provides in a natural manner examples of systems separable into three parts.

Motion of a Charged Particle under Lorentz Force
Neglecting relativistic effects, the evolution of a particle of mass and charge in a given electromagnetic field is described by the Lorentz force as where and denote the electric and magnetic field, respectively. In terms of position and velocity, the equation of motion (30) can be restated aṡ where = − ∕ is the local cyclotron frequency, = ‖ ‖ and = ∕ is the unit vector in the direction of the magnetic field. For simplicity, we assume that both and only depend on the position . System (31) can be split into three parts in such a way that (a) each subpart is explicitly solvable and (b) the volume form in the space ( , ) is exactly preserved [9,27]: The corresponding flows with initial condition ( , 0 ) are given by As in [9], we consider a static, non-uniform electromagnetic field derived from the potentials = 0.01 , = 2 3 respectively, in cylindrical coordinates ( , , ) and with the appropriate normalization. Then, it can be shown that both the angular momentum and energy are invariants of the problem [9]. With = −1, = 1 and starting from the initial position 0 = (0, −1, 0) with initial velocity 0 = (0.10, 0.01, 0), we integrate with the different numerical schemes until the final time = 200 and compute the error in energy and angular momentum along the integration interval. As reference solution we take the output generated by the standard routine DOP853 based on a Runge-Kutta method of order 8 with local error estimation and step size control (with a very stringent tolerance) [28]. In this way, we obtain Figure 1 (top and bottom, respectively), where this error is depicted in terms of the number of the computed sub-flows (by taking different time-steps). For clarity, here and in the sequel, in the left panel we include the results attained by the most efficient XA method, whereas the right panel corresponds to the XB schemes. For reference and comparison, we include in all cases the splitting method (9) proposed in [12] (denoted here as ABC 21 ) and the scheme S 6 , whose coefficients are collected in Table 2. We notice that applying the composition methods proposed here leads to more accurate results than the direct approach based on the splitting methods of Section 2 with the same computational cost, and that the new scheme XB 6 is slightly more efficient that the the splitting scheme S 6 (the remaining composition methods of Tables 2 and 3 provide results between ABC 21 and the best composition depicted here).
In Figure 2 we show the corresponding results obtained by each method for the error in the ( , ) space. One should notice that, although this system is Hamiltonian, the Hamiltonian function is not separable into kinetic plus potential energy, and thus general symplectic Runge-Kutta methods cannot be explicit [27]. In order to use explicit methods, one has to split the system into three parts. On the other hand, all the methods tested here are volume-preserving in the ( , ) space, just as the exact flow.

Disordered Discrete Nonlinear Schrödinger Equation
The Hamiltonian of the disordered discrete nonlinear Schrödinger equation (DDNLS) describes a one-dimensional chain of couples nonlinear oscillators [7]. Here the sum extends over oscillators, are complex variables, ≥ 0 stands for the nonlinearity strength and the random energies are chosen uniformly from the interval [− ∕2, ∕2], where is related with the disorder strength. This model has two invariants: the energy (35) and the norm and has been used to determine how the energy spreads in disordered systems [29]. Rather than analyzing the rich dynamics this system possesses, our interest here is to use (35) as a non-trivial test bench for the integrators we presented in previous sections. By introducing the new (real) generalized coordinates and momenta ( , ) related with through the Hamiltonian function (35) can be written as in such a way that is the sum of three explicitly solvable parts, = + + , with The corresponding flows are given, respectively, by To compare the performance of the numerical integrators previously considered, we take a lattice of = 1000 sites and fixed boundary conditions, 0 = 0 = +1 = +1 = 0. As in [7,30], we excite, at the initial time = 0, 21 central sites by taking the at random in the interval [0, 1] and the respective in such a way that each site has the same constant norm 1, so that the total norm of the system is = 21. Moreover, = 0.72, = 4 and the random disorder parameters are chosen so that the total energy is ≈ −29.63. As in the previous example, we integrate until the final time = 10 and compute the maximum relative error in energy and in norm along the integration interval. The results are depicted in Figure 3, with the top diagrams corresponding to the error in energy and the bottom to the error in norm. The same notation has been used for the tested methods. Finally, in Figure 4 we collect the error in the phase space. As before, the reference solution is obtained with the DOP853 routine. Notice that for this non trivial example the new schemes XA 4 and especially XB 6 show a better efficiency than S 6 , and not only with respect to the preservation of the invariants, but also in the computation of trajectories.

Concluding Remarks
In this work we have presented two different families of fourth-order composition methods especially designed for problems that can be separated into three parts in such a way that each part is explicitly solvable. In addition to the usual optimization criteria applied in the literature to choose the free parameters in the composition, we have introduced another one especially oriented to problems where the energy is a constant of motion. The schemes constructed in this way show an improved behavior, and in fact one of the methods exhibits a superior performance to the familiar scheme 6 of Table 2 on the tested examples. Other relevant examples include certain nonlinear oscillators, Poisson-Maxwell equations arising in plasma physics, and the treatment of non separable Hamiltonian dynamical systems [26].
Although only problems separable into three parts have been considered here, it is clear that the schemes we have introduced can also be applied to differential equations split into any number of pieces ≥ 3. The only modification one requires is to formulate the corresponding first order scheme ℎ and its adjoint * ℎ . One should be aware, however, that augmenting the number leads to evaluating an increasingly large number of flows for methods with large values of , with the subsequent deterioration in performance.
An important topic not addressed in this study concerns the stability of the proposed methods. Typically, for a given method there exists a critical step size ℎ * such that it will be unstable for |ℎ| > ℎ * . Of course, one is interested in methods with ℎ * as large as possible. The linear stability of splitting methods has been analyzed in particular in [31,32], where highly efficient schemes with optimal stability polynomials have presented for numerically approximating the evolution of linear problems. In the nonlinear case, however, the situation is more involved. In [19], a crude measure of the nonlinear stability of a given time symmetric scheme of order is proposed, taking into account the error terms of orders + 1 and + 3. The stability of splitting methods in the particular setting of (semidiscretized) partial differential equations with stiff terms have been considered, in particular, in [33,34]. A theorem is presented [34] concerning the stability of operator-splitting methods applied to linear reaction-diffusion equations with indefinite reaction terms which controls both low and high wave number instabilities. In any case, this result only affects methods up to order 2 with real and positive coefficients, whereas the application of splitting and composition methods of higher order with real coefficients in this setting leads to severe instabilities due to the existence of negative coefficients. The methods we have presented here are aimed at non-stiff problems, and they do not exhibit, at least for the examples we have considered, special step size restrictions in comparison with other splitting methods from the literature.
Finally, it is worth remarking that the local error estimators for composition methods proposed in [35] based on the construction of lower order schemes obtained at each step as a linear combination of the intermediate stages of the main integrator, can also be used in this setting. As a consequence, it is quite straightforward to implement the methods presented here with a variable step size strategy if necessary.