Composition Methods for Dynamical Systems Separable into Three Parts
CComposition Methods for Dynamical Systems Separable intoThree Parts
Fernando Casas ∗ Alejandro Escorihuela-Tomàs † March 8, 2020
Abstract
New families of fourth-order composition methods for the numerical integration of initialvalue problems defined by ordinary differential equations are proposed. They are designed whenthe problem can be separated into three parts in such a way that each part is explicitly solv-able. The methods are obtained by applying different optimization criteria and preserve geomet-ric properties of the continuous problem by construction. Different numerical examples exhibittheir improved performance with respect to previous splitting methods in the literature.Institut de Matemàtiques i Aplicacions de Castelló (IMAC) and Departament de Matemàtiques,Universitat Jaume I, E-12071 Castellón, Spain.
Splitting methods are particularly useful for the numerical integration of ordinary differential equa-tions (ODEs) ̇𝑥 ≡ 𝑑𝑥𝑑𝑡 = 𝑓 ( 𝑥 ) , 𝑥 ( 𝑡 ) = 𝑥 ∈ ℝ 𝐷 (1)when the vector field 𝑓 can be written as 𝑓 ( 𝑥 ) = ∑ 𝑛𝑖 =1 𝑓 𝑖 ( 𝑥 ) , so that each subproblem ̇𝑥 = 𝑓 𝑖 ( 𝑥 ) , 𝑥 ( 𝑡 ) = 𝑥 , 𝑖 = 1 , … , 𝑛 is explicitly solvable, with solution 𝑥 ( 𝑡 ) = 𝜑 [ 𝑖 ] 𝑡 ( 𝑥 ) . Then, by composing the different flows withappropriate 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 numer-ical 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, andin fact some of the most efficient geometric integrators are based on the related ideas of splitting andcomposition [3].In GNI the goal is to construct numerical integrators in such a way that the approximations theyfurnish share one or several qualitative (often, geometric) properties with the exact solution of thedifferential 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 general-purpose methods. In this sense, symplectic integration algorithms for Hamiltonian systems constitutea paradigmatic example of geometric integrators [5, 6]. Splitting and composition methods are widely ∗ Email:
[email protected] † Email: [email protected] a r X i v : . [ m a t h . NA ] J un sed in GNI because the composition of symplectic (or volume preserving, orthogonal, etc.) trans-formations is again symplectic (volume preserving, orthogonal, etc., respectively). In compositionmethods the numerical scheme is constructed as the composition of several simpler integrators forthe problem at hand, so as to improve their accuracy.When 𝑓 in (1) can be separated into two parts, very efficient splitting schemes have been designedand applied to solve a wide variety of problems arising in several fields, ranging from HamiltonianMonte Carlo techniques to the evolution of the 𝑁 -body gravitational problem in Celestial Mechanics(see [3, 4] and references therein).There are, however, relevant problems in applications where 𝑓 has to be decomposed into threeor more parts in order to have subproblems that are explicitly solvable. Examples include the disor-dered discrete nonlinear Schrödinger equation [7], Vlasov–Maxwell equations in plasma physics [8],the motion of a charged particle in an electromagnetic field according with the Lorentz force law [9]and problems in molecular dynamics [10]. In that case, although in principle methods of any order ofaccuracy can be built, the resulting algorithms involve such a large number of maps that they are notcompetitive in practice. It is the purpose of this paper to present an alternative class of efficient meth-ods for the problem at hand and compare their performance on some non-trivial physical examplesthan can be split into three parts.The paper is structured as follows. We first review how splitting methods can be directly appliedto get numerical solutions (Section 2). Then the attention is turned to the application of compositionmethods, and we get a family of 4th-order schemes obtained by applying a standard optimizationprocedure (Section 3). In Section 4 we show how standard splitting methods, when formulated as acomposition scheme, lead to very competitive integrators, and also propose a different optimizationcriterion for systems possessing invariant quantities. This allows us to get a new family of 4th-orderschemes. All these integration algorithms are subsequently tested in Section 5 on a pair of numericalexamples. Finally, Section 6 contains some concluding remarks. In what follows we assume that the vector field 𝑓 in (1) can be split into three parts, 𝑓 ( 𝑥 ) = 𝑓 𝑎 ( 𝑥 ) + 𝑓 𝑏 ( 𝑥 ) + 𝑓 𝑐 ( 𝑥 ) (2)in such a way that the exact ℎ -flows 𝜑 [ 𝑎 ] ℎ , 𝜑 [ 𝑏 ] ℎ , 𝜑 [ 𝑐 ] ℎ , corresponding to 𝑓 𝑎 , 𝑓 𝑏 , 𝑓 𝑐 , respectively, can becomputed exactly.It is clear that the composition 𝜒 ℎ = 𝜑 [ 𝑎 ] ℎ ◦ 𝜑 [ 𝑏 ] ℎ ◦ 𝜑 [ 𝑐 ] ℎ (3)(or any other permutation of the sub-flows) provides a first-order approximation to the exact solution 𝑥 ( ℎ ) = 𝜑 ℎ ( 𝑥 ) of (1), i.e., 𝜒 ℎ ( 𝑥 ) = 𝜑 ℎ ( 𝑥 ) + ( ℎ ) , whereas the so-called Strang splitting [2] ℎ = 𝜑 [ 𝑎 ] ℎ ∕2 ◦ 𝜑 [ 𝑏 ] ℎ ∕2 ◦ 𝜑 [ 𝑐 ] ℎ ◦ 𝜑 [ 𝑏 ] ℎ ∕2 ◦ 𝜑 [ 𝑎 ] ℎ ∕2 (4)leads to a second-order approximation.Higher order approximations to the exact solution of (1) can be obtained by generalizing (4), i.e.,by considering splitting schemes of the form 𝜓 [ 𝑟 ] ℎ = 𝜑 [ 𝑐 ] 𝑐 𝑠 ℎ ◦ 𝜑 [ 𝑏 ] 𝑏 𝑠 ℎ ◦ 𝜑 [ 𝑎 ] 𝑎 𝑠 ℎ ◦ ⋯ ◦ 𝜑 [ 𝑐 ] 𝑐 ℎ ◦ 𝜑 [ 𝑏 ] 𝑏 ℎ ◦ 𝜑 [ 𝑎 ] 𝑎 ℎ , (5)where the coefficients 𝑎 𝑖 , 𝑏 𝑖 , 𝑐 𝑖 , 𝑖 = 1 , … , 𝑠 , are chosen to achieve a prescribed order of accuracy, say, 𝑟 , 𝜓 [ 𝑟 ] ℎ ( 𝑥 ) = 𝜑 ℎ ( 𝑥 ) + ( ℎ 𝑟 +1 ) as ℎ → . (6)2equirement (6) leads to a set of polynomial equations (the so-called order conditions), whose numberand complexity grows enormously with the order. In particular, if 𝑟 = 1 (i.e., for a consistencymethod) one has 𝑠 ∑ 𝑖 =1 𝑎 𝑖 = 1 , 𝑠 ∑ 𝑖 =1 𝑏 𝑖 = 1 , 𝑠 ∑ 𝑖 =1 𝑐 𝑖 = 1 . The specific number of order conditions is determined in fact by the dimension 𝑐 𝑘 of the homogeneoussubspace of grade 𝑘 , ≤ 𝑘 ≤ 𝑟 , of the free Lie algebra ( 𝐴, 𝐵, 𝐶 ) generated by the Lie derivatives 𝐴, 𝐵, 𝐶 corresponding to 𝑓 [ 𝑎 ] , 𝑓 [ 𝑏 ] , 𝑓 [ 𝑐 ] , respectively [1]. These dimensions are collected Table 1 for ≤ 𝑘 ≤ .Table 1: Number of order conditions to be satisfied by a splitting method of the form (5) at eachorder 𝑘 . Grade 𝑘 𝑐 𝑘 order conditionsand therefore the evaluation of at least a similar number of sub-flows to have as many parametersas equations. This number can be reduced by considering time-symmetric methods, i.e., schemesverifying 𝜓 [ 𝑟 ] ℎ ◦ 𝜓 [ 𝑟 ]− ℎ = id , (7)where id is the identity map. Condition (7) is verified by left-right palindromic compositions, i.e., if 𝑎 𝑠 +1− 𝑖 = 𝑎 𝑖 , 𝑏 𝑠 +1− 𝑖 = 𝑏 𝑖 , 𝑐 𝑠 +1− 𝑖 = 𝑐 𝑖 , 𝑖 = 1 , , … in (5). Then all the conditions at even order are automatically satisfied. Thus, a symmetric methodof order 4 requires solving order conditions (instead of 32). Still, within this approach, one has tosolve 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 shownthat if one aims to get schemes (5) of order 2 with the minimum number of maps, then the Strangsplitting (4) is recovered. With respect to order 4, the following scheme was presented: 𝜓 [4] 𝜏 = 𝜑 [ 𝑐 ] 𝑐 𝜏 ◦ 𝜑 [ 𝑏 ] 𝑏 𝜏 ◦ 𝜑 [ 𝑎 ] 𝑎 𝜏 ◦ 𝜑 [ 𝑏 ] 𝑏 𝜏 ◦ 𝜑 [ 𝑐 ] 𝑐 𝜏 ◦ 𝜑 [ 𝑏 ] 𝑏 𝜏 ◦ 𝜑 [ 𝑎 ] 𝑎 𝜏 ◦ 𝜑 [ 𝑏 ] 𝑏 𝜏 ◦ 𝜑 [ 𝑐 ] 𝑐 𝜏 ◦ 𝜑 [ 𝑏 ] 𝑏 𝜏 ◦ 𝜑 [ 𝑎 ] 𝑎 𝜏 ◦ 𝜑 [ 𝑏 ] 𝑏 𝜏 ◦ 𝜑 [ 𝑐 ] 𝑐 𝜏 (8)with 𝑎 = 𝑤 , 𝑎 = 𝑤 , 𝑏 = 𝑏 = 𝑤 , 𝑏 = 𝑤 , 𝑐 = 𝑤 , 𝑐 = 𝑤 + 𝑤 and 𝑤 = 12 − 2 , 𝑤 = 1 − 2 𝑤 . In fact, 13 is the minimum number of maps required. More efficient schemes involving 17 and 25maps can also be found in [11]. For simplicity, we denote method (8) as ( 𝑐 𝑏 𝑎 𝑏 𝑐 𝑏 𝑎 𝑏 𝑐 𝑏 𝑎 𝑏 𝑐 ) .More recently, in [12] a method involving 21 maps of the form ( 𝑎 𝑏 𝑐 𝑎 𝑏 𝑐 𝑎 𝑏 𝑐 𝑎 𝑏 𝑎 𝑐 𝑏 𝑎 𝑐 𝑏 𝑎 𝑐 𝑏 𝑎 ) (9)has also been proposed and tested on several numerical examples.3 Second Approach: Composition Methods
As it is clear from the previous considerations, constructing high order splitting methods for systemsseparable into three parts requires solving a large number of polynomial equations involving the coef-ficients, and this is a very challenging task in general. For this reason, we turn our attention to anotherstrategy based on compositions of the first order 𝜒 ℎ = 𝜑 [ 𝑎 ] ℎ ◦ 𝜑 [ 𝑏 ] ℎ ◦ 𝜑 [ 𝑐 ] ℎ and its adjoint, 𝜒 ∗ ℎ ∶= ( 𝜒 − ℎ ) −1 = 𝜑 [ 𝑐 ] ℎ ◦ 𝜑 [ 𝑏 ] ℎ ◦ 𝜑 [ 𝑎 ] ℎ with appropriately chosen weights. In other words, we look for integrators of the form 𝜓 ℎ = 𝜒 𝛼 𝑠 ℎ ◦ 𝜒 ∗ 𝛼 𝑠 −1 ℎ ◦ ⋯ ◦ 𝜒 𝛼 ℎ ◦ 𝜒 ∗ 𝛼 ℎ , with ( 𝛼 , … , 𝛼 𝑠 ) ∈ ℝ 𝑠 (10)verifying in addition the time-symmetry condition 𝛼 𝑠 +1− 𝑖 = 𝛼 𝑖 for all 𝑖 . Remark 1
Methods of the form 𝜓 ℎ = [2] 𝛼 𝑚 ℎ ◦ [2] 𝛼 𝑚 −1 ℎ ◦ ⋯ ◦ [2] 𝛼 ℎ ◦ [2] 𝛼 ℎ with ( 𝛼 , … , 𝛼 𝑚 ) ∈ ℝ 𝑚 (11) and 𝛼 𝑚 +1− 𝑖 = 𝛼 𝑖 (commonly referred in the literature as symmetric compositions of symmetric meth-ods [1]) verify a much reduced number of order conditions and allows one to construct very efficienthigh-order schemes [3]. Notice that, since the Strang splitting (4) verifies [2] ℎ = 𝜒 ℎ ∕2 ◦ 𝜒 ∗ ℎ ∕2 , then itis clear that when analyzing methods (10) we also recover schemes of the form (11). The analysis of the composition methods considered here can be conveniently done by considering theLie operators associated with the vector fields involved and the graded free Lie algebra they generate.As is well known, for each infinitely differentiable map 𝑔 ∶ ℝ 𝐷 ⟶ ℝ , the function 𝑔 ( 𝜑 ℎ ( 𝑥 )) admits an expansion of the form [13, 5] 𝑔 ( 𝜑 ℎ ( 𝑥 )) = exp( ℎ𝐹 )[ 𝑔 ]( 𝑥 ) = 𝑔 ( 𝑥 ) + ∑ 𝑘 ≥ ℎ 𝑘 𝑘 ! 𝐹 𝑘 [ 𝑔 ]( 𝑥 ) , 𝑥 ∈ ℝ 𝐷 , where 𝐹 is the Lie derivative associated with 𝑓 , 𝐹 ≡ 𝐿 𝑓 = 𝐷 ∑ 𝑖 =1 𝑓 𝑖 ( 𝑥 ) 𝜕𝜕𝑥 𝑖 . (12)Analogously, for the basic method 𝜒 ℎ one can associate a series of linear operators so that [14] 𝑔 ( 𝜒 ℎ ( 𝑥 )) = exp( 𝑌 ( ℎ ))[ 𝑔 ]( 𝑥 ) , with 𝑌 ( ℎ ) = ∑ 𝑘 ≥ ℎ 𝑘 𝑌 𝑘 for all functions 𝑔 , whereas for its adjoint one has 𝑔 ( 𝜒 ∗ ℎ ( 𝑥 )) = exp ( − 𝑌 (− ℎ ) ) [ 𝑔 ]( 𝑥 ) . Then the operator series associated with the integrator (10) is Ψ( ℎ ) = exp(− 𝑌 (− ℎ𝛼 )) exp( 𝑌 ( ℎ𝛼 )) ⋯ exp(− 𝑌 (− ℎ𝛼 𝑠 −1 )) exp( 𝑌 ( ℎ𝛼 𝑠 )) . Notice that the order of the operators is the reverse of the maps in (10) ([3] p. 88). Now, byrepeated application of the Baker–Campbell–Hausdorff formula [4] we can express formally Ψ( ℎ ) asthe exponential of an operator ̃𝐹 ( ℎ ) , Ψ( ℎ ) = exp( ̃𝐹 ( ℎ )) , with ̃𝐹 ( ℎ ) = ∑ 𝑘 ≥ ℎ 𝑘 𝐹 𝑘 , (13)4 𝑘 𝐹 𝑘 ∈ 𝑘 for each 𝑘 ≥ and = ⨁ 𝑘 ≥ 𝑘 is the graded free Lie algebra generated by the operators { ℎ𝑌 , ℎ 𝑌 , ℎ 𝑌 , …} , where, by consistency, 𝑌 = 𝐹 . One has explicitly 𝑌 ( ℎ𝛼 𝑖 ) = ℎ𝛼 𝑖 𝑌 + ( ℎ𝛼 𝑖 ) 𝑌 + ( ℎ𝛼 𝑖 ) 𝑌 + ⋯ − 𝑌 (− ℎ𝛼 𝑖 ) = ℎ𝛼 𝑖 𝑌 − ( ℎ𝛼 𝑖 ) 𝑌 + ( ℎ𝛼 𝑖 ) 𝑌 − ⋯ so that ̃𝐹 ( ℎ ) = ℎ𝑤 𝑌 + ℎ 𝑤 𝑌 + ℎ ( 𝑤 𝑌 + 𝑤 [ 𝑌 , 𝑌 ])+ ℎ ( 𝑤 𝑌 + 𝑤 [ 𝑌 , 𝑌 ] + 𝑤 [ 𝑌 , [ 𝑌 , 𝑌 ]]) (14) + ℎ ( 𝑤 𝑌 + 𝑤 [ 𝑌 , 𝑌 ] + 𝑤 [ 𝑌 , 𝑌 , 𝑌 ]+ 𝑤 [ 𝑌 , 𝑌 , 𝑌 , 𝑌 ] + 𝑤 [ 𝑌 , 𝑌 ] + 𝑤 [ 𝑌 , 𝑌 , 𝑌 ] ) + ( ℎ ) , where [ 𝑌 , 𝑌 , 𝑌 ] ≡ [ 𝑌 , [ 𝑌 , 𝑌 ]] , etc, [ ⋅ , ⋅ ] refers to the usual Lie bracket and 𝑤 , 𝑤 , … are polyno-mials in the coefficients 𝛼 𝑖 . In particular, one has 𝑤 = 𝑠 ∑ 𝑖 =1 𝛼 𝑖 , 𝑤 = 𝑠 ∑ 𝑖 =1 (−1) 𝑖 𝛼 𝑖 ,𝑤 = 𝑠 ∑ 𝑖 =1 𝛼 𝑖 , 𝑤 = 𝑠 ∑ 𝑖 =1 (−1) 𝑖 𝛼 𝑖 ,𝑤 = 12 ( 𝑠 −1 ∑ 𝑖 =1 (−1) 𝑖 +1 𝛼 𝑖 𝑠 ∑ 𝑗 = 𝑖 +1 𝛼 𝑗 + 𝑠 −1 ∑ 𝑖 =1 𝛼 𝑖 𝑠 ∑ 𝑗 = 𝑖 +1 (−1) 𝑗 𝛼 𝑗 ) . (15)Thus, a time-symmetric 4th-order method has to satisfy only consistency ( 𝑤 = 1 ) and the orderconditions at order three, 𝑤 = 𝑤 = 0 . Notice, then, that the minimum number of maps to beconsidered is 𝑠 = 3 . In that case the integrator reads 𝜓 ℎ = 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 (16)and the unique (real) solution is given by 𝛼 = 𝛼 = 12(2 − 2 ) , 𝛼 = 12 − 2 𝛼 . This scheme corresponds to the familiar triple-jump integrator [15] 𝜓 ℎ = [2] 𝛼ℎ ∕2 ◦ [2] 𝛽ℎ ◦ [2] 𝛼ℎ ∕2 with 𝛼 = 1∕(2 − 2 ) . (17)If 𝜒 ℎ = 𝜑 [ 𝑎 ] ℎ ◦ 𝜑 [ 𝑏 ] ℎ ◦ 𝜑 [ 𝑐 ] ℎ , then 𝜓 ℎ involves 13 maps (the minimum number) and corresponds pre-cisely to the splitting method (8).It is worth remarking that the order conditions (15) are general for any composition method of theform (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 Eulermethod as 𝜒 ℎ and the implicit Euler method as 𝜒 ∗ ℎ , and also a symplectic semi-implicit method andits adjoint, leading to the symplectic partitioned Runge–Kutta schemes considered in [16]. 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 con-sists in adding more stages in the composition and determine the extra free parameters thus introduced5ccording with some optimization criteria. Although assessing the quality of a given integrationmethod applied to all initial value problem is by no means obvious (the dominant error terms are notnecessarily the same for different problems), several strategies have been proposed along the years tofix these free parameters in the composition method (10). Thus, in particular, one looks for solutionssuch that the absolute value of the coefficients, i.e., 𝐸 ( 𝜶 ) = 𝑠 ∑ 𝑖 =1 | 𝛼 𝑖 | (18)is as small as possible, the logic being that higher order terms in the expansion (14) involve powersof these coefficients. In fact, methods with small values of 𝐸 ( 𝜶 ) usually have large stability domainsand small error terms [1]. In addition, for a number of problems, the dominant error term is preciselythe coefficient 𝑤 multiplying 𝑌 in the expansion (14), so that it makes sense to minimize 𝐸 ( 𝜶 ) = 2 𝑠 ||| 𝑠 ∑ 𝑖 =1 𝛼 𝑖 ||| , (19)for a given composition to take also into account the computational effort measured as the number 𝑠 of basic schemes considered. Here, as in [17], we construct symmetric methods with small values of 𝐸 which, in addition, have also small values of 𝐸 . For future reference, the corresponding valuesof the objective functions for the triple-jump (17) are 𝐸 = 4 . and 𝐸 = 4 . , respectively.Next we collect the most efficient schemes we have obtained with 𝑠 = 4 , , by applying thisstrategy. 𝑠 = 4 stages. The composition is 𝜓 ℎ = 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 , (20)and involves 17 maps when the basic scheme 𝜒 ℎ is given by (3). Now we have a free parameter,which we take as 𝛼 . The minima of both 𝐸 and 𝐸 are achieved at approximately 𝛼 = 0 . , andthe resulting coefficients are collected in Table 2 as method XA . In that case, 𝐸 = 2 . and 𝐸 = 3 . . 𝑠 = 5 stages. The resulting composition 𝜓 ℎ = 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 involves 21 maps when applied to a system separable into three parts. Minimum values for 𝐸 and 𝐸 are achieved when 𝛼 = 𝛼 = 𝛼 = 𝛼 = 12(4 − 4 ) , 𝛼 = 12 − 4 𝛼 . In consequence, the method can be written as 𝜓 ℎ = [2] 𝛼ℎ ◦ [2] 𝛼ℎ ◦ [2] 𝛽ℎ ◦ [2] 𝛼ℎ ◦ [2] 𝛼ℎ with 𝛼 = 2 𝛼 , 𝛽 = 2 𝛼 . Then 𝐸 = 2 . and 𝐸 = 2 . . This method, denoted XA , was firstproposed in [18] and analyzed in detail in [19]. 6 = 6 stages. Analogously we have considered a composition involving three free parameters (and25 maps when 𝜒 ℎ is given by (3)): 𝜓 ℎ = 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 ◦ 𝜒 𝛼 ◦ 𝜒 ∗ 𝛼 . (21)The proposed solution is collected in Table 2 as method XA leading to 𝐸 = 2 . , 𝐸 =2 . . Notice how, by increasing the number of stages, it is possible to reduce the value of 𝐸 and 𝐸 as a measure of the efficiency of the schemes. This integrator has been tested in the numericalintegration of the so-called reduced Vlasov–Maxwell system [20].We could of course increase the number of stages. It turns out, however, that with 𝑠 = 7 onehas the sufficient number of parameters to satisfy all the order conditions up to order 6, resulting in amethod of the form (11) [15] involving 29 maps. More efficient 6th-order schemes can be obtainedindeed by increasing the number of stages. Thus, in particular, with 𝑠 = 9 and 𝑠 = 11 one has themethods designed in [21] (37 maps) and [22] (45 maps), respectively, when the basic scheme is givenby (3).Table 2: Fourth-order composition methods XA 𝑠 with 𝑠 stages minimizing 𝐸 and 𝐸 . Method S correspondsto the splitting method of ([23] Table 2) expressed as a composition scheme. XA 𝛼 = 0 . 𝛼 = −0 . 𝛼 = 0 . 𝛼 = 0 . XA 𝛼 = 𝛼 = 𝛼 = 𝛼 = 12(4 − 4 ) 𝛼 = − 4 𝛼 XA 𝛼 = 0 . 𝛼 = 0 . 𝛼 = 0 . 𝛼 = −0 . 𝛼 = 0 . 𝛼 = 0 . S 𝛼 = 0 . 𝛼 = 0 . 𝛼 = 0 . 𝛼 = −0 . 𝛼 = 0 . 𝛼 = 0 . 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]. Letus assume that 𝑓 in the ODE (1) can be split into two parts, ̇𝑥 = 𝑓 𝑎 ( 𝑥 ) + 𝑓 𝑏 ( 𝑥 ) , which each partexplicitly solvable, and take 𝜒 ℎ = 𝜑 [ 𝑏 ] ℎ ◦ 𝜑 [ 𝑎 ] ℎ . Then, the adjoint method reads 𝜒 ∗ ℎ = 𝜑 [ 𝑎 ] ℎ ◦ 𝜑 [ 𝑏 ] ℎ and thecomposition (10) adopts the form 𝜓 ℎ = ( 𝜑 [ 𝑏 ] 𝛼 𝑠 ℎ ◦ 𝜑 [ 𝑎 ] 𝛼 𝑠 ℎ ) ◦ ( 𝜑 [ 𝑎 ] 𝛼 𝑠 −1 ℎ ◦ 𝜑 [ 𝑏 ] 𝛼 𝑠 −1 ℎ ) ◦ ⋯ ◦ ( 𝜑 [ 𝑏 ] 𝛼 ℎ ◦ 𝜑 [ 𝑎 ] 𝛼 ℎ ) ◦ ( 𝜑 [ 𝑎 ] 𝛼 ℎ ◦ 𝜑 [ 𝑏 ] 𝛼 ℎ ) . (22)Since 𝜑 [ 𝑖 ] ℎ , 𝑖 = 𝑎, 𝑏 are exact flows, then they verify 𝜑 [ 𝑖 ] 𝛽ℎ ◦ 𝜑 [ 𝑖 ] 𝛿ℎ = 𝜑 [ 𝑖 ]( 𝛽 + 𝛿 ) ℎ , and (22) can be rewrittenas the splitting scheme 𝜓 ℎ = 𝜑 [ 𝑏 ] 𝑏 𝑠 +1 ℎ ◦ 𝜑 [ 𝑎 ] 𝑎 𝑠 ℎ ◦ 𝜑 [ 𝑏 ] 𝑏 𝑠 ℎ ◦ ⋯ ◦ 𝜑 [ 𝑏 ] 𝑏 ℎ ◦ 𝜑 [ 𝑎 ] 𝑎 ℎ ◦ 𝜑 [ 𝑏 ] 𝑏 ℎ (23)7f 𝑏 = 𝛼 and 𝑎 𝑗 = 𝛼 𝑗 + 𝛼 𝑗 −1 , 𝑏 𝑗 +1 = 𝛼 𝑗 +1 + 𝛼 𝑗 , 𝑗 = 1 , … , 𝑠 (24)(with 𝛼 𝑠 +1 = 0 ). Conversely, any integrator of the form (23) with ∑ 𝑠𝑖 =1 𝑎 𝑖 = ∑ 𝑠 +1 𝑖 =1 𝑏 𝑖 can be expressedin the form (10) with 𝜒 ℎ = 𝜑 [ 𝑏 ] ℎ ◦ 𝜑 [ 𝑎 ] ℎ and 𝛼 𝑠 = 𝑏 𝑠 +1 ,𝛼 𝑗 −1 = 𝑎 𝑗 − 𝛼 𝑗 , 𝛼 𝑗 −2 = 𝑏 𝑗 − 𝛼 𝑗 −1 , 𝑗 = 𝑠, 𝑠 − 1 , … , , with 𝛼 = 0 for consistency. In consequence, any splitting method in principle designed for sys-tems of the form ̇𝑥 = 𝑓 𝑎 ( 𝑥 ) + 𝑓 𝑏 ( 𝑥 ) with no further restrictions on 𝑓 𝑎 or 𝑓 𝑏 can be formulated asa composition (10) which, in turn, can also be applied when 𝑓 is split into three (or more) pieces, 𝑓 = 𝑓 𝑎 + 𝑓 𝑏 + 𝑓 𝑐 , by taking 𝜒 ℎ = 𝜑 [ 𝑎 ] ℎ ◦ 𝜑 [ 𝑏 ] ℎ ◦ 𝜑 [ 𝑐 ] ℎ . 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 partshas been presented in ([23] Table 2) (method S ) and will be used in our numerical tests. It is atime-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 theform (21), with coefficients 𝛼 𝑖 listed in Table 2. For comparison, the corresponding values of 𝐸 and 𝐸 are 𝐸 = 2 . and 𝐸 = 3 . . 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 as ̇𝑞 𝑖 = 𝜕𝐻𝜕𝑝 𝑖 , ̇𝑝 𝑖 = − 𝜕𝐻𝜕𝑞 𝑖 , 𝑖 = 1 , … , 𝑑 (25)so that 𝑥 = ( 𝑞, 𝑝 ) 𝑇 , 𝑓 = (∇ 𝑝 𝐻, −∇ 𝑞 𝐻 ) 𝑇 ≡ 𝑋 𝐻 and 𝐻 ( 𝑞, 𝑝 ) is the Hamiltonian. The Lie derivativeassociated with 𝑋 𝐻 verifies, for any function 𝐺 ∶ 𝐷 ⊂ ℝ 𝑑 ⟶ ℝ , 𝐿 𝑋 𝐻 𝐺 = −{ 𝐻, 𝐺 } = − 𝑑 ∑ 𝑗 =1 ( 𝜕𝐻𝜕𝑞 𝑗 𝜕𝐺𝜕𝑝 𝑗 − 𝜕𝐺𝜕𝑞 𝑗 𝜕𝐻𝜕𝑝 𝑗 ) . In other words, { 𝐻, 𝐺 } is the Poisson bracket of 𝐻 and 𝐺 . In this context, then, the Lie bracketof 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 thetotal energy of the system. If 𝐻 can be split as 𝐻 = 𝐴 + 𝐵 , then 𝑓 [ 𝑎 ] = 𝐿 𝑋 𝐴 , 𝑓 [ 𝑏 ] = 𝐿 𝑋 𝐵 and thesplitting method (23) is also symplectic. Important as it is that the method shares this feature with theexact flow, one would like in addition that the energy be preserved as accurately as possible (since anumerical scheme cannot preserve both the symplectic form and the energy). A possible optimizationcriterion would be then to select the free parameters in such a way that the error in the energy (or morein 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-nian ̃𝐻 ℎ in the limit ℎ → for a 4th-order splitting method (23). A straightforward calculation showsthat ̃𝐻 ℎ = 𝐻 + ℎ 𝑘 , { 𝐴, 𝐴, 𝐴, 𝐴, 𝐵 } + ℎ 𝑘 , { 𝐵, 𝐴, 𝐴, 𝐴, 𝐵 } + ℎ 𝑘 , { 𝐴, 𝐴, 𝐵, 𝐴, 𝐵 }+ ℎ 𝑘 , { 𝐴, 𝐵, 𝐵, 𝐴, 𝐵 } + ℎ 𝑘 , { 𝐵, 𝐴, 𝐵, 𝐴, 𝐵 } + ℎ 𝑘 , { 𝐵, 𝐵, 𝐵, 𝐴, 𝐵 }++ ℎ ∑ 𝑗 =1 𝑘 ,𝑗 𝐸 ,𝑗 + ( ℎ ) , (26)8here 𝑘 𝑖,𝑗 are polynomials in the coefficients 𝑎 𝑗 , 𝑏 𝑗 , { 𝐴, 𝐴, 𝐴, 𝐴, 𝐵 } refers to the iterated Poissonbracket { 𝐴, { 𝐴, { 𝐴, { 𝐴, 𝐵 }}}} , and 𝐸 ,𝑗 are (independent) Poisson brackets involving 6 functions 𝐴 and 𝐵 .Now the Lie formalism allows one to get the Taylor expansion of the energy after one time-step([5], Section 12.2) as 𝐻 ( 𝑞 𝑖 +1 , 𝑝 𝑖 +1 ) = exp(− ℎ ̃𝐻 ℎ ) 𝐻 ( 𝑞 𝑖 , 𝑝 𝑖 ) = 𝐻 ( 𝑞 𝑖 , 𝑝 𝑖 ) − ℎ ̃𝐻 ℎ 𝐻 ( 𝑞 𝑖 , 𝑝 𝑖 ) + 12 ℎ ̃𝐻 ℎ 𝐻 ( 𝑞 𝑖 , 𝑝 𝑖 ) + ⋯ , where ̃𝐻 ℎ ( ⋅ ) = { ̃𝐻 ℎ , ⋅ } .An elementary calculation shows that 𝐻 ( 𝑞 𝑖 +1 , 𝑝 𝑖 +1 ) − 𝐻 ( 𝑞 𝑖 , 𝑝 𝑖 ) = ℎ ( 𝑘 𝐸 + ( 𝑘 − 𝑘 ) 𝐸 + ( 𝑘 − 𝑘 ) 𝐸 + 𝑘 𝐸 + ( 𝑘 − 13 𝑘 ) 𝐸 + ( 𝑘 − 13 𝑘 ) 𝐸 + ( 𝑘 + 𝑘 ) 𝐸 + ( 𝑘 − 𝑘 ) 𝐸 + 𝑘 𝐸 ) + ( ℎ ) . Thus, for small ℎ , Δ ≡ 𝑘 + ( 𝑘 − 𝑘 ) + ( 𝑘 − 𝑘 ) + 𝑘 + ( 𝑘 − 13 𝑘 ) + ( 𝑘 − 13 𝑘 ) + ( 𝑘 + 𝑘 ) + ( 𝑘 − 𝑘 ) + 𝑘 (27)can be taken as a measure of the energy error, and consequently, 𝐸 = 2 𝑠 Δ (28)constitutes a possible objective function to minimize. The previous analysis can be also carried outfor a composition method (10), resulting in Δ = 𝑤 + 𝑤 + 𝑤 + 𝑤 + 𝑤 + 𝑤 . (29)The 𝑠 -stage methods XB 𝑠 whose coefficients are collected in Table 3 have been obtained by mini-mizing 𝐸 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 sys-tems. As a matter of fact, both classes of schemes XA 𝑠 and XB 𝑠 can be used with any first-orderbasic method and its adjoint. Their efficiency may depend, of course, of the type of problem one isapproximating and the particular basic scheme taken to form the composition. Moreover, due to theclose relationship between symplectic and composition methods, these schemes can also be seen assymplectic partitioned Runge–Kutta methods that, in contrast to splitting schemes, do not require theknowledge of the solution of the elementary flows.9able 3: Fourth-order composition methods XB 𝑠 with 𝑠 stages minimizing 𝐸 . XB 𝛼 = 0 . 𝛼 = 0 . 𝛼 = −0 . 𝛼 = 0 . XB 𝛼 = 0 . 𝛼 = 0 . 𝛼 = 0 . 𝛼 = −0 . 𝛼 = 0 . XB 𝛼 = 120 𝛼 = 71660 𝛼 = 47330 𝛼 = 37165 𝛼 = − 313660 𝛼 = 511 Although optimization criteria based on the objective functions 𝐸 , 𝐸 and 𝐸 allow one in principleto construct efficient composition schemes, it is clear that their overall performance depends verymuch on the particular problem considered, the initial conditions, etc. It is, then, worth consideringsome illustrative numerical examples to test the methods proposed here with respect to other integra-tors 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 in the sequel) andthe splitting scheme of [23] considered as a composition (10) (referred as S in Table 2).When a specific composition method (10) is applied to a particular problem of the form ̇𝑥 = 𝑓 𝑎 + 𝑓 𝑏 + 𝑓 𝑐 and the first-order method is 𝜒 ℎ = 𝜑 [ 𝑎 ] ℎ ◦ 𝜑 [ 𝑏 ] ℎ ◦ 𝜑 [ 𝑐 ] ℎ , the implementation is in fact verysimilar as for a splitting method of the form (9). Thus, in particular, for the integrator (21) one has toapply the following procedure for the time step 𝑥 𝑛 ⟼ 𝑥 𝑛 +1 , where one has to take into account thesymmetry of the coefficients: 𝛼 = 𝛼 , etc. and 𝑠 = 6 : 𝑦 = 𝑥 𝑛 do 𝑗 = 1 ∶ 6 𝑦 = 𝜑 [ 𝑎 ] 𝛼 𝑗 −1 ℎ 𝑦𝑦 = 𝜑 [ 𝑏 ] 𝛼 𝑗 −1 ℎ 𝑦̃𝛼 = 𝛼 𝑗 −1 + 𝛼 𝑗 𝑦 = 𝜑 [ 𝑐 ] ̃𝛼ℎ 𝑦𝑦 = 𝜑 [ 𝑏 ] 𝛼 𝑗 ℎ 𝑦𝑦 = 𝜑 [ 𝑎 ] 𝛼 𝑗 ℎ 𝑦 end 𝑥 𝑛 +1 = 𝑦 It is worth remarking that the examples considered here have been chosen because they admit anstraightforward separation into three parts that are explicitly solvable and thus may be used as a kindof testing bench to illustrate the main features of the proposed algorithms. Of course, many othersystems 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 explicit10ymplectic approximations of non-separable Hamiltonians provides in a natural manner examples ofsystems separable into three parts.
Neglecting relativistic effects, the evolution of a particle of mass 𝑚 and charge 𝑞 in a given electro-magnetic field is described by the Lorentz force as 𝑚 ̈ 𝐱 = 𝑞 ( 𝐄 + ̇ 𝐱 × 𝐁 ) , (30)where 𝐄 and 𝐁 denote the electric and magnetic field, respectively. In terms of position and velocity,the equation of motion (30) can be restated as ̇ 𝐱 = 𝐯 ̇ 𝐯 = 𝑞𝑚 𝐄 + 𝜔 𝐛 × 𝐯 (31)where 𝜔 = − 𝑞𝐵 ∕ 𝑚 is the local cyclotron frequency, 𝐵 = ‖ 𝐁 ‖ and 𝐛 = 𝐁 ∕ 𝐵 is the unit vector inthe direction of the magnetic field. For simplicity, we assume that both 𝐄 and 𝐁 only depend on theposition 𝐱 .System (31) can be split into three parts in such a way that (a) each subpart is explicitly solvableand (b) the volume form in the space ( 𝐱 , 𝐯 ) is exactly preserved [9, 27]: 𝑑𝑑𝑡 ( 𝐱𝐯 ) = ( 𝐯 ) + ( 𝑞𝑚 𝐄 ( 𝐱 ) ) + ( 𝜔 ( 𝐱 ) 𝐛 ( 𝐱 ) × 𝐯 ) = 𝑓 [ 𝑎 ] ( 𝐱 , 𝐯 ) + 𝑓 [ 𝑏 ] ( 𝐱 , 𝐯 ) + 𝑓 [ 𝑐 ] ( 𝐱 , 𝐯 ) . (32)The corresponding flows with initial condition ( 𝐱 , 𝐯 ) are given by 𝜑 [ 𝑎 ] 𝑡 ∶ { 𝐱 ( 𝑡 ) = 𝐱 + 𝑡 𝐯 𝐯 ( 𝑡 ) = 𝐯 , 𝜑 [ 𝑏 ] 𝑡 ∶ { 𝐱 ( 𝑡 ) = 𝐱 𝐯 ( 𝑡 ) = 𝐯 + 𝑡 𝑞𝑚 𝐄 ( 𝐱 ) 𝜑 [ 𝑐 ] 𝑡 ∶ { 𝐱 ( 𝑡 ) = 𝐱 𝐯 ( 𝑡 ) = exp( 𝑡𝜔 ( 𝐱 ) ̂ 𝐛 ) 𝐯 (33)where ̂ 𝐛 ≡ ̂ 𝐛 ( 𝐱 ) is the skew-symmetric matrix ̂ 𝐛 ( 𝐱 ) = ⎛⎜⎜⎝ 𝑏 ( 𝐱 ) 𝑏 ( 𝐱 ) 𝑏 ( 𝐱 ) 0 − 𝑏 ( 𝐱 )− 𝑏 ( 𝐱 ) 𝑏 ( 𝐱 ) 0 ⎞⎟⎟⎠ associated with 𝐛 ( 𝐱 ) = ( 𝑏 ( 𝐱 ) , 𝑏 ( 𝐱 ) , 𝑏 ( 𝐱 )) 𝑇 .As in [9], we consider a static, non-uniform electromagnetic field 𝐄 = −∇ 𝑉 = 0 . 𝑟 ( 𝑥 𝐞 𝑥 + 𝑦 𝐞 𝑦 ) , 𝐁 = ∇ × 𝐀 = 𝑟 𝐞 𝑧 (34)derived from the potentials 𝑉 = 0 . 𝑟 , 𝐀 = 𝑟 𝐞 𝜃 respectively, in cylindrical coordinates ( 𝑟, 𝜃, 𝑧 ) and with the appropriate normalization. Then, it canbe shown that both the angular momentum and energy 𝐿 = 𝑟 ̇𝜃 + 𝑟 , 𝐻 = 12 ‖ 𝐯 ‖ + 0 . 𝑟 𝑞 = −1 , 𝑚 = 1 and starting from the initial position 𝐱 = (0 , −1 , 𝑇 with initial velocity 𝐯 = (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 refer-ence solution we take the output generated by the standard routine DOP853 based on a Runge–Kuttamethod 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 interms of the number of the computed sub-flows (by taking different time-steps). For clarity, hereand 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 includein all cases the splitting method (9) proposed in [12] (denoted here as ABC ) and the scheme S ,whose coefficients are collected in Table 2. .
25 3 .
50 3 .
75 4 .
00 4 .
25 4 .
50 4 . ( N eval ) − − − − − − − l og (cid:18) m a x ( | H − H | ) H (cid:19) ABC S XA (a) .
25 3 .
50 3 .
75 4 .
00 4 .
25 4 .
50 4 . ( N eval ) − − − − − − − l og (cid:18) m a x ( | H − H | ) H (cid:19) ABC S XB (b) .
25 3 .
50 3 .
75 4 .
00 4 .
25 4 .
50 4 . ( N eval ) − − − − − − − l og (cid:18) m a x ( | L − L | ) L (cid:19) ABC S XA (c) .
25 3 .
50 3 .
75 4 .
00 4 .
25 4 .
50 4 . ( N eval ) − − − − − − − − l og (cid:18) m a x ( | L − L | ) L (cid:19) ABC S XB (d) Figure 1:
Relative error in conserved quantities due to each of the best numerical methods tested for chargedparticle under Lorentz force. ( a ) Relative error in energy for XA compared to ABC and S . ( b ) Relative errorin energy for XB compared to ABC and S . ( c )Relative error in angular momentum for XA compared to ABC and S . ( d ) Relative error in angular momentum for XB compared to ABC and S . We notice that applying the composition methods proposed here leads to more accurate resultsthan the direct approach based on the splitting methods of Section 2 with the same computational cost,and that the new scheme XB is slightly more efficient that the the splitting scheme S (the remainingcomposition methods of Tables 2 and 3 provide results between ABC and the best compositiondepicted here).In Figure 2 we show the corresponding results obtained by each method for the error in the ( 𝐱 , 𝐯 ) space. 12 .
25 3 .
50 3 .
75 4 .
00 4 .
25 4 .
50 4 . ( N eval ) − − − − − − − − l og (cid:18) | ~ z − ~ z e x | | ~ z e x | (cid:19) ABC S XA (a) .
25 3 .
50 3 .
75 4 .
00 4 .
25 4 .
50 4 . ( N eval ) − − − − − − − − l og (cid:18) | ~ z − ~ z e x | | ~ z e x | (cid:19) ABC S XB XB (b) Figure 2:
Relative error in the ( 𝐱 , 𝐯 ) space for charged particle under Lorentz force. The notation is the sameas in Figure 1. ( a ) Relative error in the ( 𝐱 , 𝐯 ) space for XA compared to ABC and S . ( b ) Relative error inthe ( 𝐱 , 𝐯 ) space for XB and XB compared to ABC and S . One should notice that, although this system is Hamiltonian, the Hamiltonian function is notseparable into kinetic plus potential energy, and thus general symplectic Runge–Kutta methods cannotbe explicit [27]. In order to use explicit methods, one has to split the system into three parts. On theother hand, all the methods tested here are volume-preserving in the ( 𝐱 , 𝐯 ) space, just as the exactflow. The Hamiltonian of the disordered discrete nonlinear Schrödinger equation (DDNLS) = ∑ 𝑗 ( 𝜖 𝑗 | 𝜓 𝑗 | + 𝛽 | 𝜓 𝑗 | − ( 𝜓 𝑗 +1 𝜓 𝑗 + 𝜓 𝑗 +1 𝜓 𝑗 ) ) (35)describes a one-dimensional chain of couples nonlinear oscillators [7]. Here the sum extends over 𝑁 oscillators, 𝜓 𝑗 are complex variables, 𝛽 ≥ stands for the nonlinearity strength and the random ener-gies 𝜖 𝑗 are chosen uniformly from the interval [− 𝑊 ∕2 , 𝑊 ∕2] , where 𝑊 is related with the disorderstrength. 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 thananalyzing the rich dynamics this system possesses, our interest here is to use (35) as a non-trivial testbench for the integrators we presented in previous sections. By introducing the new (real) generalizedcoordinates and momenta ( 𝑞 𝑗 , 𝑝 𝑗 ) related with 𝜓 𝑗 through 𝜓 𝑗 = 1 √ 𝑞 𝑗 + 𝑖𝑝 𝑗 ) , 𝜓 𝑗 = 1 √ 𝑞 𝑗 − 𝑖𝑝 𝑗 ) , the Hamiltonian function (35) can be written as 𝐻 = 𝑁 ∑ 𝑗 =1 ( 𝜖 𝑗 𝑞 𝑗 + 𝑝 𝑗 ) + 𝛽 𝑞 𝑗 + 𝑝 𝑗 ) − 𝑝 𝑗 +1 𝑝 𝑗 − 𝑞 𝑗 +1 𝑞 𝑗 ) (36)13n such a way that is the sum of three explicitly solvable parts, 𝐻 = 𝐴 + 𝐵 + 𝐶 , with 𝐴 = 𝑁 ∑ 𝑗 =1 ( 𝜖 𝑗 𝑞 𝑗 + 𝑝 𝑗 ) + 𝛽 𝑞 𝑗 + 𝑝 𝑗 ) ) , 𝐵 = − 𝑁 ∑ 𝑗 =1 𝑝 𝑗 +1 𝑝 𝑗 , 𝐶 = − 𝑁 ∑ 𝑗 =1 𝑞 𝑗 +1 𝑞 𝑗 . The corresponding flows are given, respectively, by 𝜑 [ 𝑎 ] 𝑡 ∶ { 𝑞 𝑗 ( 𝑡 ) = 𝑞 𝑗 ( 𝑡 ) cos( 𝑎 𝑗 𝑡 ) + 𝑝 𝑗 ( 𝑡 ) sin( 𝑎 𝑗 𝑡 ) 𝑝 𝑗 ( 𝑡 ) = − 𝑞 𝑗 ( 𝑡 ) sin( 𝑎 𝑗 𝑡 ) + 𝑝 𝑗 ( 𝑡 ) cos( 𝑎 𝑗 𝑡 ) ,𝜑 [ 𝑏 ] 𝑡 ∶ { 𝑞 𝑗 ( 𝑡 ) = 𝑞 𝑗 ( 𝑡 ) − 𝑡 ( 𝑝 𝑗 −1 ( 𝑡 ) + 𝑝 𝑗 +1 ( 𝑡 )) 𝑝 𝑗 ( 𝑡 ) = 𝑝 𝑗 ( 𝑡 ) (37) 𝜑 [ 𝑐 ] 𝑡 ∶ { 𝑞 𝑗 ( 𝑡 ) = 𝑞 𝑗 ( 𝑡 ) 𝑝 𝑗 ( 𝑡 ) = 𝑝 𝑗 ( 𝑡 ) + 𝑡 ( 𝑞 𝑗 −1 ( 𝑡 ) + 𝑞 𝑗 +1 ( 𝑡 )) with 𝑎 𝑗 = 𝜖 𝑗 + 𝛽 ( 𝑞 𝑗 + 𝑝 𝑗 )∕2 . . . . . . . ( N eval ) − − − − − − l og (cid:18) m a x ( | H − H | ) H (cid:19) ABC S XA (a) . . . . . . ( N eval ) − − − − − − l og (cid:18) m a x ( | H − H | ) H (cid:19) ABC S XB (b) . . . . . . ( N eval ) − − − − − − l og (cid:18) m a x ( | S − S | ) S (cid:19) ABC S XA (c) . . . . . . ( N eval ) − − − − − − l og (cid:18) m a x ( | S − S | ) S (cid:19) ABC S XB (d) Figure 3:
Relative error in conserved quantities for the DDNLS system due to each of the best numericalmethods tested. ( a ) Relative error in energy for XA compared to ABC and S . ( b ) Relative error in energyfor XB compared to ABC and S . ( c ) Relative error in norm for XA compared to ABC and S . ( d ) Relativeerror in norm for XB compared to ABC and S . To compare the performance of the numerical integrators previously considered, we take a latticeof 𝑁 = 1000 sites and fixed boundary conditions, 𝑞 = 𝑝 = 𝑞 𝑁 +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 , and14he respective 𝑝 𝑖 in such a way that each site has the same constant norm 1, so that the total norm ofthe system is 𝑆 = 21 . Moreover, 𝛽 = 0 . , 𝑊 = 4 and the random disorder parameters 𝜖 𝑗 are chosenso that the total energy is 𝐻 ≈ −29 . . 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 andthe bottom to the error in norm. The same notation has been used for the tested methods. Finally, inFigure 4 we collect the error in the phase space. As before, the reference solution is obtained with theDOP853 routine. Notice that for this non trivial example the new schemes XA and especially XB show a better efficiency than S , and not only with respect to the preservation of the invariants, butalso in the computation of trajectories. . . . . . . ( N eval ) − − − − − − l og (cid:18) | ~ z − ~ z e x | | ~ z e x | (cid:19) ABC S XA (a) . . . . . . ( N eval ) − − − − − − l og (cid:18) | ~ z − ~ z e x | | ~ z e x | (cid:19) ABC S XB (b) Figure 4:
Relative error of trajectories for the DDNLS system. Same as in Figure 3. ( a ) Relative error oftrajectories for XA compared to ABC and S . ( b ) Relative error of trajectories for XB compared to ABC and S . In this work we have presented two different families of fourth-order composition methods especiallydesigned for problems that can be separated into three parts in such a way that each part is explicitlysolvable. In addition to the usual optimization criteria applied in the literature to choose the free pa-rameters in the composition, we have introduced another one especially oriented to problems wherethe 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 𝑆 of Table 2 onthe tested examples. Other relevant examples include certain nonlinear oscillators, Poisson–Maxwellequations arising in plasma physics, and the treatment of non separable Hamiltonian dynamical sys-tems [26].Although only problems separable into three parts have been considered here, it is clear thatthe schemes we have introduced can also be applied to differential equations split into any numberof pieces 𝑛 ≥ . The only modification one requires is to formulate the corresponding first orderscheme 𝜒 ℎ and its adjoint 𝜒 ∗ ℎ . One should be aware, however, that augmenting the number 𝑛 leadsto evaluating an increasingly large number of flows for methods with large values of 𝑠 , with thesubsequent 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 stabil-ity of splitting methods has been analyzed in particular in [31, 32], where highly efficient schemes15ith optimal stability polynomials have presented for numerically approximating the evolution of lin-ear problems. In the nonlinear case, however, the situation is more involved. In [19], a crude measureof the nonlinear stability of a given time symmetric scheme of order 𝑟 is proposed, taking into accountthe 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 appliedto linear reaction-diffusion equations with indefinite reaction terms which controls both low and highwave number instabilities. In any case, this result only affects methods up to order 2 with real andpositive coefficients, whereas the application of splitting and composition methods of higher orderwith real coefficients in this setting leads to severe instabilities due to the existence of negative coeffi-cients. 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 othersplitting 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 combinationof 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 strategyif necessary. Acknowledgements
This work has been funded by Ministerio de Economía, Industria y Competitividad (Spain) throughproject MTM2016-77660-P (AEI/FEDER, UE) and by Universitat Jaume I (projects UJI-B2019-17and GACUJI/2020/05). A.E.-T. has been additionally supported by the predoctoral contract BES-2017-079697 (Spain).
References [1] McLachlan, R.; Quispel, R. Splitting methods.
Acta Numer. , , 341–434.[2] Glowinski, R.; Osher, S.; Yin, W. (Eds.) Splitting Methods in Communication, Imaging, Sci-ence, and Engineering ; Springer: Berlin, Germany, 2016.[3] Hairer, E.; Lubich, C.; Wanner, G.
Geometric Numerical Integration. Structure-PreservingAlgorithms for Ordinary Differential Equations , 2nd ed.; Springer: Berlin, Germany, 2006.[4] Blanes, S.; Casas, F.
A Concise Introduction to Geometric Numerical Integration ; CRC Press:Boca Raton, FL, USA, 2016.[5] Sanz-Serna, J.; Calvo, M.
Numerical Hamiltonian Problems ; Chapman & Hall: London, UK,1994.[6] Leimkuhler, B.; Reich, S.
Simulating Hamiltonian Dynamics ; Cambridge University Press:Cambridge, UK, 2004.[7] Skokos, C.; Gerlach, E.; Bodyfelt, J.; Papamikos, G.; Eggl, S. High order three part split sym-plectic integrators: Efficient techniques for the long time simulation of the disordered discretenonlinear Schrödinger equation.
Phys. Lett. A , , 1809–1815.[8] Crouseilles, N.; Einkemmer, L.; Faou, E. Hamiltonian splitting for the Vlasov–Maxwell equa-tions. J. Comput. Phys. , , 224–240.[9] He, Y.; Sun, Y.; Liu, J.; Qin, H. Volume-preserving algorithms for charged particle dynamics. J. Comput. Phys. , , 135–147. 1610] Shang, X.; Kroger, M.; Leimkuhler, B. Assessing numerical methods for molecular and parti-cle simulation. Soft Matter , , 8565–8578.[11] Koseleff, P.V. Exhaustive search of symplectic integrators using computer algebra. In Inte-gration Algorithms and Classical Mechanics ; Marsden, J., Patrick, G., Shadwick, W., Eds.;American Mathematical Society: Providence, RI, USA, 1996.[12] Auzinger, W.; Hofstätter, H.; Ketcheson, D.; Koch, O. Practical splitting methods for the adap-tive integration of nonlinear evolution equations. Part I: Construction of optimized schemes andpairs of schemes.
BIT Numer. Math , , 55–74.[13] Arnold, V. Mathematical Methods of Classical Mechanics , 2nd ed.; Springer: Berlin, Ger-many, 1989.[14] Blanes, S.; Casas, F.; Murua, A. Splitting and composition methods in the numerical integra-tion of differential equations.
Bol. Soc. Esp. Mat. Apl. , , 89–145.[15] Yoshida, H. Construction of higher order symplectic integrators. Phys. Lett. A , , 262–268.[16] Diele, F.; Marangi, C. Explicit symplectic partitioned Runge–Kutta–Nyström methods fornon-autonomous dynamics. Appl. Numer. Math. , , 832–843.[17] Blanes, S.; Casas, F.; Murua, A. Composition methods for differential equations with process-ing. SIAM J. Sci. Comput. , , 1817–1843.[18] Suzuki, M. Fractal decomposition of exponential operators with applications to many-bodytheories and Monte Carlo simulations. Phys. Lett. A , , 319–323.[19] McLachlan, R. Families of High-Order Composition Methods. Numer. Algorithms , , 233–246.[20] Bernier, J.; Casas, F.; Crouseilles, N. Splitting methods for rotations: application to Vlasovequations. SIAM J. Sci. Comput. . accepted for publication.[21] Kahan, W.; Li, R. Composition constants for raising the order of unconventional schemes forordinary differential equations.
Math. Comput. , , 1089–1099.[22] Sofroniou, M.; Spaletta, G. Derivation of symmetric composition constants for symmetricintegrators. Optim. Method. Softw. , , 597–613.[23] Blanes, S.; Moan, P. Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nyströmmethods. J. Comput. Appl. Math. , , 313–330.[24] McLachlan, R. On the Numerical Integration of ODE’s by Symmetric Composition Methods. SIAM J. Sci. Comput. , , 151–168.[25] Blanes, S.; Casas, F.; Sanz-Serna, J. Numerical integrators for the Hybrid Monte Carlo method. SIAM J. Sci. Comput. , , A1556–A1580.[26] Tao, M. Explicit symplectic approximation of nonseparable Hamiltonians: algorithm and longtime performance. Phys. Rev. E , , 043303.[27] He, Y.; Sun, Y.; Liu, J.; Qin, H. Higher order volume-preserving schemes for charged particledynamics. J. Comput. Phys. , , 172–184.[28] Hairer, E.; Nørsett, S.; Wanner, G. Solving Ordinary Differential Equations I, Nonstiff Prob-lems , 2nd ed.; Springer: Berlin, Germany, 1993.1729] Kopidakis, G.; Komineas, S.; Flach, S.; Aubry, S. Absence of wave packet diffusion in disor-dered nonlinear systems.
Phys. Rev. Lett. , , 084103.[30] Danieli, C.; Manda, B.; Mithun, T.; Skokos, C. Computational efficiency of numerical inte-gration methods for the tangent dynamics of many-body Hamiltonian systems in one and twospatial dimensions. Math. Eng. , , 447–488.[31] McLachlan, R.; Gray, S. Optimal stability polynomials for splitting methods, with applicationsto the time-dependent Schrödinger equation. Appl. Numer. Math. , , 275–286.[32] Blanes, S.; Casas, F.; Murua, A. On the linear stability of splitting methods. Found. Comp.Math. , , 357–393.[33] Hundsdorfer, W.; Verwer, J. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations ; Springer: Berlin, Germany, 2003.[34] Ropp, D.; Shadid, J. Stability of operator splitting methods for systems with indefinite opera-tors: reaction-diffusion systems.
J. Comput. Phys. , , 449–466.[35] Blanes, S.; Casas, F.; Thalhammer, M. Splitting and composition methods with embeddederror estimators. Appl. Numer. Math. ,146