# 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.

## Authors

• 7 publications
• 4 publications
• ### Polynomial Neural Networks and Taylor maps for Dynamical Systems Simulation and Learning

The connection of Taylor maps and polynomial neural networks (PNN) to so...
12/19/2019 ∙ by Andrei Ivanov, et al. ∙ 0

• ### On symmetric-conjugate composition methods in the numerical integration of differential equations

We analyze composition methods with complex coefficients exhibiting the ...
01/11/2021 ∙ by Sergio Blanes, et al. ∙ 0

• ### Qualitative properties of numerical methods for the inhomogeneous geometric Brownian motion

We introduce the inhomogeneous geometric Brownian motion (IGBM) as a tes...
03/23/2020 ∙ by Irene Tubikanec, et al. ∙ 0

• ### On the Bayesian Solution of Differential Equations

The interpretation of numerical methods, such as finite difference metho...
05/18/2018 ∙ by Junyang Wang, et al. ∙ 0

• ### Gradient Flows and Accelerated Proximal Splitting Methods

Proximal based methods are well-suited to nonsmooth optimization problem...
08/02/2019 ∙ by Guilherme França, et al. ∙ 3

• ### Symplectic Geometric Methods for Matrix Differential Equations Arising from Inertial Navigation Problems

02/11/2020 ∙ by Xin-Long Luo, et al. ∙ 0

• ### Positivity-Preserving Adaptive Runge-Kutta Methods

Many important differential equations model quantities whose value must ...
05/13/2020 ∙ by Stephan Nüßlein, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Splitting methods are particularly useful for the numerical integration of ordinary differential equations (ODEs)

 ˙x≡dxdt=f(x),x(t0)=x0∈RD (1)

when the vector field

can be written as , so that each subproblem

 ˙x=fi(x),x(t0)=x0,i=1,…,n

is explicitly solvable, with solution . 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 general-purpose 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 used in GNI because the composition of symplectic (or volume preserving, orthogonal, etc.) transformations is again symplectic (volume preserving, orthogonal, etc., respectively). In composition methods the numerical scheme is constructed as the composition of several simpler integrators for the problem at hand, so as to improve their accuracy.

When in (1) can be separated into two parts, very efficient splitting schemes have been designed and applied to solve a wide variety of problems arising in several fields, ranging from Hamiltonian Monte 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 three or more parts in order to have subproblems that are explicitly solvable. Examples include the disordered 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 of accuracy can be built, the resulting algorithms involve such a large number of maps that they are not competitive in practice. It is the purpose of this paper to present an alternative class of efficient methods for the problem at hand and compare their performance on some non-trivial physical examples than can be split into three parts.

The paper is structured as follows. We first review how splitting methods can be directly applied to get numerical solutions (Section 2). Then the attention is turned to the application of composition methods, and we get a family of 4th-order schemes obtained by applying a standard optimization procedure (Section 3). In Section 4 we show how standard splitting methods, when formulated as a composition scheme, lead to very competitive integrators, and also propose a different optimization criterion for systems possessing invariant quantities. This allows us to get a new family of 4th-order schemes. All these integration algorithms are subsequently tested in Section 5 on a pair of numerical examples. Finally, Section 6 contains some concluding remarks.

## 2 First Approach: Splitting Methods

In what follows we assume that the vector field in (1) can be split into three parts,

 f(x)=fa(x)+fb(x)+fc(x) (2)

in such a way that the exact -flows , , , corresponding to , , , respectively, can be computed exactly.

It is clear that the composition

 χh=φ[a]h∘φ[b]h∘φ[c]h (3)

(or any other permutation of the sub-flows) provides a first-order approximation to the exact solution of (1), i.e.,

 χh(x0)=φh(x0)+O(h2),

whereas the so-called Strang splitting

 S[2]h=φ[a]h/2∘φ[b]h/2∘φ[c]h∘φ[b]h/2∘φ[a]h/2 (4)

Higher order approximations to the exact solution of (1) can be obtained by generalizing (4), i.e., by considering splitting schemes of the form

 ψ[r]h=φ[c]csh∘φ[b]bsh∘φ[a]ash∘⋯∘φ[c]c1h∘φ[b]b1h∘φ[a]a1h, (5)

where the coefficients , , are chosen to achieve a prescribed order of accuracy, say, ,

 ψ[r]h(x0)=φh(x0)+O(hr+1) as h→0. (6)

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 (i.e., for a consistency method) one has

 s∑i=1ai=1,s∑i=1bi=1,s∑i=1ci=1.

The specific number of order conditions is determined in fact by the dimension of the homogeneous subspace of grade , , of the free Lie algebra generated by the Lie derivatives corresponding to , , , respectively [1]. These dimensions are collected Table 1 for .

Thus, a splitting method (5) of order 4 requires solving order conditions and therefore the evaluation of at least a similar number of sub-flows to have as many parameters as equations. This number can be reduced by considering time-symmetric methods, i.e., schemes verifying

 ψ[r]h∘ψ[r]−h=id, (7)

where is the identity map. Condition (7) is verified by left-right palindromic compositions, i.e., if

 as+1−i=ai,bs+1−i=bi,cs+1−i=ci,i=1,2,…

in (5). Then all the conditions at even order are automatically satisfied. Thus, a symmetric method of order 4 requires solving 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:

 ψ[4]τ=φ[c]c1τ∘φ[b]b1τ∘φ[a]a1τ∘φ[b]b2τ∘φ[c]c2τ∘φ[b]b3τ∘φ[a]a2τ∘φ[b]b3τ∘φ[c]c2τ∘φ[b]b2τ∘φ[a]a1τ∘φ[b]b1τ∘φ[c]c1τ (8)

with

 a1=w1,a2=w0,b1=b2=w12,b3=w02,c1=w12,c2=w0+w12

and

 w1=12−21/3,w0=1−2w1.

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 . More recently, in [12] a method involving 21 maps of the form

 (a1b1c1a2b2c2a3b3c3a4b4a4c3b3a3c2b2a2c1b1a1) (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 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 and its adjoint,

 χ∗h:=(χ−h)−1=φ[c]h∘φ[b]h∘φ[a]h

with appropriately chosen weights. In other words, we look for integrators of the form

 ψh=χα2sh∘χ∗α2s−1h∘⋯∘χα2h∘χ∗α1h, with (α1,…,α2s)∈R2s (10)

verifying in addition the time-symmetry condition for all .

###### Remark 1

Methods of the form

 ψh=S[2]αmh∘S[2]αm−1h∘⋯∘S[2]α2h∘S[2]α1h with (α1,…,αm)∈Rm (11)

and (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 , then it is clear that when analyzing methods (10) we also recover schemes of the form (11).

### 3.1 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.

As is well known, for each infinitely differentiable map , the function admits an expansion of the form [13, 5]

 g(φh(x))=exp(hF)[g](x)=g(x)+∑k≥1hkk!Fk[g](x),x∈RD,

where is the Lie derivative associated with ,

 F≡Lf=D∑i=1fi(x)∂∂xi. (12)

Analogously, for the basic method one can associate a series of linear operators so that [14]

 g(χh(x))=exp(Y(h))[g](x), with Y(h)=∑k≥1hkYk

for all functions , whereas for its adjoint one has

 g(χ∗h(x))=exp(−Y(−h))[g](x).

Then the operator series associated with the integrator (10) is

 Ψ(h) = exp(−Y(−hα1))exp(Y(hα2))⋯exp(−Y(−hα2s−1))exp(Y(hα2s)).

Notice that the order of the operators is the reverse of the maps in (10) ([3] p. 88). Now, by repeated application of the Baker–Campbell–Hausdorff formula [4] we can express formally as the exponential of an operator ,

 Ψ(h)=exp(~F(h)), with ~F(h)=∑k≥1hkFk, (13)

for each and is the graded free Lie algebra generated by the operators , where, by consistency, . One has explicitly

 Y(hαi) = hαiY1+(hαi)2Y2+(hαi)3Y3+⋯ −Y(−hαi) = hαiY1−(hαi)2Y2+(hαi)3Y3−⋯

so that

 ~F(h) = hw1Y1+h2w2Y2+h3(w3Y3+w12[Y1,Y2]) +h4(w4Y4+w13[Y1,Y3]+w112[Y1,[Y1,Y2]]) +h5(w5Y5+w14[Y1,Y4]+w113[Y1,Y1,Y3] +w1112[Y1,Y1,Y1,Y2]+w23[Y2,Y3]+w212[Y2,Y1,Y2])+O(h6),

where , etc, refers to the usual Lie bracket and are polynomials in the coefficients . In particular, one has

 w1=2s∑i=1αi,w2=2s∑i=1(−1)iα2i, (15) w3=2s∑i=1α3i,w4=2s∑i=1(−1)iα4i, w12=12(2s−1∑i=1(−1)i+1α2i2s∑j=i+1αj+2s−1∑i=1αi2s∑j=i+1(−1)jα2j).

Thus, a time-symmetric 4th-order method has to satisfy only consistency () and the order conditions at order three, . Notice, then, that the minimum number of maps to be considered is . In that case the integrator reads

 ψh=χα1∘χ∗α2∘χα3∘χ∗α3∘χα2∘χ∗α1 (16)

and the unique (real) solution is given by

 α1=α2=12(2−21/3),α3=12−2α1.

This scheme corresponds to the familiar triple-jump integrator [15]

 ψh=S[2]αh/2∘S[2]βh∘S[2]αh/2 with α=1/(2−21/3). (17)

If , then involves 13 maps (the minimum number) and corresponds precisely to the splitting method (8).

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].

### 3.2 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.,

 E1(\boldmathα)=2s∑i=1|αi| (18)

is as small as possible, the logic being that higher order terms in the expansion (3.1) involve powers of these coefficients. In fact, methods with small values of 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 multiplying in the expansion (3.1), so that it makes sense to minimize

 (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 values of the objective functions for the triple-jump (17) are and , respectively.

Next we collect the most efficient schemes we have obtained with by applying this strategy.

#### s=4 stages.

The composition is

 ψh=χα1∘χ∗α2∘χα3∘χ∗α4∘χα4∘χ∗α3∘χα2∘χ∗α1, (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 , and the resulting coefficients are collected in Table 2 as method XA. In that case, and .

#### s=5 stages.

The resulting composition

 ψh=χα1∘χ∗α2∘χα3∘χ∗α4∘χα5∘χ∗α5∘χα4∘χ∗α3∘χα2∘χ∗α1

involves 21 maps when applied to a system separable into three parts. Minimum values for and are achieved when

 α1=α2=α3=α4=12(4−41/3),α5=12−4α1.

In consequence, the method can be written as

 ψh=S[2]αh∘S[2]αh∘S[2]βh∘S[2]αh∘S[2]αh

with , . Then and . This method, denoted XA, was first proposed in [18] and analyzed in detail in [19].

#### s=6 stages.

Analogously we have considered a composition involving three free parameters (and 25 maps when is given by (3)):

 ψh=χα1∘χ∗α2∘χα3∘χ∗α4∘χα5∘χ∗α6∘χα6∘χ∗α5∘χα4∘χ∗α3∘χα2∘χ∗α1. (21)

The proposed solution is collected in Table 2 as method XA leading to , . 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 numerical integration of the so-called reduced Vlasov–Maxwell system [20].

We could of course increase the number of stages. It turns out, however, that with 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 and one has the methods designed in [21] (37 maps) and [22] (45 maps), respectively, when the basic scheme is given by (3).

## 4 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) can be split into two parts, , which each part explicitly solvable, and take . Then, the adjoint method reads and the composition (10) adopts the form

 ψh=(φ[b]α2sh∘φ[a]α2sh)∘(φ[a]α2s−1h∘φ[b]α2s−1h)∘⋯∘(φ[b]α2h∘φ[a]α2h)∘(φ[a]α1h∘φ[b]α1h). (22)

Since , are exact flows, then they verify and (22) can be rewritten as the splitting scheme

 ψh=φ[b]bs+1h∘φ[a]ash∘φ[b]bsh∘⋯∘φ[b]b2h∘φ[a]a1h∘φ[b]b1h (23)

if and

 aj=α2j+α2j−1,bj+1=α2j+1+α2j,j=1,…,s (24)

(with ). Conversely, any integrator of the form (23) with can be expressed in the form (10) with and

 α2s=bs+1,α2j−1=aj−α2j,α2j−2=bj−α2j−1,j=s,s−1,…,1,

with for consistency. In consequence, any splitting method in principle designed for systems of the form 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, , 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 parts has been presented in ([23] Table 2) (method S) 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 and are and .

### 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

 ˙qi=∂H∂pi,˙pi=−∂H∂qi,i=1,…,d (25)

so that , and is the Hamiltonian. The Lie derivative associated with verifies, for any function ,

 LXHG=−{H,G}=−d∑j=1(∂H∂qj∂G∂pj−∂G∂qj∂H∂pj).

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 Hamiltonian in the limit for a 4th-order splitting method (23). A straightforward calculation shows that

 ~Hh=H+h4k5,1{A,A,A,A,B}+h4k5,2{B,A,A,A,B}+h4k5,3{A,A,B,A,B} (26) +h4k5,4{A,B,B,A,B}+h4k5,5{B,A,B,A,B}+h4k5,6{B,B,B,A,B}+ +h59∑j=1k6,jE6,j+O(h6),

where are polynomials in the coefficients , , refers to the iterated Poisson bracket , 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

 H(qi+1,pi+1)=exp(−hL~Hh)H(qi,pi)=H(qi,pi)−hL~HhH(qi,pi)+12h2L2~HhH(qi,pi)+⋯,

where .

An elementary calculation shows that

 H(qi+1,pi+1)−H(qi,pi)=h5(k51E61+(k51−k53)E62+(k52−k53)E63+k54E64 +(k52−13k53)E65+(k55−13k53)E66+(k55+k54)E67+(k56−k54)E68+k56E69) +O(h6).

Thus, for small ,

 Δ≡k251+(k51−k53)2+(k52−k53)2+k254+(k52−13k53)2 (27) +(k55−13k53)2+(k55+k54)2+(k56−k54)2+k256

can be taken as a measure of the energy error, and consequently,

 E3=2sΔ1/4 (28)

constitutes a possible objective function to minimize. The previous analysis can be also carried out for a composition method (10), resulting in

 Δ=w25+w214+w2113+w21112+w223+w2212. (29)

The -stage methods XB whose coefficients are collected in Table 3 have been obtained by minimizing 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.

## 5 Numerical Examples

Although optimization criteria based on the objective functions , and 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 in the sequel) and the 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 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 , where one has to take into account the symmetry of the coefficients: , etc. and :

 y=xndo j=1:6y=φ[a]α2j−1hyy=φ[b]α2j−1hy˜α=α2j−1+α2jy=φ[c]˜αhyy=φ[b]α2jhyy=φ[a]α2jhyend xn+1=y

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.

### 5.1 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

 m¨x=q(E+˙x×B), (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

 ˙x=v (31) ˙v=qmE+ωb×v

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]:

 ddt(xv) = (v0)+(0qmE(x))+(0ω(x)b(x)×v) (32) = f[a](x,v)+f[b](x,v)+f[c](x,v).

The corresponding flows with initial condition are given by

 φ[c]t:{x(t)=x0v(t)=exp(tω(x0)^b0)v0 (33)

where

is the skew-symmetric matrix

 ^b(x)=⎛⎜⎝0−b3(x)b2(x)b3(x)0−b1(x)−b2(x)b1(x)0⎞⎟⎠

associated with .

As in [9], we consider a static, non-uniform electromagnetic field

 E=−∇V=0.01r3(xex+y ey),B=∇×A=rez (34)

derived from the potentials

 V=0.01r,A=r23eθ

respectively, in cylindrical coordinates and with the appropriate normalization. Then, it can be shown that both the angular momentum and energy

 L=r2˙θ+r33,H=12∥v∥2+0.01r