CMlargesymbols”00 CMlargesymbols”01 Splitting methods for differential
equations–References

# Splitting methods for differential equations

Sergio Blanes
Instituto de Matemática Multidisciplinar
Universitat Politècnica de València
46022-Valencia  Spain
E-mail<EMAIL_ADDRESS>
Fernando Casas
Departament de Matemàtiques and IMAC  Universitat Jaume I
12071-Castellón  Spain
E-mail<EMAIL_ADDRESS>
Ander Murua
Konputazio Zientziak eta A.A. Saila  Informatika Fakultatea  EHU/UPV
Donostia/San Sebastián  Spain
E-mail<EMAIL_ADDRESS>
###### Abstract

This overview is devoted to splitting methods, a class of numerical
integrators intended for differential equations that can be subdivided into
different problems easier to solve than the original system. Closely connected
with this class of integrators are composition methods, in which one or
several low-order schemes are composed to construct higher-order numerical
approximations to the exact solution. We analyze in detail the order
conditions that have to be satisfied by these classes of methods to achieve a
given order, and provide some insight about their qualitative properties in
connection with geometric numerical integration and the treatment of highly
oscillatory problems. Since splitting methods have received considerable
attention in the realm of partial differential equations, we also cover this
subject in the present survey, with special attention to parabolic equations
and their problems. An exhaustive list of methods of different orders is
collected and tested on simple examples. Finally, some applications of
splitting methods in different areas, ranging from celestial mechanics to
statistics, are also provided.

###### doi:

XXXXXXXX

28th August 2024

###### Contents

  1. 1 Introduction
    1. 1.1 Lie–Trotter and Strang methods
    2. 1.2 Flows and differential operators
    3. 1.3 Adjoint method, conjugate method
    4. 1.4 The mathematical pendulum
    5. 1.5 The gravitational N-body problem
    6. 1.6 The time-dependent Schrödinger equation
    7. 1.7 Splitting methods as geometric numerical integrators
    8. 1.8 Relevance of splitting methods
    9. 1.9 Some historical remarks
    10. 1.10 Plan of the paper
  2. 2 High order splitting and composition methods
    1. 2.1 Raising the order by composition
      1. 2.1.1 Composition of Strang maps
      2. 2.1.2 Composition of Lie-Trotter maps
    2. 2.2 Order conditions I: splitting schemes with the BCH formula
    3. 2.3 Order conditions II: splitting schemes with Lyndon words
      1. 2.3.1 Basic expansions and necessary order conditions
      2. 2.3.2 Integral representation and remainders
      3. 2.3.3 Iterated integrals and shuffle relations
    4. 2.4 Order conditions III: splitting methods with Lyndon multi-indices
    5. 2.5 Order conditions IV: composition methods with Lyndon multi-indices
      1. 2.5.1 Order conditions of compositions of a basic method and its adjoint
      2. 2.5.2 Explicit characterization of the order conditions of scheme (48)
    6. 2.6 Negative time steps
  3. 3 Splitting methods for special problems
    1. 3.1 RKN splitting methods
    2. 3.2 Methods with commutators
    3. 3.3 Near-integrable systems
    4. 3.4 Splitting methods for linear systems
    5. 3.5 Splitting methods for non-autonomous systems
  4. 4 Qualitative properties of splitting methods
    1. 4.1 Changes of variables and differential equations on smooth manifolds
    2. 4.2 Stability
    3. 4.3 Modified equations
    4. 4.4 Modified equations and long-term behavior
    5. 4.5 Processing and long-term precision
  5. 5 Highly oscillatory problems
    1. 5.1 Expansion of the exact solution
    2. 5.2 Expansion of the discrete solution given by Strang splitting
    3. 5.3 Composition formulae
    4. 5.4 Expansions for arbitrary splitting methods
    5. 5.5 Modified ODE for the discrete flow furnished by a splitting method
    6. 5.6 Splitting methods with processing
    7. 5.7 A processed Strang scheme
    8. 5.8 More general assumptions
  6. 6 Splitting methods for PDEs
    1. 6.1 Splitting, LOD and ADI methods
    2. 6.2 IMEX methods
    3. 6.3 Schrödinger equations
    4. 6.4 Parabolic evolution equations
  7. 7 Splitting methods with complex coefficients
    1. 7.1 Compositions
    2. 7.2 Symmetric-conjugate methods and unitary problems
    3. 7.3 Projecting on the real axis
  8. 8 A collection of splitting methods
    1. 8.1 Symmetric compositions of time-symmetric 2nd-order schemes
    2. 8.2 Splitting into two parts / composition of a basic first-order method and its adjoint
    3. 8.3 Runge–Kutta–Nyström methods
    4. 8.4 Methods for near-integrable systems
  9. 9 Splitting everywhere: some relevant applications
    1. 9.1 Splitting and extrapolation
      1. 9.1.1 The Gragg/GBS method
      2. 9.1.2 Multi-product expansions
    2. 9.2 Crouch–Grossman methods
      1. 9.2.1 Positivity-preserving splitting methods
    3. 9.3 Exact splitting
    4. 9.4 Gravitational $N$-body problem
    5. 9.5 Molecular simulations
    6. 9.6 Hamiltonian Monte Carlo
    7. 9.7 Quantum statistical mechanics
    8. 9.8 Vlasov–Poisson equations
    9. 9.9 Quantum simulation of quantum systems
    10. 9.10 Other topics
  10. A Testing splitting methods on matrices

##  1 Introduction

###  1.1 Lie–Trotter and Strang methods

If, as has sometimes been argued, there are only ten big ideas in numerical
analysis and all the rest are merely variations on those themes, _splitting_
is undoubtedly one of them MacNamara and Strang (2016). In fact, one could say
that, at least since Descartes stated his four rules of logic in the _Discours
de la méthode_ 111 _Le second, de diviser chacune des difficultés que
j’examinarais, en autant de parcelles qu’il se pourrait, et qu’il serait
requis pour les mieux résoudre. Le troisième, de conduire par ordre mes
pensées, en commençant par les objets les plus simples et les plus aisés à
connaître, pour monter peu à peu, comme par degrés jusques à la connaissance
des plus composés; et supposant même de l’ordre entre ceux qui ne se précèdent
point naturellement les uns les autres._ (René Descartes, _Discours de la
méthode_ , Seconde partie)., the notion of sub-dividing a complicated problem
into its simpler constituent parts, solving each one of them separately and
combining those separated solutions in a controlled way to get a solution to
the original overall problem constitutes a guiding principle in all areas of
science and philosophy.

In the realm of numerical analysis of differential equations, this basic
principle can be stated as follows. Suppose one has the abstract initial value
problem

$x^{\prime}\equiv\frac{dx}{dt}=f(x),\qquad x(0)=x_{0}$ (1)

associated to an ordinary differential equation (ODE) or a partial
differential equation (PDE), in which case $f$ is a certain spatial partial
differential operator. Furthermore, suppose that $f$ does not depend
explicitly on time and can be decomposed as

$f=f_{1}+\cdots+f_{m},\qquad m\geq 2,$ (2)

so that each initial value problem $x^{\prime}=f_{j}(x)$, $x(0)=x_{0}$ is
easier to solve than (1). Most commonly, they can be integrated exactly in
closed form. Then, it is possible to take advantage of the decomposition (2)
to get accurate approximations to the solution of (1) by means of _splitting
methods_.

The best example to start with is perhaps the case of a linear differential
equation defined in $\mathbb{R}^{D}$ and $m=2$, namely

$x^{\prime}=f_{1}(x)+f_{2}(x)=F_{1}\,x+F_{2}\,x,\qquad x(0)=x_{0},$ (3)

where $F_{1}$ and $F_{2}$ are $D\times D$ matrices and $x\in\mathbb{R}^{D}$.
The solution reads

$x(t)=\mathrm{e}^{t(F_{1}+F_{2})}x_{0},$

so that by computing this matrix exponential directly we would have solved (3)
without requiring splitting methods. Associated with (3) one has the matrix
differential equation

$\frac{dX}{dt}=(F_{1}+F_{2})X,\qquad X(0)=I,$ (4)

in the sense that

$x(t)=X(t)\,x_{0}\qquad\mbox{ and }\qquad X(t)=\mathrm{e}^{t(F_{1}+F_{2})}.$

This is useful (mainly) for theoretical purposes, since one usually tries to
compute $\mathrm{e}^{t(F_{1}+F_{2})}x_{0}$ directly instead of first computing
the matrix exponential and then multiplying it by $x_{0}$.

It often happens that evaluating the action of $X(t)$ on $x_{0}$ is difficult
and/or computationally expensive. If, however, this is not the case for each
$\mathrm{e}^{tF_{j}}$ separately, then one may use the well known Lie product
formula (Reed and Simon 1980, p. 295)

$\mathrm{e}^{t(F_{1}+F_{2})}=\lim_{n\rightarrow\infty}\left(\mathrm{e}^{\frac{t}{n}F_{2}}\,\mathrm{e}^{\frac{t}{n}F_{1}}\right)^{n}.$
(5)

To get an approximate solution of (3) at the final time $t=t_{f}$, one
subdivides the interval $[0,t_{f}]$ into $N$ steps of length $h$, with
$Nh=t_{f}$, and computes the sequence

$x_{n+1}=\mathrm{e}^{hF_{2}}\,\mathrm{e}^{hF_{1}}\,x_{n},\qquad n\geq 0,$ (6)

so that $x_{n+1}\approx x(t_{n+1}=(n+1)h)$. This is the so-called _Lie–Trotter
scheme_. When the matrices commute, the sequence produces the exact solution.
To put it another way, if the commutator $[F_{1},F_{2}]\equiv
F_{1}F_{2}-F_{2}F_{1}=0$, then
$\exp(h(F_{1}+F_{2}))=\exp(hF_{2})\exp(hF_{1})$. Otherwise, a direct
calculation shows that

$\mathrm{e}^{h(F_{1}+F_{2})}-\mathrm{e}^{hF_{2}}\mathrm{e}^{hF_{1}}=\frac{1}{2}h^{2}[F_{1},F_{2}]+\mathcal{O}(h^{3})$

as $h\rightarrow 0$, and hence the previous approximation is only of first
order of accuracy. Another version of the method is possible of course by
reversing the order of $F_{1}$ and $F_{2}$, namely

$x_{n+1}=\mathrm{e}^{hF_{1}}\,\mathrm{e}^{hF_{2}}\,x_{n},\qquad n\geq 0$ (7)

has the same order of accuracy and properties as (6). Needless to say, for any
$m>2$, it results in

$x_{n+1}=\mathrm{e}^{hF_{m}}\,\cdots\,\mathrm{e}^{hF_{2}}\,\mathrm{e}^{hF_{1}}\,x_{n}$

(or any other permutation of the matrices $F_{j}$).

A higher order approximation can be achieved by considering a symmetrized
version of (6),

$\mathrm{e}^{h(F_{1}+F_{2})}-\mathrm{e}^{\frac{1}{2}hF_{1}}\mathrm{e}^{hF_{2}}\mathrm{e}^{\frac{1}{2}hF_{1}}=Ch^{3}+\mathcal{O}(h^{4}),$
(8)

where the constant $C$ can be obtained either by comparing Taylor series or by
applying the Baker–Campbell–Hausdorff (BCH) formula as Varadarajan (1984)
$C=\frac{1}{24}([F_{1},[F_{1},F_{2}]]+2[F_{2},[F_{1},F_{2}]])$. Therefore, the
sequence

$x_{n+1}=\mathrm{e}^{\frac{1}{2}hF_{1}}\mathrm{e}^{hF_{2}}\mathrm{e}^{\frac{1}{2}hF_{1}}\,x_{n},\qquad
n\geq 0$ (9)

produces a second-order approximation for the solution of (3). This
corresponds to the _Strang_ (splitting) _scheme_. Again, if the role of
$F_{1}$ and $F_{2}$ is interchanged, one has another version of the Strang
splitting scheme.

Simple generalizations to the case $m>2$ include in particular the product

$\mathrm{e}^{\frac{1}{2}hF_{1}}\,\mathrm{e}^{\frac{1}{2}hF_{2}}\,\cdots\,\mathrm{e}^{hF_{m}}\,\cdots\,\mathrm{e}^{\frac{1}{2}hF_{2}}\,\mathrm{e}^{\frac{1}{2}hF_{1}}.$

Higher order splitting methods could in principle be constructed by including
more exponentials with their corresponding coefficients in a time step, namely

$\Psi(h)=\mathrm{e}^{a_{s+1}hF_{1}}\,\mathrm{e}^{b_{s}hF_{2}}\,\mathrm{e}^{a_{s}hF_{1}}\,\cdots\,\mathrm{e}^{b_{1}hF_{2}}\,\mathrm{e}^{a_{1}hF_{1}}.$
(10)

The number $s$ as well as the coefficients $a_{j}$, $b_{j}$ are chosen so that

$\Psi(h)=\mathrm{e}^{h(F_{1}+F_{2})}+\mathcal{O}(h^{r+1})$

as $h\rightarrow 0$ for a given order $r$. In (10), the first and last
exponentials correspond to $F_{1}$. This format is convenient for
implementation, since the last exponential in one step can be concatenated
with the first one at the next step, thus reducing the number of evaluations
by one. This corresponds to the well known FSAL (First Same As Last) property.
The situation when one has an exponential of $F_{2}$ as the first and last
term is recovered by taking $a_{1}=a_{s+1}=0$.

###  1.2 Flows and differential operators

The Lie–Trotter scheme can be easily generalized to any system (1)-(2) when
the solution is no longer given by exponentials, as in the linear case. If
$m=2$, it is equivalent to the following.

###### Algorithm 1.1

(Lie–Trotter). Starting from $x_{0}=x(0)$, for $n\geq 0$

  * •

solve $y_{1}^{\prime}=f_{1}(y_{1})$, $\ y_{1}(t_{n})=x_{n}$, in
$[t_{n},t_{n+1}]$;

  * •

set $y_{n+1/2}=y_{1}(t_{n+1})$;

  * •

solve $y_{2}^{\prime}=f_{2}(y_{2})$, $\ y_{2}(t_{n})=y_{n+1/2}$, in
$[t_{n},t_{n+1}]$;

  * •

finally, set $x_{n+1}=y_{2}(t_{n+1})$.

Alternatively, if we denote the solution of equation (1) for each
$t\in\mathbb{R}$ as

$x(t)=\varphi_{t}^{[f]}(x_{0}),$

then Algorithm 1.1 can be formally expressed as

$x_{n+1}=\chi_{h}(x_{n})\equiv\varphi_{h}^{[2]}\big{(}\varphi_{h}^{[1]}(x_{n})\big{)}=\big{(}\varphi_{h}^{[2]}\circ\varphi_{h}^{[1]}\big{)}(x_{n}),$
(11)

where we have used the simplified notation $\varphi_{t}^{[j]}(x_{0})$ for the
solutions $y(t)=\varphi_{t}^{[f_{j}]}(x_{0})$ of the subproblems
$y^{\prime}=f_{j}(y)$, $y(0)=x_{0}$. Analogously, the Strang splitting (9) is
generalized as follows.

###### Algorithm 1.2

(Strang). From $x_{0}=x(0)$, for $n\geq 0$,

  * •

solve $y_{1}^{\prime}=f_{1}(y_{1})$, $\ y_{1}(t_{n})=x_{n}$, in
$[t_{n},t_{n+\frac{1}{2}}]$, with $t_{n+\frac{1}{2}}=(n+\frac{1}{2})h$;

  * •

set $y_{n+\frac{1}{2}}=y_{1}(t_{n+\frac{1}{2}})$;

  * •

solve $y_{2}^{\prime}=f_{2}(y_{2})$, $\ y_{2}(t_{n})=y_{n+\frac{1}{2}}$, in
$[t_{n},t_{n+1}]$;

  * •

set $\hat{y}_{n+\frac{1}{2}}=y_{2}(t_{n+1})$;

  * •

solve $y_{1}^{\prime}=f_{1}(y_{1})$, $\
y_{1}(t_{n+\frac{1}{2}})=\hat{y}_{n+\frac{1}{2}}$, in
$[t_{n+\frac{1}{2}},t_{n+1}]$;

  * •

finally, set $x_{n+1}=y_{1}(t_{n+1})$

or, in short,

$x_{n+1}=S_{h}^{[2]}(x_{n})\equiv\big{(}\varphi_{h/2}^{[1]}\circ\varphi_{h}^{[2]}\circ\varphi_{h/2}^{[1]}\big{)}(x_{n}).$
(12)

If equation (1) corresponds to a (nonlinear) ordinary differential equation
evolving in $\mathbb{R}^{D}$,

$x^{\prime}=f(x),\qquad x(0)=x_{0}\in\mathbb{R}^{D},$ (13)

$f$ is called the _vector field_. If (13) admits for each
$x_{0}\in\mathbb{R}^{D}$ a unique solution $x(t)$ defined for all
$t\in\mathbb{R}$, the map

$\displaystyle\varphi_{t}^{[f]}:$
$\displaystyle\mathbb{R}^{D}\longrightarrow\mathbb{R}^{D}$ $\displaystyle
x_{0}\longmapsto x(t)=\varphi_{t}^{[f]}(x_{0})$

is referred to as the $t$-_flow_ Arnold (1989). Thus, for each value of the
real parameter $t$, $\varphi_{t}^{[f]}$ maps $\mathbb{R}^{D}$ in
$\mathbb{R}^{D}$ in such a way that $\varphi_{t}^{[f]}(z)$ is the value at
time $t$ of the solution of the system with initial value $z$ at time $0$,
whereas, for fixed $x_{0}$ and varying $t$, $\varphi_{t}^{[f]}(x_{0})$ is the
solution of the initial value problem (13).

It is worth mentioning that the solution $x(t)$ of (1) is in general defined
for a maximal time-interval $(t_{\min}(x_{0}),t_{\max}(x_{0}))$ (with
$-\infty\leq t_{\min}(x_{0})\leq 0\leq t_{\max}(x_{0})\leq+\infty$).
Furthermore, the vector field $f$ of many systems of ordinary differential
equations is singular (or undefined) for some $x\in\mathbb{R}^{D}$. Thus, in
general, $f$ is defined for some open set $\mathcal{U}\subset\mathbb{R}^{D}$.
Hence, for a given $t\in\mathbb{R}$, the $t$-flow $\varphi_{t}^{[f]}$ is a map
from

$\mathcal{D}_{t}=\\{x_{0}\in\mathcal{U}\ :\
t\in(t_{\min}(x_{0}),t_{\max}(x_{0}))\\}$

to $\mathcal{U}$. In general, one may have different domains of definition
$\mathcal{U}$ (resp. $\mathcal{D}_{t}$) for each $f_{j}$ (resp.
$\varphi_{t}^{[f_{j}]}$). In this general situation the compositions in (11)
and (12) are not well defined for all $x_{0}$, $h$ and $n$. In order to avoid
these technicalities, we will assume in what follows that
$\mathcal{U}=\mathbb{R}^{D}$, $t_{\min}=-\infty$, and $t_{\max}=+\infty$, for
each vector field $f_{j}$.

Associated with the vector field $f$ is the Lie derivative or _Lie operator_
$F$ Arnold (1989), mapping smooth functions
$g:\mathbb{R}^{D}\rightarrow\mathbb{R}$ into the real-valued function
$F\,g:\mathbb{R}^{D}\rightarrow\mathbb{R}$ such that, for
$x\in\mathbb{R}^{D}$,

$(F\,g)(x)=\left.\frac{d}{dt}\right|_{t=0}g\big{(}\varphi_{t}^{[f]}(x)\big{)},$

that is,

$(F\,g)(x)=f(x)\cdot\nabla g(x).$

Then, the flow of (13) verifies Sanz-Serna and Calvo (1994), Hairer et al.
(2006)

$g\big{(}\varphi_{t}^{[f]}(x)\big{)}=\big{(}\mathrm{e}^{tF}g\big{)}(x)\equiv\left(\sum_{k\geq
0}\frac{t^{k}}{k!}F^{k}g\right)(x).$

The operator $X(t)\equiv\mathrm{e}^{tF}$ is called _Lie transformation_ , and
can be seen as the formal solution of the operator equation

$\frac{dX}{dt}=X\,F,\qquad X(0)=I.$ (14)

This can be seen as follows: on the one hand,

$\frac{d}{dt}g\big{(}\varphi_{t}^{[f]}(x)\big{)}=\frac{d}{dt}\big{(}X(t)g\big{)}(x)=\left(\frac{dX(t)}{dt}g\right)(x),$

and on the other hand

$\frac{d}{dt}g\big{(}\varphi_{t}^{[f]}(x)\big{)}=(Fg)\big{(}\varphi_{t}^{[f]}(x)\big{)}=X(t)(Fg)(x).$

Lie operators satisfy some remarkable properties Arnold (1989). In particular,
although they do not commute, their commutator is nevertheless a first-order
linear differential operator. Specifically, let $F$ and $G$ be the Lie
operators associated with $f$ and $g$, respectively, and
$u:\mathbb{R}^{D}\rightarrow\mathbb{R}$ a given smooth function. Then,

$[F,G]u=(FG-GF)\,u=\sum_{i,j=1}^{D}\left(f_{j}\frac{\partial g_{i}}{\partial
x_{j}}-g_{j}\frac{\partial f_{i}}{\partial x_{j}}\right)\frac{\partial
u}{\partial x_{i}}$

and it is possible to associate a new vector field to this differential
operator, $w=(f,g)$, with components

$w_{i}=(f,g)_{i}=\sum_{j=1}^{D}\left(f_{j}\frac{\partial g_{i}}{\partial
x_{j}}-g_{j}\frac{\partial f_{i}}{\partial x_{j}}\right).$

It is called the _Lie–Poisson bracket_ of $f$ and $g$, and its Lie operator
$W$ satisfies $W=[F,G]$.

Now suppose that $f(x)=f_{1}(x)+f_{2}(x)$, so that each part
$x^{\prime}=f_{j}(x)$ is exactly solvable (or can be numerically solved up to
round-off accuracy) with flow $x(t)=\varphi^{[j]}_{t}(x_{0})$. Denoting by
$F_{1}$ and $F_{2}$ the Lie operators associated with $f_{1}$ and $f_{2}$,
respectively, it holds that

$g\big{(}\varphi^{[1]}_{t}(x)\big{)}=\left(\mathrm{e}^{tF_{1}}g\right)(x),\qquad
g\big{(}\varphi^{[2]}_{t}(x)\big{)}=\left(\mathrm{e}^{tF_{2}}g\right)(x).$

Then, for the first-order approximation
$\chi_{h}=\varphi^{[2]}_{h}\circ\varphi^{[1]}_{h}$ furnished by Algorithm 1.1
one has $g\big{(}\chi_{h}(x)\big{)}=\big{(}\Psi(h)g\big{)}(x)$, where
$\Psi(h)$ is a series of linear differential operators defined as

$\Psi(h)=\mathrm{e}^{hF_{1}}\,\mathrm{e}^{hF_{2}}.$ (15)

Notice that the exponentials of Lie derivatives in (15) appear in reverse
order with respect to the maps in the integrator (Hairer et al. 2006, p. 88).
Of course, the same procedure can be applied to the Strang splitting,
resulting in the product

$\Psi(h)=\mathrm{e}^{\frac{h}{2}F_{1}}\,\mathrm{e}^{hF_{2}}\,\mathrm{e}^{\frac{h}{2}F_{1}}.$
(16)

These considerations show that: (i) splitting methods for the problem (1)-(2)
defined in a certain function space (as for instance with partial differential
equations) can be formulated in terms of the solution of each subproblem
(either exact or approximate) by means of Algorithms 1.1 and 1.2, and (ii)
splitting methods applied to nonlinear ODEs evolving in $\mathbb{R}^{D}$ can
also be formally expressed as products of exponentials of differential
operators, since it is possible to transform the original nonlinear problem
into a linear one with the Lie formalism. This observation is very useful when
analyzing the order conditions for a method to be of a given order. In
particular, one has the same order conditions for linear and nonlinear ODEs
(see Section 2).

The previous integrators are sometimes called _multiplicative operator-
splitting methods_ , especially in the literature concerning the numerical
treatment of partial differential equations. In that area, one has still to
specify how to solve each initial value sub-problem in Algorithms 1.1-1.2 as
well as the boundary conditions. Moreover, one should take into account that,
for a given differential equation, different ways to carry out the splitting
in fact lead to different integrators.

###  1.3 Adjoint method, conjugate method

The flow $\varphi_{t}^{[f]}$ of (13) verifies
$\big{(}\varphi_{-t}^{[f]}\big{)}^{-1}=\varphi_{t}^{[f]}$, but this property
is not shared by many numerical integrators, and in particular by the map
$\chi_{h}$ corresponding to the Lie–Trotter scheme.

In general, if $\psi_{h}(x)$ represents a numerical method of order at least
one, i.e., $\psi_{h}(x)=x+h\,f(x)+\mathcal{O}(h^{2})$, then
$(\psi_{-h})^{-1}(x)=x+h\,f(x)+\mathcal{O}(h^{2})$, so that

$\psi_{h}^{*}\equiv(\psi_{-h})^{-1}$

is also a numerical method of order at least one. It is called the adjoint
method of $\psi_{h}$ Sanz-Serna and Calvo (1994). In words, stepping forwards
with the given method $\psi_{h}$ is the same as stepping backwards with the
inverse of its adjoint $\psi_{h}^{*}$. If
$\psi_{h}\equiv\chi_{h}=\varphi_{h}^{[2]}\circ\varphi_{h}^{[1]}$, then clearly
$\chi_{h}^{*}=\varphi_{h}^{[1]}\circ\varphi_{h}^{[2]}$. Additional examples
are the explicit and implicit Euler methods:

$x_{n+1}=\psi_{h}^{e}(x_{n})=x_{n}+hf(x_{n}),\qquad
x_{n+1}=\psi_{h}^{i}(x_{n})=x_{n}+hf(x_{n+1}),$

since $\psi_{h}^{i}=(\psi_{h}^{e})^{*}$.

Whenever an integrator satisfies

$\psi_{h}=\psi_{h}^{*}=(\psi_{-h})^{-1},$

it is called a time-symmetric or self-adjoint method. Alternatively,
$x_{n+1}=\psi_{h}(x_{n})$ is time-symmetric if and only if exchanging
$h\leftrightarrow-h$ and $x_{n}\leftrightarrow x_{n+1}$ one gets the same
expression, i.e., $\psi_{-h}(x_{n+1})=x_{n}$. The Strang scheme (12) is an
example of a time-symmetric method.

It is in fact straightforward to construct time-symmetric methods using the
adjoint: given an arbitrary method $\psi_{h}$ of order $r\geq 1$, then the
compositions

$\psi_{h/2}\circ\psi_{h/2}^{*}\quad\mbox{ and
}\quad\psi_{h/2}^{*}\circ\psi_{h/2}$ (17)

are time-symmetric methods of order $r\geq 2$ Sanz-Serna and Calvo (1994).
Moreover, symmetric methods are necessarily of even order, as we will show in
Section 2. Notice that the Strang method (12) is simply

$S_{h}^{[2]}=\chi_{h/2}^{*}\circ\chi_{h/2},$

where $\chi_{h}$ is given by eq. (11). Additional examples are the trapezoidal
rule $\psi_{h}^{t}=\psi_{h/2}^{i}\circ\psi_{h/2}^{e}$ and the midpoint rule
$\psi_{h}^{m}=\psi_{h/2}^{e}\circ\psi_{h/2}^{i}$.

The Strang scheme can also be expressed as

$\displaystyle S_{h}^{[2]}$
$\displaystyle=\varphi_{-h/2}^{[1]}\circ(\varphi_{h}^{[1]}\circ\varphi_{h}^{[2]})\circ\varphi_{h/2}^{[1]}$
(18) $\displaystyle=\pi_{h}^{-1}\circ\chi_{h}^{*}\circ\pi_{h},$

with $\pi_{h}=\varphi_{h/2}^{[1]}$. In the terminology of dynamical systems,
the Strang and Lie–Trotter schemes are said to be conjugate to each other by
the ($\mathcal{O}(h)$-near to identity) map $\pi_{h}=\varphi_{h/2}^{[1]}$,
which can be considered as a change of coordinates. Furthermore, the result of
$n$ applications of the Strang scheme can be recovered from $n$ applications
of Lie–Trotter by carrying out just an initial transformation at the initial
step and its _inverse_ at the final step.

Since many dynamical properties are invariant under changes of coordinates,
conjugate methods provide the same characterization of these properties. Other
examples of conjugate methods are the trapezoidal and midpoint rules, the map
$\pi_{h}$ being in this case the implicit Euler method Hairer et al. (2006).
We will treat conjugate methods in detail in Sections 4 and 5.

###  1.4 The mathematical pendulum

The simple mathematical pendulum constitutes a standard example of a nonlinear
Hamiltonian system. In appropriate units, the corresponding Hamiltonian
function reads

$H(q,p)=\frac{1}{2}p^{2}+(1-\cos q),$ (19)

where $q$ denotes the angle from the vertical suspension point and $p$ is the
associated momentum.

As is well known, the equations of motion of a generic Hamiltonian system with
Hamiltonian $H(q,p)$, and $q,p\in\mathbb{R}^{d}$, are given by Goldstein
(1980)

$\frac{dq}{dt}=\nabla_{p}H,\qquad\quad\frac{dp}{dt}=-\nabla_{q}H,$ (20)

the function $H(q,p)$ remains constant along the evolution, and the
corresponding $t$-flow, denoted as $\varphi_{t}^{[H]}$, is a _symplectic_
transformation Arnold (1989): its Jacobian matrix $\varphi_{t}^{\prime[H]}$
verifies the identity

$(\varphi_{t}^{\prime[H]})^{T}\,J\,\varphi_{t}^{\prime[H]}=J\qquad\mbox{ for
}\;t\geq 0,$

where $J$ is the basic canonical matrix

$J=\left(\begin{array}[]{rr}0_{d}&I_{d}\\\ -I_{d}&0_{d}\end{array}\right).$
(21)

In the particular case of (19), $d=1$, and the equations of motion are

$\frac{dq}{dt}=p,\qquad\qquad\frac{dp}{dt}=-\sin q.$

Given a Hamiltonian $H(q,p)$ that can be decomposed as

$H(q,p)=H_{1}(q,p)+H_{2}(q,p),$ (22)

it makes sense to split the equations of motion (20) as

$\frac{d}{dt}\left(\begin{matrix}q\\\
p\end{matrix}\right)=\left(\begin{matrix}\nabla_{p}H_{1}(q,p)\\\
-\nabla_{q}H_{1}(q,p)\end{matrix}\right)+\left(\begin{matrix}\nabla_{p}H_{2}(q,p)\\\
-\nabla_{q}H_{2}(q,p)\end{matrix}\right),$

so that each subsystem is itself Hamiltonian. In that case, we can then apply
Algorithm 1.1 by composing both maps and form the first-order scheme

$x_{n+1}=\chi_{h}(x_{n})\equiv\big{(}\varphi_{h}^{[H_{2}]}\circ\varphi_{h}^{[H_{1}]}\big{)}(x_{n}),\quad
n=0,1,2,\ldots,$ (23)

where $x_{n}=(q_{n},p_{n})^{T}$. Similarly, Algorithm 1.2 gives the second-
order scheme

$x_{n+1}=S_{h}^{[2]}(x_{n})\equiv\big{(}\varphi_{h/2}^{[H_{1}]}\circ\varphi_{h}^{[H_{2}]}\circ\varphi_{h/2}^{[H_{1}]}\big{)}(x_{n}),\quad
n=0,1,2,\ldots.$ (24)

Notice that, since both $\chi_{h}$ and $S_{h}^{[2]}$ are defined as
compositions of flows of Hamiltonian systems and the composition of symplectic
maps is also symplectic Arnold (1989), then, both (23) and (24) are
_symplectic integrators_ Sanz-Serna and Calvo (1994).

The fact that schemes (23) and (24) share the symplectic property with the
exact flow has remarkable consequences in practice concerning the preservation
of properties and the error propagation on long time integrations, as we will
shortly illustrate.

For the Hamiltonian (19) describing the pendulum (and in fact, for many other
mechanical systems), one can separate the contributions of the kinetic energy
$T(p)=\frac{1}{2}p^{2}$ and the potential energy $V(q)=1-\cos q$, so that a
natural splitting of the form (22) is then

$H(q,p)=T(p)+V(q).$ (25)

This corresponds to splitting the equations of motion (20) into the subsystems

$\left(\begin{array}[]{c}q^{\prime}\\\
p^{\prime}\end{array}\right)=\left(\begin{array}[]{c}\nabla_{p}T(p)\\\
0\end{array}\right)\qquad\quad\mbox{ and
}\qquad\quad\left(\begin{array}[]{c}q^{\prime}\\\
p^{\prime}\end{array}\right)=\left(\begin{array}[]{c}0\\\
-\nabla_{q}V(q)\end{array}\right),$ (26)

which in turn implies that

$\varphi_{t}^{[T]}:\;\left(\begin{matrix}q_{0}\\\
p_{0}\end{matrix}\right)\longmapsto\left(\begin{matrix}q_{0}+t\,\nabla
T_{p}(p_{0})\\\ p_{0}\end{matrix}\right)$ (27)

and

$\varphi_{t}^{[V]}:\;\left(\begin{matrix}q_{0}\\\
p_{0}\end{matrix}\right)\longmapsto\left(\begin{matrix}q_{0}\\\
p_{0}-t\,\nabla V_{q}(q_{0})\end{matrix}\right).$ (28)

Then, the first order scheme (23) reduces to

$\displaystyle q_{n+1}$ $\displaystyle=q_{n}+h\,\nabla_{p}T(p_{n}),\quad
p_{n+1}=p_{n}-h\,\nabla_{q}V(q_{n+1}),\qquad n=0,1,2\ldots$ (29)

Compared to the explicit Euler method

$q_{n+1}=q_{n}+h\,\nabla_{p}T(p_{n}),\qquad
p_{n+1}=p_{n}-h\,\nabla_{q}V(q_{n}),$

it only differs in that $\nabla_{q}V$ is evaluated at the updated value
$q_{n+1}$ instead of $q_{n}$. It makes sense, then, to call scheme (29) the
symplectic Euler-VT method: one first computes the gradient of the kinetic
energy $T$ and then computes the gradient of the potential energy $V$.

In accordance with our treatment in subsection 1.3, the adjoint of (29)
corresponds to composing the maps $\varphi_{t}^{[T]}$ and $\varphi_{t}^{[V]}$
in reverse order,

$\displaystyle p_{n+1}$ $\displaystyle=p_{n}-h\,\nabla_{q}V(q_{n}),\qquad
q_{n+1}=q_{n}+h\,\nabla T_{p}(p_{n+1}),$ (30)

so that we call it the symplectic Euler-TV method.

Obviously, our discussion of schemes (29) and (30) above applies for any
Hamiltonian system whose Hamiltonian function can be written in the (so-called
separable) form (25). That is, the Lie–Trotter scheme leads to the two
variants of the symplectic Euler method when it is applied to separable
Hamiltonian systems.

As for the Strang splitting scheme, described in Algorithm 1.2 in general, and
in (24) for Hamiltonian problems, when $H(q,p)=T(p)+V(q)$ it reduces to the
much celebrated _Störmer–Verlet method_ Hairer et al. (2003). Specifically,
depending on the order in which both parts are evaluated, one has the
following two variants:

###### Algorithm 1.3

(Störmer–Verlet-VTV). From $(q_{0},p_{0})=(q(0),p(0))$, for $n\geq 0$,

  * •

$p_{n+1/2}=p_{n}-\frac{h}{2}\nabla_{q}V(q_{n})$

  * •

$q_{n+1}=q_{n}+h\nabla_{p}T(p_{n+1/2})$

  * •

$p_{n+1}=p_{n+1/2}-\frac{h}{2}\nabla_{q}V(q_{n+1})$,

and

###### Algorithm 1.4

(Störmer–Verlet-TVT). From $(q_{0},p_{0})=(q(0),p(0))$, for $n\geq 0$,

  * •

$q_{n+1/2}=q_{n}+\frac{h}{2}\nabla_{p}T(p_{n})$

  * •

$p_{n+1}=p_{n}-h\nabla_{q}V(q_{n+1/2})$

  * •

$q_{n+1}=q_{n+1/2}+\frac{h}{2}\nabla_{p}T(p_{n+1})$,

respectively. Clearly, Algorithms 1.3 and 1.4 correspond to _time-symmetric_
methods and can be obtained by composing the Euler-TV method and its adjoint.
Specifically, if $\chi_{h}$ corresponds to method (30), then

$S_{h}^{[2]}\equiv\chi_{h/2}^{*}\circ\chi_{h/2}=\varphi_{h/2}^{[V]}\circ\varphi_{h}^{[T]}\circ\varphi_{h/2}^{[V]},$
(31)

recovers Störmer–Verlet-VTV, whereas the TVT version corresponds to
$\varphi_{h/2}^{[T]}\circ\varphi_{h}^{[V]}\circ\varphi_{h/2}^{[T]}$.

The top panel of Figure 1.1 shows trajectories of the pendulum (19), starting
from three different initial conditions:
$(q_{0},p_{0})=(-5,\frac{5}{2}),(1,1),(\frac{1}{10},0)$, corresponding to
different regions of the phase space: for $(q_{0},p_{0})=(-5,\frac{5}{2})$ it
holds that $T(p(t))>V(q(t))$, whereas for $(q_{0},p_{0})=(1,1)$ one has
$T(p(t))\simeq V(q(t))$ on average. Finally, for
$(q_{0},p_{0})=(\frac{1}{10},0)$, the system can be seen as a slightly
perturbed harmonic oscillator. In this case one could consider the following
decomposition:

$H=\frac{1}{2}(p^{2}+q^{2})+(1-\frac{1}{2}q^{2}-\cos q)\equiv
H_{1}(q,p)+H_{2}(q),$ (32)

where $H_{1}(q,p)=\frac{1}{2}(p^{2}+q^{2})$ corresponds to the harmonic
oscillator, whose exact solution is known (a rotation in phase space).
Moreover, $|H_{2}(q)|=\varepsilon|H_{1}(q,p)|$, with $\varepsilon\approx
10^{-3}$ along the orbit originated in $(\frac{1}{10},0)$. With splitting
(32), the map $\varphi_{h}^{[H_{1}]}\circ\varphi_{h}^{[H_{2}]}$ reads

$\left(\begin{array}[]{c}q_{n+1}\\\
p_{n+1}\end{array}\right)=R(h)\left(\begin{array}[]{c}q_{n}\\\
p_{n}+h(q_{n}-\sin q_{n})\end{array}\right),\quad\mbox{with}\quad
R(h)=\left(\begin{array}[]{cc}\cos h&\sin h\\\ -\sin h&\cos
h\end{array}\right),$ (33)

whereas the second-order scheme
$\varphi_{h/2}^{[H_{1}]}\circ\varphi_{h}^{[H_{2}]}\circ\varphi_{h/2}^{[H_{1}]}$
can be formulated as

$\displaystyle\left(\begin{array}[]{c}q_{n+1/2}\\\
p_{n+1/2}\end{array}\right)=R(h/2)\left(\begin{array}[]{c}q_{n}\\\
p_{n}\end{array}\right)$ (34) $\displaystyle
p_{n+1/2}^{*}=p_{n+1/2}+h(q_{n+1/2}-\sin q_{n+1/2})$
$\displaystyle\left(\begin{array}[]{c}q_{n+1}\\\
p_{n+1}\end{array}\right)=R(h/2)\left(\begin{array}[]{c}q_{n+1/2}\\\
p_{n+1/2}^{*}\end{array}\right).$

To illustrate the different splittings, we take $(\frac{1}{10},0)$ as the
initial condition, integrate until the final time $t_{f}=500$ and measure the
relative error in energy, $|H(q_{n},p_{n})-H(q_{0},p_{0})|/|H(q_{0},p_{0})|$,
along the trajectory. The step size is taken so that all the methods tested
require the same number of evaluations of the potential (and thus essentially
involve the same computational cost): 1200 evaluations (bottom left panel) and
2400 evaluations (bottom right panel). The schemes we test are the following.
On the one hand, the Störmer–Verlet algorithm (31) (denoted as $S_{2}$ in the
graphs) and the 4th-order Runge–Kutta–Nyström splitting method proposed in
Blanes and Moan (2002) (RKN${}_{6}4$), as representatives of the $T+V$
splitting. On the other hand, the specially adapted schemes (34), denoted
$(2,2)$ and the $(10,6,4)$ integrator presented in Blanes et al. (2013b), both
for the case where $H=H_{1}+H_{2}$, with $H_{2}$ small compared with $H_{1}$.
The notation $(10,6,4)$ refers to the fact that the local error of the method
is of order $\mathcal{O}(\varepsilon
h^{11}+\varepsilon^{2}h^{7}+\varepsilon^{3}h^{5})$ if $H_{2}$ is $\varepsilon$
times smaller than $H_{1}$. For analogy, we label Strang method applied to the
perturbed harmonic oscillator (32) as $(2,2)$. Notice that the splitting (32)
is more advantageous for this initial condition (the simple method (34)
behaves better than the 4th-order scheme), and that the improvement with
respect to the $T+V$ splitting is approximately of the size of $\varepsilon$.
We therefore see that for this type of problem it is possible to construct
integrators providing much more accurate results with the same computational
effort.

Figure 1.1: Simple pendulum. Top: phase space and three trajectories with
initial conditions $(q_{0},p_{0})=(-5,\frac{5}{2}),(1,1),(\frac{1}{10},0)$.
Bottom: relative error in energy committed by different splitting methods
along the solution with initial condition $(q_{0},p_{0})=(\frac{1}{10},0)$ in
the interval $t\in[0,500]$ with (left) 1200 evaluations and (right) 2400
evaluations of the potential.

Since all the schemes are symplectic integrators and the evolution is taking
place in a compact domain, the error in energy is bounded (in contrast with
standard non-symplectic methods, whose error in energy usually grows linearly
with $t$), whereas the error in phase space $(q,p)$ grows linearly when
applied to near-integrable Hamiltonian systems Hairer et al. (2006). One
should notice, however, that if the scheme is conjugate to another more
accurate integrator, then the global error in phase space will remain bounded
for some time interval before it starts growing linearly (see Subsection 4.5
for more details). This feature is illustrated in Figure 1.2: we integrate the
system starting with the same initial condition and final time $t_{f}=500$
with $h=\frac{5}{12}$ and compute the relative error

$\|(q(t_{n}),p(t_{n}))-(q_{n},p_{n})\|/\|(q(t_{n}),p(t_{n}))\|,$

with different schemes. The reference solution $(q(t_{n}),p(t_{n}))$ is
computed numerically with very high accuracy. Specifically, we test the
following integrators: the Lie–Trotter method, eq. (23) (denoted LT in the
graph), and the Störmer–Verlet ($S_{2}$) method for the splitting $H=T+V$, and
scheme $(2,2)$ and (33) (LTpert), which corresponds to the Lie–Trotter method
applied to the perturbed harmonic oscillator (32). We observe that, since LT
and $S_{2}$ are conjugate to each other (see eq. (18)), their errors are quite
similar after some time interval. On the other hand, $(2,2)$ and LTpert (which
are also conjugate to each other) show a different behavior: no linear growth
is visible, and the error of LTpert is larger than that of (2,2) by
approximately the same factor for the whole time interval considered in Figure
1.2. This is related to the fact that $(2,2)$ and LTpert are conjugate to
another method with a local error essentially
$\mathcal{O}(h^{3}\varepsilon^{2})$ (see Subsections 4.5 and 5.7). More
comments on these observations along with additional explanations will be
given in Subsection 4.5.

Figure 1.2: Pendulum. Relative error in phase space for different splitting
methods along the solution with initial condition
$(q_{0},p_{0})=(\frac{1}{10},0)$ in the interval $t\in[0,500]$ and time step
$h=\frac{5}{12}$.

###  1.5 The gravitational N-body problem

Another popular example to illustrate the behavior and performance of
splitting methods corresponds to the important problem in Classical Mechanics
of a planetary system modeled as $N$ bodies (a massive star and $N-1$ planets)
under mutual gravitational Newtonian interaction. This is also a Hamiltonian
system with

$H(q,p)=\sum_{i=0}^{N-1}\frac{1}{2m_{i}}p_{i}^{T}p_{i}-G\sum_{i=1}^{N-1}\sum_{j=0}^{i-1}\frac{m_{i}m_{j}}{\|q_{i}-q_{j}\|}.$
(35)

Here $(q,p)$ denote the ‘supervectors’ composed by the positions
$q_{i}\in\mathbb{R}^{3}$ and momenta $p_{i}\in\mathbb{R}^{3}$ of the ‘Sun’
($i=0$) and the $N-1$ planets ($i=1,\ldots,N-1$) in some Cartesian coordinate
system: $q=(q_{0},q_{1},\ldots,q_{N-1})^{T}$,
$p=(p_{0},p_{1},\ldots,p_{N-1})^{T}$. In (35), $m_{i}$ is the mass of the
$i$-th body, and $G$ is the universal gravitational constant.

Now the equations of motion (20) read

$\begin{split}\frac{dq_{i}}{dt}&=\nabla_{p_{i}}H=\frac{1}{m_{i}}p_{i},\qquad
i=0,\ldots,N-1,\\\ \frac{dp_{i}}{dt}&=-\nabla_{q_{i}}H=-G\sum_{j\neq
i}\frac{m_{i}m_{j}}{\|q_{i}-q_{j}\|^{3}}(q_{i}-q_{j}),\qquad
i,j=0,\ldots,N-1.\end{split}$ (36)

Since the kinetic energy $T(p)$ and the potential energy $V(q)$ are in this
case

$T(p)=\sum_{i=0}^{N-1}\frac{1}{2m_{i}}p_{i}^{T}p_{i},\qquad\quad
V(q)=-G\sum_{i=1}^{N-1}\sum_{j=0}^{i-1}\frac{m_{i}m_{j}}{\|q_{i}-q_{j}\|},$
(37)

it also makes sense to separate the Hamiltonian (35) as $H(q,p)=T(p)+V(q)$, so
that the symplectic Euler and the Stömer–Verlet schemes can be applied to
arbitrary configurations of the bodies. However, this choice is suboptimal for
planetary systems, where planets describe near-Keplerian orbits around the
central star. An alternative procedure taking advantage of the hierarchical
nature of the motion of the planets around the central massive body was first
proposed in Wisdom and Holman (1991), and is known in the literature as the
Wisdom–Holman integration map. It essentially consists in changing coordinates
so that the transformed $H$ can be written as an integrable part $H_{1}$
(corresponding to the Keplerian motion of the planets), and a small
perturbation $H_{2}$ (that accounts for the gravitational interaction of the
planets among themselves), and then applying the second-order scheme (24) to
this new Hamiltonian.

Specifically, Wisdom and Holman (1991) consider a linear canonical change of
variables to rewrite (35) in the so-called Jacobi coordinates
$(\hat{q}_{i},\hat{p}_{i})$, $i=0,1,\ldots,N-1$. Here $\hat{q}_{0}$ is the
position of the center of mass of the system, $\hat{q}_{1}$ is the relative
position of the first planet with respect to the central star, and for
$i=2,\ldots,N-1$, $\hat{q}_{i}$ is the position of the $i$-th planet relative
to the center of mass of the central star and the planets with lower indices.
That is,

$\hat{q}_{0}=\frac{1}{M_{N-1}}\sum_{j=0}^{N-1}m_{j}\,q_{j},\qquad\hat{q}_{i}=q_{i}-\frac{1}{M_{i}}\sum_{j=0}^{i-1}m_{j}\,q_{j},\quad
i=1,\ldots,N-1,$ (38)

where $M_{i}=\sum_{j=0}^{i}m_{j}$ for $i=0,\ldots,N-1$. This can be written in
a more compact way as $\hat{q}=Aq$, where $A$ is an $N\times N$ invertible
matrix with mass-dependent entries,
$\hat{q}=(\hat{q}_{0},\ldots,\hat{q}_{N-1})^{T}$ and
$q=(q_{0},\ldots,q_{N-1})^{T}$.

The conjugate momenta $\hat{p}_{i}$, $i=0,1,\ldots,N-1$, are uniquely
determined by the requirement that the change of variables be canonical (i.e.,
$p=A^{T}\hat{p}$), so that the transformed Hamiltonian function is obtained by
rewriting (35) in the new variables: $\hat{H}(\hat{q},\hat{p})\equiv
H(A^{-1}\hat{q},A^{T}\hat{p})$.

It is straightforward to check that the kinetic energy, expressed in terms of
$\hat{p}$, has the same diagonal structure as in $p$:

$T=\sum_{i=0}^{N-1}\frac{1}{2m_{i}}p_{i}^{T}p_{i}=\sum_{i=0}^{N-1}\frac{1}{2\hat{m}_{i}}\hat{p}_{i}^{T}\hat{p}_{i},$

where

$\hat{m}_{0}=M_{N-1},\qquad\hat{m}_{i}=\frac{M_{i-1}}{M_{i}}m_{i},\qquad
i=1,\ldots,N-1.$

The Hamiltonian in the new variables $\hat{H}(\hat{q},\hat{p})$ can be split
as $\hat{H}(\hat{q},\hat{p})=H_{1}(\hat{q},\hat{p})+H_{2}(\hat{q})$, where

$\displaystyle H_{1}(\hat{q},\hat{p})$
$\displaystyle=\frac{1}{2\hat{m}_{0}}\hat{p}_{0}^{T}\hat{p}_{0}+\sum_{i=1}^{N-1}\big{(}\frac{1}{2\hat{m}_{i}}\hat{p}_{i}^{T}\hat{p}_{i}-\frac{m_{0}m_{i}}{\|\hat{q}_{i}\|}\big{)},$
(39) $\displaystyle H_{2}(\hat{q})$
$\displaystyle=V(q)+G\sum_{i=1}^{N-1}\frac{m_{0}m_{i}}{\|\hat{q}_{i}\|}$
$\displaystyle=Gm_{0}\sum_{i=1}^{N-1}m_{i}\big{(}\frac{1}{\|\hat{q}_{i}\|}-\frac{1}{\|q_{i}-q_{0}\|}\big{)}-G\sum_{i=2}^{N-1}\sum_{j=1}^{i-1}\frac{m_{i}m_{j}}{\|q_{i}-q_{j}\|}$

and $q$ has to be expressed in terms of $\hat{q}$ according to
$q=A^{-1}\hat{q}$.

Observe that the potential energy does not depend on $\hat{q}_{0}$, so
$\hat{p}_{0}$ is constant (in fact, it is the linear momentum of the system)
and therefore we can remove it from $\hat{H}$ if we assume that the center of
the mass is at rest.

Clearly, for fixed $\hat{q}$ and varying mass ratios,

$\hat{q}_{i}=q_{i}-q_{0}+\mathcal{O}(\varepsilon),\quad
H_{2}(\hat{q})=\mathcal{O}(\varepsilon)\quad\mbox{as}\quad\varepsilon\equiv\frac{1}{m_{0}}\max_{1\leq
i\leq N-1}m_{i}\to 0.$

Hamiltonian $H_{1}$ can then be considered as a collection of $N-1$ two-body
problems, and $H_{2}$ as a perturbation. It turns out that the flow
$\varphi_{h}^{[H_{1}]}$ can be computed with the algorithm proposed e.g. in
(Danby 1988, p. 165), whilst $H_{2}$ only depends on $\hat{q}$, and thus its
flow $\varphi_{h}^{[H_{2}]}$ can be explicitly evaluated in an efficient way.
Notice that the number of terms in $H_{1}$ grows linearly with the number of
bodies $N$, whereas the number of terms in $H_{2}$ grows quadratically.

As an illustrative example, we next consider the outer Solar System modeled as
a six-body system with the inner Solar System ($i=0$), the four giant planets
($i=1,2,3,4$), and Pluto ($i=5$), all considered as point masses. The initial
conditions for each planet are taken at Julian time (TDB) 2440400.5 (28 June
1969), obtained from the DE430 ephemerides Folkner et al. (2014) and
normalized so that the center of mass of the system is at rest. A schematic
diagram of the trajectories is shown in Figure 1.3 (top) with the initial
(circles) and final (stars) positions of each object after 200000 days. For
this problem we test the same methods as for the pendulum in Figure 1.1. On
the one hand, Störmer–Verlet (31), $S_{2}$, and the 4th-order
Runge–Kutta–Nyström splitting method of Blanes and Moan (2002), RKN${}_{6}4$
when the Hamiltonian (35) is separated into kinetic and potential energy (37).
On the other hand, the specially adapted schemes (34) (called in this setting
the Wisdom-Holman integrator, and denoted (2,2) as before) and the $(10,6,4)$
integrator presented in Blanes et al. (2013b) when $H$ is expressed as
$H(\hat{q},\hat{p})=H_{1}(\hat{q},\hat{p})+H_{2}(\hat{q})$, with (39) in
Jacobi coordinates. We integrate for a relatively short time interval,
$t_{f}=200000$ days (or approximately 46 periods of Jupyter and two periods of
Pluto) and compute the relative error in the energy with each integrator with
a step size so that all of them require the same number of force evaluations.
The results are displayed in Figure 1.3 with 1200 (bottom left) and 2400
evaluations (bottom right). As in the previous example of the pendulum, the
error in energy remains bounded in all cases, and scheme $(2,2)$ provides an
error energy almost $1000$ times smaller than $S_{2}$. This illustrates the
fact that, by taking a splitting adapted to the structure of the problem and
designing integrators taking these specific features into account, it is
possible to greatly improve the efficiency.

Figure 1.3: Top: Trajectories of the six-body system modeling the outer Solar
System. Bottom: relative error in energy as a function of time for an interval
of 200000 days obtained with different splitting methods with (left) 1200
evaluations and (right) 2400 evaluations of the force.

We display in Figure 1.4 the relative errors in positions
$q=(q_{0},\ldots,q_{5})\in\mathbb{R}^{18}$ along the time integration for the
same methods as in Figure 1.2. As in the pendulum problem, the error for LT is
larger than for $S_{2}$ at the beginning of the integration interval, but they
become similar after some time. Here $(2,2)$ is also more accurate than
$S_{2}$ by a factor significantly smaller than $\varepsilon$. Compared to the
pendulum problem, the error grows linearly right from the beginning for
schemes LTpert and $(2,2)$ (although with a smaller slope than LT and
$S_{2}$). The curve labelled by pLTpert corresponds to the relative error
obtained by applying LTpert with initial conditions
$(\bar{q}_{0},\bar{p}_{0})=\varphi_{h/2}^{[H_{1}]}(q_{0},p_{0})$. Notice that
the phase errors of pLTpert are very similar to those of $(2,2)$. All these
observations will be accounted for in Subsection 4.5.

Figure 1.4: Outer Solar System. Relative error in position as a function of
$t$ for a time interval $[0,t_{f}=200000]$ days with $h=t_{f}/1200$ obtained
with different splitting methods.

###  1.6 The time-dependent Schrödinger equation

The basic object to study the time evolution of a system in Quantum Mechanics
is the (time-dependent) Schrödinger equation. In the case of one particle of
unit mass in a potential $V(x)$, it reads Messiah (1999)

$i\hbar\frac{\partial}{\partial
t}\psi(x,t)=-\frac{1}{2}\Delta\psi(x,t)+V(x)\psi(x,t),\qquad\psi(x,0)=\psi_{0}(x).$
(40)

Here
$\psi:\Omega\subset\mathbb{R}^{3}\times\mathbb{R}\longrightarrow\mathbb{C}$ is
the wave function, representing the state of the system and $\hbar$ is the
reduced Planck constant. The quantity $|\psi(x,t)|^{2}$ represents a
probability density for the position of the particle, in the sense that the
probability of the particle to be located in $S\subset\Omega$ at time $t$ is
$\int_{S}|\psi(x,t)|^{2}dx$. The equation is then defined in the Hilbert space
$L^{2}(\Omega,\mathbb{C})$.

If we introduce the self-adjoint operators $\hat{T}$, $\hat{V}$ acting on
$\psi\in L^{2}(\Omega,\mathbb{C})$ as

$\hat{T}\psi=-\frac{1}{2}\Delta\psi,\qquad\quad\hat{V}\psi=V(x)\psi,$

then, a straightforward calculation shows that
$[\hat{V},[\hat{T},\hat{V}]]\psi=|\nabla V|^{2}\psi$, and therefore

$[\hat{V},[\hat{V},[\hat{T},\hat{V}]]]\psi=0.$ (41)

A standard approach for applying splitting methods in this setting consists in
discretizing first the equation in space. If we consider the one-dimensional
case for simplicity and if the wave function is negligible outside a space
interval $[a,b]$ on the time interval of interest, then we can limit ourselves
to the study of the equation on that finite interval with periodic boundary
conditions Lubich (2008). After rescaling, the periodic interval can always be
restricted to $[-\pi,\pi]$. In this way, the original problem is transformed
into ($\hbar=1$)

$i\frac{\partial}{\partial
t}\psi(x,t)=-\frac{1}{2}\frac{\partial^{2}\psi}{\partial
x^{2}}(x,t)+V(x)\psi(x,t),\qquad x\in[-\pi,\pi],$ (42)

with $\psi(-\pi,t)=\psi(\pi,t)$ for all $t$.

The wave function is then approximated by a trigonometric polynomial $u(x,t)$
whose coefficients are obtained by requiring that the approximation satisfies
(42) in a grid of $M$ equispaced points $x_{j}=-\pi+j\cdot 2\pi/(M-1)$ on the
interval $[-\pi,\pi]$. The vector
$u=(u_{0},\ldots,u_{M-1})^{T}\in\mathbb{C}^{M}$ formed by the grid values
$u_{j}\approx\psi(x_{j},t)$, $j=0,1,\ldots,M-1$, then verifies the
$M$-dimensional linear ODE

$i\frac{d}{dt}u(t)=H\,u(t)=(T+V)u(t),\qquad u(0)=u_{0}\in\mathbb{C}^{M}.$ (43)

Here $V=\mbox{diag}(V(x_{j}))$ and $T=-\frac{1}{2}D$, where $D$ is the second-
order periodic spectral differentiation matrix Trefethen (2000). As is well
known, $Tu=\mathcal{F}^{-1}D_{T}\mathcal{F}u$, where $\mathcal{F}$ and
$\mathcal{F}^{-1}$ are the forward and backward discrete Fourier transform,
and $D_{T}$ is also diagonal. The transformations $\mathcal{F}$ and
$\mathcal{F}^{-1}$ are computed with the fast Fourier transform (FFT)
algorithm, requiring $\mathcal{O}(M\log M)$ operations.

Notice that solving equations $iu^{\prime}=Tu$ and $iu^{\prime}=Vu$ is done
trivially by using exponentials of diagonal matrices and FFTs, namely,

$\left(\mathrm{e}^{\tau V}u\right)_{j}=\mathrm{e}^{\tau
V(x_{j})}u_{j},\qquad\quad\mathrm{e}^{\tau
T}u=\mathcal{F}^{-1}\mathrm{e}^{\tau D_{T}}\mathcal{F}\,u,$

for a time step $h$, with $\tau=-ih$. Therefore, splitting methods constitute
a valid alternative to approximate the solution $u(t)=\mathrm{e}^{\tau
H}u_{0}=\mathrm{e}^{\tau(T+V)}u_{0}$, which may be prohibitively expensive to
evaluate for large values of $M$. Thus, the Lie–Trotter scheme reads

$\mathrm{e}^{\tau(T+V)}=\mathrm{e}^{\tau T}\,\mathrm{e}^{\tau
V}+\mathcal{O}(\tau^{2}),$ (44)

whereas the 2nd-order Strang splitting constructs the numerical approximation
$u_{n+1}$ at time $t_{n+1}=t_{n}+\Delta t$ by

$u_{n+1}=\mathrm{e}^{\tau/2V}\,\mathrm{e}^{\tau
T}\,\mathrm{e}^{\tau/2V}\,u_{n}.$ (45)

The resulting scheme is called the _split-step Fourier_ method in the chemical
literature Feit et al. (1982), and has some remarkable properties. In
particular, it is unitary and time-symmetric Lubich (2008), as is the exact
solution $\mathrm{e}^{\tau H}$.

Relation (41) still holds for the matrices $T$ and $V$ if the number of points
$M$ in the space discretization is sufficiently large, and in fact $[V,[T,V]]$
is diagonal if the derivatives of the potential are computed first and then
evaluated at the corresponding space grid.

We next illustrate the procedure with the one-dimensional double-well
potential

$V(x)=\frac{1}{80}(x^{2}-20)^{2},$ (46)

and the initial wave function
$\psi(x,0)=\psi_{0}(x)=\sigma\cos^{2}(x)\mathrm{e}^{-\frac{1}{2}(x-1)^{2}}$,
where $\sigma$ is an appropriate normalizing constant. We take $M=256$
discretization points on the interval $x\in[-13,13]$ and integrate the
resulting linear ODE (43) in the interval $0\leq t\leq t_{f}=Nh=10$. Figure
1.5 shows in the left panel $|\psi_{0}(x)|^{2}$, $|\psi(x,t_{f})|^{2}$ and the
potential $V(x)$, whereas the right panel shows an efficiency diagram.
Specifically, we display the error in energy measured at the final time,
$|u_{N}^{T}Hu_{N}-u_{0}^{T}Hu_{0}|$, as a function of the number of FFT calls
(and its inverse) as an estimate of the computational effort of each method.
The lines correspond to the Strang splitting (45), $S_{2}$, with time steps
$h=10/2^{k}$, $k=1,2,\ldots,12$, the 4th-order RKN splitting method
RKN${}_{6}4$ from Blanes and Moan (2002), already illustrated in Figures 1.1
and 1.3, and another 4th-order scheme including the double commutator
$[V,[T,V]]$ into its formulation (denoted RKNm${}_{4}4$). In this diagram the
slope of each line for sufficiently small $h$ indicates the order of the
scheme. As in the previous examples, by taking into account the specific
features of the problem at hand, it is possible to construct more accurate and
efficient numerical approximations.

Figure 1.5: Time-dependent Schrödinger equation with a double-well potential.
Left: $V(x)$, initial and final wave function with $M=256$ discretization
points. Right: relative error in energy at the final time vs. number of FFTs
for different values of the time step obtained with the Strang method $S_{2}$,
and two 4th-order splitting schemes: one involving 6 evaluations of $V$
(RKN${}_{6}4$) and another with 4 evaluations of $V$, and incorporating in
addition the double commutator $[V,[T,V]]$ (RKNm${}_{4}4$).

The error in energy also remains bounded for these unitary integrators, as in
the previous examples involving classical Hamiltonian problems, whereas the
error in the wave function grows linearly with $t$, unless the scheme is
conjugate to another more accurate one, in which case it is bounded for some
time before linear growth takes place. To illustrate this feature, in Figure
1.6 we depict how the error in the solution $\|u(t_{n})-u_{n}\|$ evolves with
time for a longer integration interval $t\in[0,1000]$, with step size
$h=\frac{1}{20}$. As usual, the reference solution is computed numerically
with sufficiently high accuracy, and the tested schemes are as follows:
Lie–Trotter, eq. (44) (LT), Strang, eq. (45), and a variant of Strang
involving a double commutator, namely

$S_{2}m:\qquad
u_{n+1}=\mathrm{e}^{\tau/2V+(\tau^{3}/48)[V,[T,V]]}\,\mathrm{e}^{\tau
T}\,\mathrm{e}^{\tau/2V+(\tau^{3}/48)[V,[T,V]]}\,u_{n}.$ (47)

This scheme does not require any additional FFTs and in addition is actually
conjugate to a method of order 4. We observe that LT approaches $S_{2}$ after
a transition time, whereas the error for $S_{2}m$ remains bounded for the
whole interval. In fact, its linear error growth only appears when the time
interval is doubled.

Figure 1.6: Time-dependent Schrödinger equation with a double-well potential:
relative error in the wave function as a function of $t$ with $h=1/20$ for
Strang ($S_{2}$), Lie-Trotter (LT) and scheme $S_{2}m$ involving a double
commutator and conjugate to a method of order 4.

###  1.7 Splitting methods as geometric numerical integrators

As the examples of subsections 1.2-1.6 illustrate, splitting schemes, even of
low order of consistency such as Lie–Trotter and Strang methods, preserve by
construction structural properties of the exact solution, such as
symplecticity (in classical Hamiltonian dynamics) and unitarity (in quantum
evolution problems). This feature gives them qualitative superiority with
respect to other standard integrators in practice, especially when long time
intervals are concerned. In this sense, splitting methods constitute an
important class of _geometric numerical integrators_.

Although the idea that numerical integrators applied to an ordinary
differential equation should preserve as many properties of the system as
possible has been implicitly assumed since the early days of numerical
analysis, it is fair to say that in the classical consistency/stability
approach the emphasis has been on other issues. In particular, the goal has
mainly been to compute the solution of (1) at time $t_{f}=Nh$ with a global
error $\|x_{N}-x(t_{f})\|$ smaller than a prescribed tolerance and as
efficiently as possible. To do that one chooses the class of method (one-step,
multistep, extrapolation, etc.), the order (fixed or adaptive) and the time
step (constant or variable). This approach has proved to be very fruitful,
giving rise to highly tuned and thoroughly tested software packages generally
available to solve a great variety of problems.

On the other hand, there are special types of problems arising in many fields
of science and applied mathematics that possess an underlying geometric
structure which influences the qualitative character of their solutions, and
so one aims naturally to construct numerical approximations that preserve this
structure. Classical Hamiltonian systems such as those illustrated previously
constitute a case in point. It turns out, however, that many numerical
integrators included in standard software packages do not take into account
these distinctive features of the equations to be solved, and the question is
whether it is possible to design, analyze and apply new schemes providing
approximate solutions that share one or several geometric properties with the
exact solution. This is precisely the realm of Geometric Numerical
Integration, a terminology introduced in Sanz-Serna (1997).

According to McLachlan and Quispel (2006), _‘Geometric integration’ is the
term used to describe numerical methods for computing the solution of
differential equations, while preserving one or more physical/mathematical
properties of the system_ exactly _(i.e., up to round-off error)_. Thus,
rather than primarily taking into account prerequisites such as consistency
and stability, the aim is to reproduce the qualitative features of the
solution of the differential equation being discretized, in particular its
geometric properties, such as the symplectic character (for Hamiltonian
systems) and unitarity (quantum mechanics), but also the phase-space volume
(for divergence-free vector fields), time-reversal symmetries, first integrals
of motion (energy, linear and angular momentum), Casimirs, Lyapunov functions,
etc. In these structure-preserving methods one tries to incorporate as many of
these properties as possible and, as a result, they exhibit improved
qualitative behavior. In addition, they typically allows for a significantly
more accurate integration for long-time intervals than with general-purpose
methods Hairer et al. (2006), Blanes and Casas (2016).

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.), it is fair to say that the interest in
splitting has revived with the advent of Geometric Numerical Integration, and
new and very efficient schemes have been put to use to solve a wide variety of
problems. The reason is clear: if the problem (3) has some property that is
deemed to preserve (symplectic, unitary, volume-preserving, etc.) and each
subproblem $x^{\prime}=f_{j}(x)$, $x(0)=x_{0}$ can be integrated exactly or by
means of a numerical method preserving these properties, then the splitting
method constructed by composing the solution of the subproblems is also
symplectic, unitary, volume-preserving, etc. In other words, splitting methods
provide by construction approximations lying in the same group of
diffeomorphisms as the system $x^{\prime}=f(x)$ McLachlan and Quispel (2002).
Here we are assuming of course that each subproblem $x^{\prime}=f_{j}(x)$
possesses the same characteristic feature as the total problem considered.

###  1.8 Relevance of splitting methods

Given a certain differential equation $x^{\prime}=f(x)$, the application of
splitting methods to solve the corresponding initial value problem involves
three main steps McLachlan and Quispel (2002):

  * (1)

Choosing the set of terms $f_{j}$ such that $f=f_{1}+\cdots+f_{m}$. Different
ways of decomposing $f$ may give rise to integrators with different
qualitative behavior and efficiency, as we have seen in the previous examples.

  * (2)

Solving each subproblem $x^{\prime}=f_{j}(x)$ either exactly or approximately.

  * (3)

Combining these solutions to get an approximation for the original overall
problem.

Being such a simple idea, it is hardly surprising that the splitting principle
can be used in so many different settings. In particular, one may

  * •

split the differential equation into linear and nonlinear parts;

  * •

in an ODE describing a Hamiltonian system with an additional small
dissipation, separate the Hamiltonian part and the dissipation;

  * •

decompose into parts describing different physical processes, as for example
diffusion and reaction in partial differential equations;

  * •

get approximations to the original problem by solving one space direction at a
time (dimensional splitting in PDEs).

In addition, splitting methods possess some advantages concerning their
implementation, in particular the following.

  * •

They are typically explicit.

  * •

Their storage requirements are quite small. The algorithms are sequential and
the solution at intermediate stages can be stored in the solution vectors.

  * •

Programming higher order schemes is not more difficult than Lie–Trotter and
Strang splitting methods, at least in the context of ODEs. Usually, a few more
lines of code is all what is required to deal with the additional stages.

  * •

As stated earlier, they can preserve a wide variety of structures possessed by
the differential equation.

They also present some disadvantages, of course. Among them, we can mention
the following.

  * •

Splitting schemes of order three or higher necessarily involve negative
coefficients. In other words, they require sub-steps that go backwards in
time, and this has severe repercussions when applying them to, for instance,
reaction-diffusion equations (see Section 6).

  * •

Although it is possible to construct accurate high order splitting methods,
stability can be an issue, in the sense that their stability interval might be
reduced to render them useless in practice. This aspect has to be taken
seriously into account when designing new methods.

  * •

Ordinary splitting does not capture the correct steady state solutions (where
these solutions exist) MacNamara and Strang (2016), in the sense that the
numerical solutions obtained converge to limits that are not steady-state
solutions, but just approximations of them. This can lead to unacceptable
errors, for instance in the simulation of combustion. In this context,
balanced splitting techniques have been introduced to correct this flaw Speth
et al. (2013).

Splitting methods constitute an important tool in different areas of science
where the evolution of systems is governed by differential equations. In
addition to Hamiltonian systems, they can be successfully applied in the
numerical study of Poisson systems, systems possessing integrals of motion
(such as energy and angular momentum) and systems with (continuous, discrete
and time-reversal) symmetries. In fact, splitting methods have been designed
(often independently) and extensively used in fields as distant as molecular
dynamics, simulation of storage rings in particle accelerators, celestial
mechanics, astronomy, quantum (statistical) mechanics, plasma physics,
hydrodynamics, and Markov Chain Monte Carlo methods.

Operator-splitting methods also appear outside the realm of differential
equations, and in particular in optimization, in a variety of different
special forms and different denominations (gradient-projection, proximal-
gradient, alternating direction method of multipliers or ADMM, split Bregman,
etc.). All of them are related with special types of splitting methods, such
as Douglas–Rachford and Peaceman–Rachford schemes. More details can be found
in several contributions collected in the comprehensive book Glowinski et al.
(2016a).

###  1.9 Some historical remarks

There is ample consensus that the beginning of splitting is related to the
product formula (5). What is not so clear is the origin of the formula itself.
Thus, (Reed and Simon 1980, p. 295) establish it as ‘the classical theorem of
Lie’, but give no exact source, whereas Chorin et al. (1978) call it ‘the 1875
formula of S. Lie’, citing the classical treatise Lie (1888), and, based on
this reference, Glowinski et al. (2016a) even ascribe to Lie himself ‘the
first operator-splitting scheme recorded in history’. The problem is that the
reference Lie (1888) is clearly not from 1875, and it is not evident (at least
to us) that this formula appears there explicitly.

On the other hand, as pointed out in Cohen et al. (1982), the result (5) can
be found in several references published during the 1950s, namely Butler and
Friedman (1955), Golden (1957), whereas it was Trotter who generalized it to
self-adjoint linear operators Trotter (1959) without mentioning Lie or these
previous references. Subsequently, formula (5), even in the matrix case, has
been attributed to Trotter (Bellman 1970, p. 181). In view of the situation,
we believe we are not committing an act of historical injustice by referring
to the approximation (6) and Algorithm 1.1 as the _Lie–Trotter_ scheme.

With respect to the splitting method (9), it first appeared in print in Strang
(1968) as an alternative way to solve multidimensional problems with one-
dimensional operators. We have already seen that, when applied to Hamiltonian
systems of the form $H(q,p)=T(p)+V(q)$, it leads to the  Algorithm 1.4 when
composing the flows associated to $T(p)$ and $V(q)$. It is called the
_Störmer–Verlet method_ since it was used by the astronomer C. Störmer in 1907
in his computations of the motion of ionized particles in the Earth’s magnetic
field Störmer (1907), and by L. Verlet in molecular dynamics Verlet (1967). It
is also referred to as the _leapfrog method_ in the context of PDEs describing
wave propagation and as the _Wisdom–Holman method_ when applied to the
splitting (39) Wisdom and Holman (1991) In fact, it can be found in several
classical references, the oldest one being perhaps Newton’s _Principia_ 222
_Philosophiae Naturalis Principia Mathematica_ , Book 1, Section 2,
Proposition 1.. For a detailed account the reader is referred to the
enlightening review Hairer et al. (2003).

As well as symplecticity when applied to Hamiltonian systems, the
Störmer–Verlet method preserves many other geometric properties of the exact
flow associated with an ordinary differential equation. This includes the
preservation of all linear first integrals (such as the linear momentum), and
quadratic first integrals of the form $I(q,p)=p^{T}Cq$ for Hamiltonian
systems, where $C$ is a symmetric matrix. In other words, $I(q,p)$, computed
along the numerical trajectory, is constant. A classical example is the
angular momentum in $N$-body problems if the forces only depend on the
distances of the particles.

All these favorable properties, in addition to its optimal stability property
and reversibility, help to understand why this method is probably the most
used splitting scheme and geometric integrator, especially in molecular
dynamics Schlick (2010), condensed matter simulations Ceperley (1995),
sampling with the hybrid Monte Carlo method Neal (2011), etc.

The convenience of designing numerical integration methods that, by
construction, preserve the symplectic structure when applied to Hamiltonian
systems was duly recognized during the 1950s in the field of accelerator
physics. Thus, in the words of an early pioneer, _if one wishes to examine
solutions to differential equations, adoption of a ‘Hamiltonian’ or
‘canonical’ integration algorithm would be reassuring_ Laslett (1986). This
was the point of view adopted by De Vogelaere in his pioneering paper de
Vogelaere (1956)333See also Skeel and Cieśliński (2020) for the context of the
work and the preprint itself, typeset in LaTeX., where he devoted himself to
the task of designing _a method of integration which, if there was no round-
off error, would give a solution with the contact transformation property_.
Here “contact transformation” has to be understood as “symplectic
transformation”. The first order schemes proposed in de Vogelaere (1956),
although implicit in general, turn out to be explicit when $H(q,p)=T(p)+V(q)$,
in which case they reproduce the symplectic Euler schemes (29) and (30).

It was another accelerator physicist, R. Ruth, who presented in 1983 what is
probably the first splitting method of order three Ruth (1983). This paper can
be considered as the actual starting point in the systematic exploration of
symplectic integrators along several parallel avenues: (i) the use of
generating functions in the context of Hamiltonian mechanics to produce
appropriate canonical transformations approximating the exact flow in each
integration step Feng and Qin (1987), Channell and Scovel (1990); (ii) the
conditions that Runge–Kutta methods have to satisfy to be symplectic Sanz-
Serna (1988), Lasagni (1988), Suris (1988); (iii) the design of explicit
symplectic methods of order 4 and higher for Hamiltonian systems that can be
split into two pieces which can be solved exactly when considered as
independent systems Neri (1988), Forest (1992), with the help of the Lie
formalism. This approach was further elaborated in Yoshida (1990). Working in
the context of the hybrid Monte Carlo algorithm for dynamical fermions, a
splitting method of order 4 was also independently proposed around the same
time in Campostrini and Rossi (1990). In parallel developments, what is now
called the Suzuki–Yoshida composition technique for increasing the order of
numerical integrators appeared in Creutz and Gocksch (1989) for Monte Carlo
simulations and in Suzuki (1990) and Yoshida (1990).

We should also mention the papers de Raedt and de Raedt (1983) and Takahashi
and Imada (1984), who pioneered the use of double commutators to get
approximations of higher order than those obtained by the Störmer–Verlet
method in path integral Monte Carlo simulations: in the first case by
constructing a 4th-order splitting scheme, and in the second, a method that it
is also of order four by conjugation. In fact, scheme (47) is conjugate to the
one proposed in Takahashi and Imada (1984). The paper Ruth (1983) also
presents a third-order method using double commutators.

That splitting and composition methods could be used to construct integrators
for problems evolving in groups other than the symplectic group was emphasized
in Forest and Ruth (1990) and further developed in Feng (1992), with the aim
of constructing schemes able to preserve different structures.

The 1990s saw a dramatic increase in the interest and applications of
splitting integrators in several fields, often with spectacular results. We
should mention in particular those achieved in Wisdom and Holman (1991),
revealing the existence of chaotic phenomena in the Solar System by
numerically integrating the planetary equations of motion over very large time
intervals.

The state of the art of splitting methods in the context of geometric
numerical integration was masterfully summarized in 2002 by McLachlan and
Quispel in their review paper McLachlan and Quispel (2002), which has greatly
influenced subsequent investigation in the field, as testified by its growing
number of citations over the years, in many different areas.

Among the huge number of published works on splitting methods, the following
surveys are worth highlighting.

  * •

The monograph Yanenko (1971) (an English translation of the Russian edition
published in 1967) is perhaps the first to be devoted to the method of
splitting (or _method of fractional steps_) for solving “complicated problems
of mathematical physics in several variables”. Those include the numerical
treatment of parabolic and hyperbolic equations, as well as boundary value
problems for the Laplace and Poisson equations, and several applications in
elasticity theory and hydrodynamics. It is based on the early contributions of
Peaceman, Rachford, Douglas, and several authors of the Soviet school
(Dyakonov, Marchuk, Samarskii, Yanenko, and others).

  * •

The exhaustive review article Marchuk (1990) included in the Volume I of the
Handbook of Numerical Analysis Ciarlet and Lions (1990). It can be seen as an
update of the previous work, with a systematic study of operator splitting and
alternating direction methods for solving linear and nonlinear partial
differential equations. It includes convergence analyses and new applications
to problems in hydrodynamics, meteorology and oceanography.

  * •

The review paper McLachlan and Quispel (2002), mainly focused on the
application of splitting methods as geometric numerical integrators for
various classes of ordinary differential equations. In that context, a
classification of ODEs and their integration methods into different categories
is carried out, and the question of how to decompose a given vector field into
much simpler vector fields is also examined, as well as the composition of
these elementary flows.

  * •

Books and monographs dealing with geometric numerical integration and
structure-preserving algorithms contain plenty of material on splitting
methods. Among them, we can cite the influential work Sanz-Serna and Calvo
(1994), the canonical reference Hairer et al. (2006), as well as Leimkuhler
and Reich (2004), Leimkuhler and Matthews (2015) and Blanes and Casas (2016).

  * •

The multi-author book Glowinski et al. (2016a) constitutes an excellent
illustration of the ample scope and wide range of applications that toay’s
operator-splitting methods are able to deal with. These include the numerical
solution of problems modeled by linear and nonlinear partial differential
equations and inequalities, problems in information sciences and image
processing, and large-scale optimization problems, among others.

###  1.10 Plan of the paper

In this paper we will focus on splitting methods applied to evolutionary
problems, mostly described by ordinary differential equations. These can
directly model the problem one is interested in, or they can result from
evolutionary PDEs previously discretized in space. Particular attention will
be addressed to problems possessing special properties, very often from a
geometric origin, that are worth preserving by the numerical methods. In so
doing, we will follow a strategy similar to that in Blanes et al. (2008b),
trying to avoid any duplication of the material already collected in the
classic references cited above, and including new results, schemes and
applications which have appeared in the literature during the last few years.

In particular, no general rule is provided here on how to split the defining
operator in (1). As mentioned earlier, this issue is further analyzed in
McLachlan and Quispel (2002), and in fact some of the open problems listed
there are related with it. We have already seen in the examples provided in
this section that several splittings of the same problem are possible, often
leading to methods with very different performances. Moreover, in certain
cases, the original system has several geometric properties that are
simultaneously preserved along the evolution, whereas different splittings
preserve different properties and it is generally difficult to find one
splitting that preserves most of them.

With these considerations in mind, the rest of the paper is organized as
follows. In Section 2 we first review the general composition technique and
then provide a detailed analysis of the order conditions required by splitting
and composition methods to achieve a given order of accuracy. There are some
relevant problems, however, whose particular structure allows one to design
specially adapted methods, and some of them are reviewed in Section 3, where
we also show how to adapt existing splitting methods to non-autonomous
problems.

We summarize in Section 4 some of the qualitative properties possessed by
splitting methods in the context of geometric numerical integration of
ordinary differential equations, with special attention to the idea of
processing, whereas Section 5 is devoted to the treatment of highly
oscillatory problems.

Splitting methods are particularly well adapted to deal with partial
differential equations whose defining operator contains contributions coming
from very different physical sources, and so they have a long history in this
area. Section 6 contains a brief survey, with special emphasis on Schrödinger
equations and general parabolic evolution equations. The existence of negative
coefficients in the methods, however, leads to an order barrier for parabolic
equations, and Section 7 reviews splitting methods with complex coefficients
as a possible way to overcome this order barrier.

In Section 8 we present an extended list of existing methods, classifying them
into different families and giving the appropriate references. Their
corresponding coefficients are also provided as supplementary material at the
website

https://www.gicas.uji.es/SplittingMethods.html

These methods are numerically tested on simple examples in the Appendix.
Finally, some relevant applications of splitting methods in different fields
are discussed in Section 9.

##  2 High order splitting and composition methods

The Lie–Trotter and Strang splitting methods, in spite of their low order of
accuracy, provide a fairly good description of the systems they are
approximately solving. In fact, for many problems, including molecular
dynamics applications and reaction-diffusion equations, Verlet and Strang
splitting are the most popular integrators, perhaps an illustration that,
according to MacNamara and Strang (2016), “it is a meta-theorem of numerical
analysis that second order methods often achieve the right balance between
accuracy and complexity”. There are other areas, however, where a higher
degree of precision is required, in addition to the preservation of
qualitative properties. A classical example is the long-term numerical
integration of the Solar System, both forwards (e.g to analyze the existence
of chaos Laskar (1989), Sussman and Wisdom (1992)) and backwards in time (to
study the insolation quantities of the Earth Laskar et al. (2004)). Thus, in
this section, after reviewing a general technique to get high order
integrators by composing low order ones, we analyze from different
perspectives the order conditions a splitting method has to satisfy to achieve
a given order. Such analysis allows one to provide complementary information
about the integrators: number of order conditions, explicit expressions,
remainders in the asymptotic expansions, etc.

###  2.1 Raising the order by composition

####  2.1.1 Composition of Strang maps

Starting from the Strang splitting
$S_{h}^{[2]}=\varphi_{h/2}^{[1]}\circ\varphi_{h}^{[2]}\circ\varphi_{h/2}^{[1]}$,
the composition

$\psi_{h}=S_{\gamma_{s}h}^{[2]}\circ S^{[2]}_{\gamma_{s-1}h}\circ\cdots\circ
S^{[2]}_{\gamma_{1}h}$ (48)

is at least of order three if

$\sum_{j=1}^{s}\gamma_{j}=1,\qquad\mbox{ and
}\qquad\sum_{j=1}^{s}\gamma_{j}^{3}=0.$ (49)

The smallest value of $s$ for which equations (49) admit real solutions is
$s=3$. In that case, by imposing the symmetry $\gamma_{1}=\gamma_{3}$, we
indeed get a method of order 4, sometimes called the _triple jump_ :

$S_{h}^{[4]}=S_{\gamma_{3}h}^{[2]}\circ S_{\gamma_{2}h}^{[2]}\circ
S_{\gamma_{1}h}^{[2]},\qquad\mbox{ with
}\qquad\gamma_{1}=\gamma_{3}=\frac{1}{2-2^{1/3}},\qquad\gamma_{2}=1-2\gamma_{1}.$
(50)

In general, the recursion

$S_{h}^{[2k]}=S_{\gamma_{1}h}^{[2k-2]}\circ S_{(1-2\gamma_{1})h}^{[2k-2]}\circ
S_{\gamma_{1}h}^{[2k-2]},\qquad\mbox{ with
}\qquad\gamma_{1}=\frac{1}{2-2^{1/(2k-1)}}$ (51)

can be used to get methods of arbitrarily high order $2k$ ($k\geq 2$) Creutz
and Gocksch (1989) starting from the Strang map $S^{[2]}_{h}$ (notice that
such methods can be written in the form (48)). The price to be paid is the
existence of large positive and negative coefficients $\gamma_{j}$ and the
great number of elementary flows in (51) for high orders. The alternative
formed by the five maps composition (_quintuple jump_)

$S_{h}^{[2k]}=S_{\gamma_{1}h}^{[2k-2]}\circ S_{\gamma_{1}h}^{[2k-2]}\circ
S_{(1-4\gamma_{1})h}^{[2k-2]}\circ S_{\gamma_{1}h}^{[2k-2]}\circ
S_{\gamma_{1}h}^{[2k-2]},\qquad\gamma_{1}=\frac{1}{4-4^{1/(2k-1)}}$ (52)

also gives methods $S_{h}^{[2k]}$ of order $2k$ of the form (48) with
relatively smaller coefficients $\gamma_{j}$ but even larger numbers of
elementary flows.

In general, other choices for the coefficients $\gamma_{j}$ in (48) are more
appropriate if one is interested in achieving orders $\geq 6$ with a lower
number of elementary flows and relatively small coefficients.

Condition (49), using the approach based on linear differential operators
discussed in Subsection 1.2, can be derived as follows: the Lie transformation
associated with the Strang map $S_{h}^{[2]}$, eq. (16), can be written as
$\mathrm{e}^{h/2F_{1}}\mathrm{e}^{hF_{2}}\mathrm{e}^{h/2F_{1}}=\mathrm{e}^{Y(h)}$,
where

$\displaystyle Y(h)$ $\displaystyle=\sum_{n=1}^{\infty}h^{2n-1}Y_{2n-1}$
$\displaystyle:=\log(\mathrm{e}^{h/2F_{1}}\mathrm{e}^{hF_{2}}\mathrm{e}^{h/2F_{1}})$
$\displaystyle=h(F_{1}+F_{2})-\frac{h^{3}}{24}[F_{1},[F_{1},F_{2}]]-\frac{h^{3}}{12}[F_{2},[F_{1},F_{2}]]+\cdots,$

that is, $Y_{1}=F_{1}+F_{2}$,
$Y_{3}=-\frac{1}{24}[F_{1},[F_{1},F_{2}]]-\frac{1}{12}[F_{2},[F_{1},F_{2}]]$,
and for each $n>2$, $Y_{2n-1}$ is certain linear combination of $(2n-1)$-fold
commutators of $F_{1}$ and $F_{2}$. In consequence, the Lie transformation
$\Psi(h)$ of (48) formally satisfies

$\Psi(h)=\mathrm{e}^{Y(\gamma_{1}h)}\cdots\,\mathrm{e}^{Y(\gamma_{s}h)},$ (53)

so that $g(\psi_{h}(x))=(\Psi(h)g)(x)$ for any $x\in\mathbb{R}^{D}$ and any
smooth function $g:\mathbb{R}^{D}\to\mathbb{R}$. It is straightforward to
check that

$\begin{split}\Psi(h)&=\mathrm{e}^{\gamma_{1}hY_{1}+\gamma_{1}^{3}h^{3}Y_{3}+\cdots}\cdots\,\mathrm{e}^{\gamma_{s}hY_{1}+\gamma_{s}^{3}h^{3}Y_{3}+\cdots}\\\
&=\mathrm{e}^{h\left(\sum_{j=1}^{s}\gamma_{j}\right)(F_{1}+F_{2})}+h^{3}\left(\sum_{j=1}^{s}\gamma_{j}^{3}\right)Y_{3}+\mathcal{O}(h^{4}).\end{split}$
(54)

This shows that the composition (48) of Strang maps is of order at least 3 if
condition (49) holds.

In fact, (54) is also true if in (48) the Strang map $S^{[2]}_{h}$ is replaced
by any second order time-symmetric map. Furthermore, the triple jump recursion
(51) (resp. the quintuple jump recursion (52)) also gives rise to $2k$th order
maps $S_{h}^{[2k]}$ starting from an arbitrary time-symmetric second order map
$S_{h}^{[2]}$. Indeed, this is a consequence of the following four statements:

  * (1)

Given an arbitrary near-identity map
$\chi_{h}:\mathbb{R}^{D}\to\mathbb{R}^{D}$ (i.e.,
$\chi_{h}(x)=x+\mathcal{O}(h)$ as $h\to 0$), there exists a series

$Y(h)=\sum_{n\geq 1}h^{n}Y_{n}$

of (first order) differential operators acting on smooth functions such that
formally, $g(\chi_{h}(x))=(\mathrm{e}^{Y(h)}g)(x)$ for each
$x\in\mathbb{R}^{D}$ and $g\in C^{\infty}(\mathbb{R}^{D},\mathbb{R})$.
Moreover, $\chi_{h}$ is an $r$th order integrator for the ODE system
$x^{\prime}=f_{1}(x)+f_{2}(x)$ if and only if

$Y_{1}=F_{1}+F_{2},\qquad Y_{n}=0\quad\mbox{for}\quad 2\leq n\leq r.$ (55)

This statement can be proved as follows. Given a basic integrator
$\chi_{h}:\mathbb{R}^{D}\to\mathbb{R}^{D}$, consider the linear differential
operators $X_{n}$ ($n\geq 1$) acting on smooth functions $g\in
C^{\infty}(\mathbb{R}^{D},\mathbb{R})$ as

$X_{n}g(y)=\left.\frac{1}{n!}\frac{d^{n}}{dh^{n}}\right|_{h=0}g(\chi_{h}(y)),\qquad
y\in\mathbb{R}^{D},$ (56)

so that formally $g(\chi_{h}(x))=(X(h)g)(x)$, where

$X(h)=I+\sum_{n\geq 1}h^{n}X_{n},$ (57)

and $I$ denotes the identity operator. Each $X_{n}$ is an $n$th-order
differential operator. Thus, the integrator $\chi_{h}$ is of order $r$ if

$X_{n}=\frac{1}{n!}(F_{1}+F_{2})^{n},\qquad 1\leq n\leq r.$

Now consider the series of differential operators

$Y(h)=\sum_{n\geq 1}h^{n}Y_{n}:=\log(X(h))=\sum_{m\geq
1}\frac{(-1)^{m+1}}{m}\left(hX_{1}+h^{2}X_{2}+\cdots\right)^{m},$

that is,

$Y_{n}=\sum_{m\geq
1}^{n}\frac{(-1)^{m+1}}{m}\sum_{j_{1}+\cdots+j_{m}=n}X_{j_{1}}\cdots
X_{j_{m}},$

so that $X(h)=\exp(Y(h))$, and formally, $g(\chi_{h}(x))=(\exp(Y(h))g(x)$. It
can be shown that each $Y_{n}$ is a first order differential operator.
Clearly, the basic integrator is of order $r$ if (55) holds.

  * (2)

The map $\chi_{h}$ is time-symmetric if and only if
$Y(h)=hY_{1}+h^{3}Y_{3}+\cdots$, which implies that time-symmetric methods
_are necessarily of even order_. Indeed, for the adjoint integrator
$\chi^{*}_{h}=\chi_{-h}^{-1}$, one obviously gets
$g(\chi^{*}_{h}(x))=\mathrm{e}^{-Y(-h)}g(x)$. Hence, $\chi_{h}$ is time-
symmetric if and only if $-Y(-h)\equiv Y(h)$.

  * (3)

If $S_{h}^{[2k-2]}$ is a time-symmetric integrator of order at least $2k-2$,
then the composition

$\psi_{h}=S_{\gamma_{s}h}^{[2k-2]}\circ
S^{[2k-2]}_{\gamma_{s-1}h}\circ\cdots\circ S^{[2k-2]}_{\gamma_{1}h}$ (58)

is of order at least $2k-1$ if

$\sum_{j=1}^{s}\gamma_{j}=1,\qquad\quad\sum_{j=1}^{s}\gamma_{j}^{2k-1}=0.$
(59)

In fact, if $\mathrm{e}^{h(F_{1}+F_{2})+h^{2k-1}Y_{2k-1}+\cdots}$ is the Lie
transform of $S_{h}^{[2k-2]}$, then $g(\psi_{h}(x))=(\Psi(h)g)(x)$ for any
$x\in\mathbb{R}^{D}$ and any smooth function $g:\mathbb{R}^{D}\to\mathbb{R}$,
where

$\displaystyle\Psi(h)$
$\displaystyle=\mathrm{e}^{\gamma_{1}h(F_{1}+F_{2})+\gamma_{1}^{2k-1}h^{2k-1}Y_{2k-1}+\mathcal{O}(h^{2k})}\cdots\,\mathrm{e}^{\gamma_{s}h(F_{1}+F_{2})+\gamma_{s}^{2k-1}h^{2k-1}Y_{2k-1}+\mathcal{O}(h^{2k})}$
$\displaystyle=\mathrm{e}^{\left(\sum_{j=1}^{s}\gamma_{j}\right)(F_{1}+F_{2})}+h^{2k-1}\left(\sum_{j=1}^{s}\gamma_{j}^{2k-1}\right)Y_{2k-1}+\mathcal{O}(h^{2k}).$

  * (4)

If $S_{h}^{[2k-2]}$ is a time-symmetric integrator and the sequence
$(\gamma_{1},\dots,\gamma_{s})$ is palindromic, in the sense that for all $j$

$\gamma_{s-j+1}=\gamma_{j},$

then clearly the composition (58) is time-symmetric, and hence of even order.

####  2.1.2 Composition of Lie-Trotter maps

One could consider an analogous composition to (48), but this time with the
Lie–Trotter scheme $\chi_{h}=\varphi_{h}^{[1]}\circ\varphi_{h}^{[2]}$ as the
basic method. In that case, the Lie transform of $\chi_{h}$ is of the form
$\mathrm{e}^{Y(h)}$ with

$Y(h)=h(F_{1}+F_{2})+h^{2}Y_{2}+h^{3}Y_{3}+\cdots,$ (60)

so that the Lie transform $\Psi(h)$ associated to the composition
$\chi_{\gamma_{s}h}\circ\cdots\circ\chi_{\gamma_{2}h}\circ\chi_{\gamma_{1}h}$
is

$\displaystyle\Psi(h)$
$\displaystyle=\mathrm{e}^{Y(\gamma_{1}h)}\cdots\mathrm{e}^{Y(\gamma_{s}h)}$
$\displaystyle=\mathrm{e}^{\left(\sum_{j=1}^{s}\gamma_{j}\right)(F_{1}+F_{2})}+h^{2}\left(\sum_{j=1}^{s}\gamma_{j}^{2}\right)Y_{2}+\mathcal{O}(h^{3}).$

This shows that such a scheme is of order 2 if

$\sum_{j=1}^{s}\gamma_{j}=1,\qquad\quad\sum_{j=1}^{s}\gamma_{j}^{2}=0.$ (61)

Obviously, such a system of equations does not admit real solutions.

The situation is different, however, if one composes $\chi_{h}$ with its
adjoint $\chi_{h}^{*}=\varphi_{h}^{[2]}\circ\varphi_{h}^{[1]}$, that is,

$\psi_{h}=\chi_{\alpha_{2s}h}\circ\chi^{*}_{\alpha_{2s-1}h}\circ\cdots\circ\chi_{\alpha_{2}h}\circ\chi^{*}_{\alpha_{1}h}.$
(62)

If $\mathrm{e}^{Y(h)}$, with $Y(h)$ given by (60), is the Lie transform of
$\chi_{h}$, then the Lie transform of its adjoint $\chi^{*}_{h}$ is
$\mathrm{e}^{-Y(-h)}$ and the Lie transform $\Psi(h)$ of (62) satisfies

$\displaystyle\Psi(h)$
$\displaystyle=\mathrm{e}^{-Y(-\alpha_{1}h)}\mathrm{e}^{Y(\alpha_{2}h)}\cdots\,\mathrm{e}^{-Y(-\alpha_{2s-1}h)}\mathrm{e}^{Y(\alpha_{2s}h)}$
(63)
$\displaystyle=\mathrm{e}^{\left(\sum_{j=1}^{2s}\alpha_{j}\right)(F_{1}+F_{2})}+h^{2}\left(\sum_{j=1}^{s}\big{(}\alpha_{2j}^{2}-\alpha_{2j-1}^{2}\big{)}\right)Y_{2}+\mathcal{O}(h^{3}),$

whence composition (62) is at least of order 2 if the coefficients satisfy

$\sum_{j=1}^{2s}\alpha_{j}=1,\qquad\quad\sum_{j=1}^{s}\big{(}\alpha_{2j}^{2}-\alpha_{2j-1}^{2}\big{)}=0.$

The argument above also holds when $\chi_{h}$ is a first order integrator
other than Lie-Trotter.

As mentioned earlier, the simplest situation corresponds to $s=1$, in which
case we recover the Strang splitting,
$S_{h}^{[2]}=\chi_{h/2}\circ\chi_{h/2}^{*}$. In fact, the general scheme (62)
can be rewritten as the splitting method

$\psi_{h}=\varphi^{[1]}_{a_{s+1}h}\circ\varphi^{[2]}_{b_{s}h}\circ\varphi^{[1]}_{a_{s}h}\circ\cdots\circ\varphi^{[1]}_{a_{2}h}\circ\varphi^{[2]}_{b_{1}h}\circ\varphi^{[1]}_{a_{1}h},$
(64)

where $a_{1}=\alpha_{1}$, and for $j=1,\ldots,s$,

$a_{j+1}=\alpha_{2j}+\alpha_{2j+1},\qquad\quad
b_{j}=\alpha_{2j-1}+\alpha_{2j}$ (65)

(with $\alpha_{2s+1}=0$). Conversely, any integrator of the form (64)
satisfying the condition $\sum_{j=1}^{s+1}a_{j}=\sum_{j=1}^{s}b_{j}$ can be
expressed in the form (62), as shown in McLachlan (1995b).

Clearly, any splitting scheme (64) with a palindromic sequence of
coefficients, that is, satisfying

$a_{s-j+2}=a_{j},\qquad b_{s-j+1}=b_{j}\qquad\mbox{ for all }j,$

is time-symmetric, and thus of even order. Written in the composition format
(62), it is time-symmetric if

$\alpha_{2s-j+1}=\alpha_{j}.$

###  2.2 Order conditions I: splitting schemes with the BCH formula

In the analysis of the order conditions for splitting methods, and without
lost of generality, we will consider the linear case (4), so that the
treatment is essentially based on matrices, or more generally, on linear
operators. Thus, the integrator (64) corresponds in this setting to the
product of exponentials

$\Psi(h)=\mathrm{e}^{a_{s+1}hF_{1}}\,\mathrm{e}^{b_{s}hF_{2}}\,\mathrm{e}^{a_{s}hF_{1}}\,\cdots\,\mathrm{e}^{a_{2}hF_{1}}\,\mathrm{e}^{b_{1}hF_{2}}\,\mathrm{e}^{a_{1}hF_{1}},$
(66)

intended to approximate $\mathrm{e}^{h(F_{1}+F_{2})}$.

In the case of systems of ODEs $x^{\prime}=f_{1}(x)+f_{2}(x)$ with (non
necessarily linear) vector fields $f_{1}$ and $f_{2}$, one may consider the
Lie operators $F_{j}$ of $f_{j}$, and compare the Lie transformation
$\mathrm{e}^{h(F_{1}+F_{2})}$ with the operator

$\Psi(h)=\mathrm{e}^{a_{1}hF_{1}}\,\mathrm{e}^{b_{1}hF_{2}}\,\mathrm{e}^{a_{2}hF_{1}}\,\cdots\,\mathrm{e}^{a_{s}hF_{1}}\,\mathrm{e}^{b_{s}hF_{2}}\,\mathrm{e}^{a_{s+1}hF_{1}}$
(67)

associated with the map $\psi_{h}$ (in the sense that
$g(\psi_{h}(x))=\Psi(h)g(x)$ for arbitrary smooth $g\in
C^{\infty}(\mathbb{R}^{D},\mathbb{R})$). Thus, all the formulas derived below
for the linear case can be applied to the Lie transformation (67) associated
with the map (64) just by reversing the order of the sequence
$(a_{1},b_{1},a_{2},\ldots,a_{s},b_{s},a_{s+1})$ of coefficients of the
splitting scheme. Alternatively, we may reverse the order of the products of
the operators in all the expressions involved.

Equivalent results could also be obtained for the ODE case using the concept
of word series and related techniques introduced in Murua and Sanz-Serna
(2017). The Lie transform approach is more general, however, as it can be
directly applied to the case of splitting methods for differential equations
on manifolds, as discussed at the end of Subsection 4.1.

Generally speaking, the order conditions for a method of order $r$ are systems
of polynomial equations in the coefficients obtained by requiring that the
Taylor expansions in the step size $h$ of both the exact and numerical
solution agree up to terms in $h^{r}$. A standard approach to obtain the order
conditions of scheme (66) consists in formally using the
Baker–Cambpell–Hausdorff (BCH) formula to express $\Psi(h)$ as one exponential
of a series of operators in powers of $h$, and finally to compare this series
with $\mathrm{e}^{h(F_{1}+F_{2})}$. In this way, one gets

$\displaystyle\log(\Psi(h))=h(w_{1}F_{1}+w_{2}F_{2})+h^{2}w_{12}F_{12}+h^{3}(w_{122}F_{122}+w_{112}F_{112})$
(68)
$\displaystyle\qquad+h^{4}(w_{1222}F_{1222}+w_{1122}F_{1122}+w_{1112}F_{1112})+\mathcal{O}(h^{5}),$

where

$\displaystyle F_{12}=[F_{1},F_{2}],\qquad F_{122}=[F_{12},F_{2}],\qquad
F_{112}=[F_{1},F_{12}]$ $\displaystyle F_{1222}=[F_{122},F_{2}],\qquad
F_{1122}=[F_{1},F_{122}],\qquad F_{1112}=[F_{1},F_{112}]$

and $w_{1}$, $w_{2}$, $w_{12}$, $w_{122},\ldots$ are polynomials of
homogeneous degree in the parameters $a_{j}$, $b_{j}$. In particular

$\displaystyle w_{1}=\sum_{i=1}^{s+1}a_{i},\qquad
w_{2}=\sum_{i=1}^{s}b_{i},\qquad w_{12}=\frac{1}{2}w_{1}w_{2}-\sum_{1\leq
i<j\leq s+1}b_{i}a_{j},$ $\displaystyle
w_{122}=\frac{1}{12}w_{1}w_{2}^{2}-\frac{1}{2}\sum_{1\leq i\leq j<k\leq
s}b_{i}a_{j}b_{k},$ $\displaystyle
w_{112}=\frac{1}{12}w_{1}^{2}w_{2}-\frac{1}{2}\sum_{1\leq i<j\leq k\leq
s+1}a_{i}b_{j}a_{k},$

with $b_{s+1}=0$. From (68), it is clear that a characterization of the order
of the splitting scheme (64) is obtained by requiring the consistency
conditions $w_{1}=w_{2}=1$, that is,

$\sum_{j=1}^{s+1}a_{j}=\sum_{j=1}^{s}b_{j}=1$ (69)

(ensuring that the scheme (64) is at least of order 1) and
$w_{12}=w_{122}=w_{112}=\cdots=0$ up to the required order. The set of order
conditions thus obtained will be independent in the general case if the
operators $F_{1}$, $F_{2}$, $F_{12}$, $F_{122}$, $F_{112},\ldots$ form a basis
of the free Lie algebra in the alphabet $\\{1,2\\}$. We have considered in
(68) the so-called _Lyndon basis_ , associated to the set of Lyndon words in
the alphabet $\\{1,2\\}$

$\\{1,2,12,122,112,1222,1122,1112,\ldots\\}.$ (70)

They are defined as follows Reutenauer (1993): a word $\ell_{1}\cdots\ell_{m}$
is a Lyndon word if $(\ell_{1}\cdots\ell_{k})\prec(\ell_{k+1}\cdots\ell_{m})$
for each $1\leq k<m$, where $\prec$ is the lexicographical order (i.e., the
order used when ordering words in the dictionary) on the set of words on the
alphabet $\\{1,2\\}$. For instance, $112$ is a Lyndon word, while neither
$211$ nor $121$ are, as $2\not\prec 11$ and as $12\not\prec 1$, respectively.

The element of the basis associated to a Lyndon word $\ell_{1}\cdots\ell_{m}$
with $m\geq 2$ is given as
$F_{\ell_{1}\cdots\ell_{m}}=[F_{\ell_{1}\cdots\ell_{n}},F_{\ell_{n+1}\cdots\ell_{m}}]$,
where $n$ is the smallest number such that both $\ell_{1}\cdots\ell_{n}$ and
$\ell_{n+1}\cdots\ell_{m}$ are themselves Lyndon words.

An efficient algorithm (based on the results in Murua (2006)) for the BCH
formula and related calculations in the Lyndon basis (and some other basis)
that allows one to obtain (68) up to terms of arbitrarily high degree is
presented in Casas and Murua (2009).

Of course, if another basis of the free Lie algebra in the alphabet
$\\{1,2\\}$ is used to expand $\log(\Psi(h))$ in (68), a different
characterization of the order conditions will be obtained, with a different
set of polynomial functions on $(a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1})$. In
any case, the number of such independent conditions arising at each order $n$
can be obtained just by determining the dimension of
$\mathcal{L}_{n}(F_{1},F_{2})$, the linear span of all commutators containing
$n$ operators $F_{1}$, $F_{2}$. This number, denoted $c_{n}$, is given in
Table 2.1 (see Munthe-Kaas and Owren (1999), McLachlan and Quispel (2002)).

Order $n$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11
---|---|---|---|---|---|---|---|---|---|---|---
$c_{n}$ | 2 | 1 | 2 | 3 | 6 | 9 | 18 | 30 | 56 | 99 | 186
$d_{n}$ | 2 | 1 | 2 | 2 | 4 | 5 | 10 | 14 | 25 | 39 | 69
$m_{n}$ | 1 | 0 | 1 | 1 | 2 | 2 | 4 | 5 | 8 | 11 | 17
Table 2.1: Number of independent order conditions for general splitting
methods, $c_{n}$, and for RKN-type splitting methods, $d_{n}$. The number
$m_{n}$ corresponds to the number of order conditions for compositions of a
2nd-order time-symmetric method.

###  2.3 Order conditions II: splitting schemes with Lyndon words

Whereas the previous characterization of the order of the splitting scheme
(66) allows one to easily get the number of order conditions, obtaining
explicit expressions for the polynomials $w_{\ell_{1}\cdots\ell_{n}}$ is much
more difficult when the considered order increases. The following alternative
characterization, based on the direct comparison of the power series
expansions of $\mathrm{e}^{h\,(F_{1}+F_{2})}$ and (66), tries to ameliorate
this difficulty.

####  2.3.1 Basic expansions and necessary order conditions

Both $\mathrm{e}^{h\,(F_{1}+F_{2})}$ and (66) admit an expansion in series
indexed by the set

$\mathcal{W}=\\{1,2,11,12,21,22,111,112,121,211,122,\cdots\\}$

of words in the alphabet $\\{1,2\\}$. More precisely,
$\mathrm{e}^{h\,(F_{1}+F_{2})}$ can be expanded as

$I+h\alpha_{1}F_{1}+h\alpha_{2}F_{2}+h^{2}\alpha_{11}F_{1}F_{1}+h^{2}\alpha_{12}F_{1}F_{2}+h^{2}\alpha_{21}F_{2}F_{1}+h^{2}\alpha_{22}F_{2}F_{2}+\cdots,$
(71)

with $\alpha_{\ell_{1}\cdots\ell_{n}}=1/n!$. As for (66), it can be expanded,
for arbitrary $s$, as (71), where for each word
$\ell_{1}\cdots\ell_{n}\in\mathcal{W}$ with $n$ letters, the corresponding
coefficient

$\alpha_{\ell_{1}\cdots\ell_{n}}=u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1})$

is a homogeneous polynomial of degree $n$ in the variables
$a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1}$.

It is straightforward to check that such polynomials satisfy the following
relations, which allows one to compute them recursively:

  * •

if $\ell_{1}=\cdots=\ell_{j}=1$ and $\ell_{j+1}\neq 1$ with $j\geq 1$,

$\displaystyle
u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1})$
$\displaystyle=\sum_{k=0}^{j}u_{\ell_{k+1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},0)\,\frac{a_{s+1}^{k}}{k!},$
$\displaystyle u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},0)$
$\displaystyle=u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s}).$

  * •

if $\ell_{1}=\cdots=\ell_{j}=2$ and $\ell_{j+1}\neq 2$ with $j\geq 1$,

$u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1})=\sum_{k=0}^{j}u_{\ell_{k+1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s})\,\frac{b_{s}^{k}}{k!},$

  * •

if $\ell_{n}=\cdots=\ell_{1}=1$,

$u_{\ell_{1}\cdots\ell_{n}}(a_{1})=\frac{a_{1}^{n}}{n!},$

  * •

if $\ell_{1}\neq 1$,

$u_{\ell_{1}\cdots\ell_{n}}(a_{1})=0.$

In this way (66) is at least of order $r$ if and only if the conditions

$u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1})=\frac{1}{n!}$
(72)

hold for each word $\ell_{1}\cdots\ell_{n}$ with $n\leq r$ letters in the
alphabet $\\{1,2\\}$. For illustration, in Table 2.2 we explicitly give these
conditions corresponding to words with up to two letters. Notice that (72) for
the single-letter words $1$ and $2$ coincide with the consistency conditions
(69).

Word |  Condition | Word |  Condition
---|---|---|---
$1$ | $\displaystyle\sum_{j=1}^{s+1}a_{j}=1$ | $2$ | $\displaystyle\sum_{j=1}^{s}b_{j}=1$
$11$ | $\displaystyle\frac{1}{2}\sum_{j=1}^{s+1}\frac{a_{j}^{2}}{2}+\sum_{1\leq i<j\leq s+1}a_{i}a_{j}=\frac{1}{2}$ | $12$ | $\displaystyle\sum_{1\leq i<j\leq s+1}b_{i}a_{j}=\frac{1}{2}$
$21$ | $\displaystyle\sum_{1\leq j\leq i\leq s}b_{i}a_{j}=\frac{1}{2}$ | $22$ | $\displaystyle\frac{1}{2}\sum_{j=1}^{s}\frac{b_{j}^{2}}{2}+\sum_{1\leq j<i\leq s}b_{j}b_{i}=\frac{1}{2}$
Table 2.2: Conditions (72) corresponding to words with up to two indices.

However, such order conditions are not all independent. For instance, from
Table 2.2, one can check that

$u_{1}^{2}=2u_{11},\quad u_{2}^{2}=2u_{22},\quad u_{1}u_{2}=u_{12}+u_{21}.$

For a consistent method, $u_{1}=u_{2}=1$, hence $u_{11}=\frac{1}{2}$,
$u_{22}=\frac{1}{2}$, and

$u_{12}-\frac{1}{2}=\frac{1}{2}-u_{21},$

which implies that, if $u_{12}=\frac{1}{2}$, then automatically
$u_{21}=\frac{1}{2}$.

A complete characterization of the relations among the order conditions (72)
will be obtained in paragraph 2.3.3 below. As a previous step, we obtain
integral representations of both $\mathrm{e}^{h\,(F_{1}+F_{2})}$ and (66) ,
which in addition give useful expressions for the remainders of their
truncated series expansions.

####  2.3.2 Integral representation and remainders

Consider the solution $Y(\tau)$ of the initial value problem

$\frac{d}{d\tau}Y(\tau)=hA(\tau)Y(\tau),\qquad Y(0)=I,$ (73)

with $A(\tau)=d_{1}(\tau)F_{1}+d_{2}(\tau)F_{2}$ and

$(d_{1}(\tau),d_{2}(\tau))=\left\\{\begin{matrix}(a_{i},0)&\mbox{ if
}&\tau\in[2i-2,2i-1],\ i\in\\{1,\ldots,s+1\\},\\\ (0,b_{i})&\mbox{ if
}&\tau\in[2i-1,2i],\ i\in\\{1,\ldots,s\\},\\\ (0,0)&\mbox{ if
}&\tau>2s+1.\end{matrix}\right.$ (74)

It is straightforward to check that

$Y(\tau)=\left\\{\begin{matrix}\mathrm{e}^{(\tau-(2i-2))ha_{i+1}F_{1}}\,\mathrm{e}^{b_{i}hF_{2}}\,\mathrm{e}^{a_{i}hF_{1}}\,\cdots\,\mathrm{e}^{a_{2}hF_{1}}\,\mathrm{e}^{b_{1}hF_{2}}\,\mathrm{e}^{a_{1}hF_{1}},&\tau\in[2i-2,2i-1],\\\
\mathrm{e}^{(\tau-(2i-1))hb_{i}F_{2}}\,\mathrm{e}^{a_{i}hF_{1}}\,\cdots\,\mathrm{e}^{a_{2}hF_{1}}\,\mathrm{e}^{b_{1}hF_{2}}\,\mathrm{e}^{a_{1}hF_{1}},&\tau\in[2i-1,2i].\end{matrix}\right.$

In particular, $\Psi(h)=Y(2s+1)$. The solution $Y(\tau)$ of (73) satisfies

$Y(\tau)=I+h\int_{0}^{\tau}A(\tau_{1})Y(\tau_{1})d\tau_{1}.$ (75)

From that, one obtains

$Y(\tau)=I+h\int_{0}^{\tau}A(\tau_{1})d\tau_{1}+h^{2}\int_{0}^{\tau}\int_{0}^{\tau_{1}}A(\tau_{1})A(\tau_{2})Y(\tau_{2})d\tau_{2}d\tau_{1}$

and more generally,

$Y(\tau)=I+\sum_{n=1}^{m}h^{n}\int_{0}^{\tau}\int_{0}^{\tau_{1}}\cdots\int_{0}^{\tau_{n-1}}A(\tau_{1})\cdots
A(\tau_{n})d\tau_{n}\cdots d\tau_{1}+h^{m+1}\mathcal{R}_{m+1}(\tau,h),$ (76)

where for each $n$, the remainder $\mathcal{R}_{n}(\tau,h)$ satisfies

$\mathcal{R}_{n}(\tau,h)=\int_{0}^{\tau}\int_{0}^{\tau_{1}}\cdots\int_{0}^{\tau_{n-1}}A(\tau_{1})\cdots
A(\tau_{n})Y(\tau_{n})d\tau_{n}\cdots d\tau_{1}.$ (77)

By substituting $A(\tau_{j})=d_{1}(\tau_{j})F_{1}+d_{2}(\tau_{j})F_{2}$, then
expanding all the products and taking constant linear operators out from
integral signs, we obtain

$Y(\tau)=I+\sum_{n=1}^{m}h^{n}\sum_{\ell_{1},\ldots,\ell_{n}\in\\{1,2\\}}\alpha_{\ell_{1}\cdots\ell_{n}}(\tau)\,F_{\ell_{1}}\cdots
F_{\ell_{n}}+h^{m+1}\mathcal{R}_{m+1}(\tau,h),$

where

$\alpha_{\ell_{1}\cdots\ell_{n}}(\tau)=\int_{0}^{\tau}\int_{0}^{\tau_{1}}\cdots\int_{0}^{\tau_{n-1}}d_{\ell_{1}}(\tau_{1})\cdots
d_{\ell_{n}}(\tau_{n})d\tau_{n}\cdots d\tau_{1}.$ (78)

In particular, we have

$\begin{split}\Psi(h)&=I+\sum_{n=1}^{m}h^{n}\sum_{\ell_{1},\ldots,\ell_{n}\in\\{1,2\\}}u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1})\,F_{\ell_{1}}\cdots
F_{\ell_{n}}\\\ &+h^{m+1}\mathcal{R}_{m+1}(2s+1,h),\end{split}$

where
$u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1})=\alpha_{\ell_{1}\cdots\ell_{n}}(2s+1)$,
that is,

$\displaystyle
u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1})=$
$\displaystyle\qquad\qquad\int_{0}^{2s+1}\int_{0}^{\tau_{1}}\cdots\int_{0}^{\tau_{n-1}}d_{\ell_{1}}(\tau_{1})\cdots
d_{\ell_{n}}(\tau_{n})d\tau_{n}\cdots d\tau_{1}.$

We now consider (73) with $d_{1}(\tau)\equiv d_{2}(\tau)\equiv 1$. Clearly, in
that case, $Y(\tau)=\mathrm{e}^{\tau h(F_{1}+F_{2})}$, and

$\displaystyle\mathrm{e}^{h(F_{1}+F_{2})}$
$\displaystyle=I+\sum_{n=1}^{m}h^{n}\int_{0}^{1}\int_{0}^{\tau_{1}}\cdots\int_{0}^{\tau_{n-1}}(F_{1}+F_{2})^{n}d\tau_{n}\cdots
d\tau_{1}+h^{m+1}\overline{\mathcal{R}}_{m+1}(1,h),$
$\displaystyle=I+\sum_{n=1}^{m}h^{n}\sum_{\ell_{1},\ldots,\ell_{n}\in\\{1,2\\}}\frac{1}{n!}\,F_{\ell_{1}}\cdots
F_{\ell_{n}}+h^{m+1}\overline{\mathcal{R}}_{m+1}(1,h),$

where for each $n$, the remainder $\overline{\mathcal{R}}_{n}(\tau,h)$ is
given by

$\overline{\mathcal{R}}_{n}(\tau,h)=\int_{0}^{\tau}\int_{0}^{\tau_{1}}\cdots\int_{0}^{\tau_{n-1}}(F_{1}+F_{2})^{n}\mathrm{e}^{\tau_{n}h(F_{1}+F_{2})}d\tau_{n}\cdots
d\tau_{1}.$

We finally arrive at the following expression for the local error of the
splitting scheme (73):

$\displaystyle\Psi(h)-\mathrm{e}^{h(F_{1}+F_{2})}=$
$\displaystyle\quad\sum_{n=1}^{m}h^{n}\sum_{\ell_{1},\ldots,\ell_{n}\in\\{1,2\\}}\left(u_{\ell_{1}\cdots\ell_{n}}(a_{1},b_{1},\ldots,a_{s},b_{s},a_{s+1})-\frac{1}{n!}\right)\,F_{\ell_{1}}\cdots
F_{\ell_{n}}$
$\displaystyle\quad+h^{m+1}\left(\mathcal{R}_{m+1}(2s+1,h)-\overline{\mathcal{R}}_{m+1}(1,h)\right).$

Hence, if the scheme is of order $r$ (that is, (72) holds for each word
$(\ell_{1},\ldots,\ell_{n})$ with $n\leq r$ letters in the alphabet
$\\{1,2\\}$), then

$\Psi(h)-\mathrm{e}^{h(F_{1}+F_{2})}=h^{n+1}\left(\mathcal{R}_{n+1}(2s+1,h)-\overline{\mathcal{R}}_{n}(1,h)\right).$

####  2.3.3 Iterated integrals and shuffle relations

Iterated integrals of the form (78) were first considered and studied in Chen
(1957). It is well known that the integration-by-parts formula gives (for an
arbitrary integrable path $(d_{1}(\tau),d_{2}(\tau))$) the relations

$\displaystyle\alpha_{\ell_{1}}(\tau)\alpha_{\ell_{2}}(\tau)$
$\displaystyle=\alpha_{\ell_{1}\ell_{2}}(\tau)+\alpha_{\ell_{2}\ell_{1}}(\tau),$
$\displaystyle\alpha_{\ell_{1}}(\tau)\alpha_{\ell_{2}\ell_{3}}(\tau)$
$\displaystyle=\alpha_{\ell_{1}\ell_{2}\ell_{3}}(\tau)+\alpha_{\ell_{2}\ell_{1}\ell_{3}}(\tau)+\alpha_{\ell_{2}\ell_{3}\ell_{1}}(\tau),$

and more generally,

$\alpha_{\ell_{1}\cdots\ell_{n}}(\tau)\alpha_{\ell_{n+1}\cdots\ell_{n+m}}(\tau)=\sum_{\sigma\in\mathrm{Sh}(n,m)}\alpha_{\ell_{\sigma(1)}\cdots\ell_{\sigma(n+m)}}(\tau),$
(79)

where $\mathrm{Sh}(n,m)$ is the set of the $(n+m)!/(n!m!)$ permutations
$\sigma$ of $(1,\ldots,n+m)$ that are obtained by interleaving $(1,\ldots,n)$
and $(n+1,\ldots,n+m)$ while preserving their respective ordering.

It will be useful to interpret the relations (79) in terms of the so-called
shuffle product of words: the shuffle product $\shuffle$ of two words
$\ell_{1}\cdots\ell_{n}$ and $\ell_{n+1}\cdots\ell_{n+m}$ is defined as the
following formal sum of words,

$\ell_{1}\cdots\ell_{n}\shuffle\ell_{n+1}\cdots\ell_{n+m}=\sum_{\sigma\in\mathrm{Sh}(n,m)}\ell_{\sigma(1)}\cdots\ell_{\sigma(n+m)}.$

By extending $\alpha_{w}(\tau)$ linearly to the case where $w$ is a linear
combination of words, the relations (79) can be interpreted as

$\alpha_{w}(\tau)\alpha_{w^{\prime}}(\tau)=\alpha_{w\shuffle w^{\prime}}$

for arbitrary words $w=\ell_{1}\cdots\ell_{n}$ and
$w^{\prime}=\ell_{n+1}\cdots\ell_{n+m}$ in the alphabet
$\mathcal{A}=\\{1,2\\}$.

The shuffle product $\shuffle$ defines a commutative algebra (the so-called
shuffle algebra) over the vector space of formal linear combinations of words
in the alphabet $\mathcal{A}$. The shuffle algebra is freely generated by the
set of Lyndon words Reutenauer (1993).

The fact that the coefficients of the series expansions of both
$\mathrm{e}^{h(F_{1}+F_{2})}$ and $\Psi(h)$ satisfy the shuffle relations,
together with the fact that the set of Lyndon words freely generate the
shuffle algebra, implies that a set of independent conditions for a consistent
splitting scheme to attain order $r$ can be obtained by considering (72) for
each Lyndon word $(\ell_{1},\ldots,\ell_{n})$ of length $n\leq r$.

###  2.4 Order conditions III: splitting methods with Lyndon multi-indices

In Blanes et al. (2013b), yet another characterization of the order conditions
is obtained in terms of explicitly given polynomial equations. We next
describe this alternative formulation. To do that, we always assume that the
consistency conditions (69) hold, so that the method is at least of first
order. In this case, the polynomial equations are expressed in terms of the
coefficients $b_{1},\ldots,b_{s}$ and the coefficients $c_{1},\ldots,c_{s}$
given by

$c_{i}=\sum_{j=1}^{i}a_{j},\qquad i=1,2,\ldots,s.$ (80)

We begin by rewriting (66) as

$\displaystyle\Psi(h)=\mathrm{e}^{a_{s+1}hF_{1}}\,\mathrm{e}^{b_{s}hF_{2}}\,\mathrm{e}^{a_{s}hF_{1}}\cdots\mathrm{e}^{b_{1}hF_{2}}\,\mathrm{e}^{a_{1}hF_{1}}$
$\displaystyle\qquad=\mathrm{e}^{hF_{1}}\,\left(\prod_{j=s}^{1}\mathrm{e}^{-c_{j}hF_{1}}\,\mathrm{e}^{b_{j}hF_{2}}\,\mathrm{e}^{c_{j}hF_{1}}\right)=\mathrm{e}^{hF_{1}}\mathrm{e}^{b_{s}hC(c_{s}h)}\cdots\mathrm{e}^{b_{1}hC(c_{1}h)},$

where

$C(h)=\mathrm{e}^{-hF_{1}}\,F_{2}\,\mathrm{e}^{hF_{1}}=\sum_{n=1}^{\infty}h^{n-1}C_{n}=C_{1}+hC_{2}+h^{2}C_{3}+h^{3}C_{4}+\cdots,$
(81)

with $C_{1}=F_{2}$, and

$C_{k}=\frac{1}{(k-1)}[C_{k-1},F_{1}]\qquad\mbox{ for }\,k>1.$

Now the order of the scheme (64) is established by comparing the expansion in
powers of $h$ of $\mathrm{e}^{-hF_{1}}\Psi(h)$ with that of
$\mathrm{e}^{-hF_{1}}\mathrm{e}^{h(F_{1}+F_{2})}$.

Clearly, $Y(\tau):=\mathrm{e}^{-h\tau F_{1}}\mathrm{e}^{\tau h(F_{1}+F_{2})}$
is the solution of (73) with $A(\tau)=C(\tau h)$. Since the solution $Y(\tau)$
of (73) admits the representation (76) with remainder (77), and
$\mathrm{e}^{-hF_{1}}\mathrm{e}^{h(F_{1}+F_{2})}=Y(1)$ we conclude that

$\begin{split}\mathrm{e}^{-hF_{1}}\mathrm{e}^{h(F_{1}+F_{2})}&=I+\sum_{k=1}^{n}h^{k}\int_{0}^{1}\int_{0}^{\tau_{k}}\cdots\int_{0}^{\tau_{2}}C(\tau_{k}h)\cdots
C(\tau_{1}h)d\tau_{1}\cdots d\tau_{k}\\\
&+h^{n+1}\mathcal{R}_{n+1}(1,h),\end{split}$

where for each $k$,

$\mathcal{R}_{k}(1,h)=\int_{0}^{1}\int_{0}^{\tau_{k}}\cdots\int_{0}^{\tau_{2}}C(\tau_{k}h)\cdots
C(\tau_{1}h)Y(\tau_{1})d\tau_{1}\cdots d\tau_{k}.$

By substitution of (81), we obtain that for each $k$,

$\displaystyle
h^{k}\int_{0}^{1}\int_{0}^{\tau_{k}}\cdots\int_{0}^{\tau_{2}}C(\tau_{k}h)\cdots
C(\tau_{1}h)d\tau_{1}\cdots d\tau_{k}$
$\displaystyle=\sum_{i_{1},\ldots,i_{k}\geq
1}h^{i_{1}+\cdots+i_{k}}\left(\int_{0}^{1}\int_{0}^{\tau_{k}}\cdots\int_{0}^{\tau_{2}}\tau_{1}^{i_{1}-1}\cdots\tau_{k}^{i_{k}-1}d\tau_{1}\cdots
d\tau_{k}\right)\,C_{i_{k}}\cdots C_{i_{1}}$
$\displaystyle=\sum_{i_{1},\ldots,i_{k}\geq
1}\frac{h^{i_{1}+\cdots+i_{k}}}{(i_{1}+\cdots+i_{k})\cdots(i_{1}+i_{2})i_{1}}\,C_{i_{k}}\cdots
C_{i_{1}}.$

As for $\mathrm{e}^{-hF_{1}}\,\Psi(h)$,

$\displaystyle\mathrm{e}^{-hF_{1}}\,\Psi(h)=\prod_{j=s}^{1}\mathrm{e}^{b_{j}hC(c_{j}h)}=\prod_{j=s}^{1}\left(I+\sum_{k\geq
1}\frac{h^{k}b_{j}^{k}}{k!}C(c_{j}h)^{k}\right)$ (82)
$\displaystyle\quad=I+\sum_{k\geq
1}\frac{h^{k}}{k!}\sum_{j=1}^{s}b_{j}^{k}C(c_{j}h)^{k}+\sum_{k_{1},k_{2}\geq
1}\frac{h^{k_{1}+k_{2}}}{k_{1}!k_{2}!}\sum_{1\leq j_{1}<j_{2}\leq
s}b_{j_{1}}^{k_{1}}b_{j_{2}}^{k_{2}}C(c_{j_{2}}h)^{k_{2}}C(c_{j_{1}}h)^{k_{1}}$
$\displaystyle\qquad+\sum_{k_{1},k_{2},k_{3}\geq
1}\frac{h^{k_{1}+k_{2}+k_{3}}}{k_{1}!k_{2}!k_{3}!}\sum_{1\leq
j_{1}<j_{2}<j_{3}\leq
s}b_{j_{1}}^{k_{1}}b_{j_{2}}^{k_{2}}b_{j_{3}}^{k_{3}}C(c_{j_{3}}h)^{k_{3}}C(c_{j_{2}}h)^{k_{2}}C(c_{j_{1}}h)^{k_{1}}+\cdots$
$\displaystyle\quad=I+\sum_{k\geq 1}h^{k}\sum_{1\leq j_{1}\leq\cdots\leq
j_{k}\leq s}\frac{b_{j_{1}}\cdots b_{j_{k}}}{\sigma(j_{1},\ldots
j_{k})}\,C(c_{j_{k}}h)\cdots C(c_{j_{1}}h)$

where

$\displaystyle\sigma(j_{1},\ldots,j_{k})=1$ if $\displaystyle
j_{1}<\cdots<j_{k},$
$\displaystyle\sigma(j_{1},\ldots,j_{k})=\ell!\,\sigma(j_{\ell+1},\ldots,j_{k})$
if $\displaystyle j_{1}=\cdots=j_{\ell}<j_{\ell+1}\leq\cdots\leq j_{k}.$

Since

$h^{k}C(c_{j_{k}}h)\cdots C(c_{j_{1}}h)=\sum_{i_{1},\ldots,i_{k}\geq
1}h^{i_{1}+\cdots+i_{k}}c_{j_{1}}^{i_{1}-1}\cdots
c_{j_{k}}^{i_{k}-1}\,C_{i_{k}}\cdots C_{i_{1}},$

we arrive at

$\mathrm{e}^{-hF_{1}}\,\Psi(h)=I+\sum_{k\geq 1}\sum_{i_{1},\ldots,i_{k}\geq
1}h^{i_{1}+\cdots+i_{k}}\,v_{i_{1},\ldots,i_{k}}(b_{1},c_{1},\ldots,b_{s},c_{s})\,C_{i_{k}}\cdots
C_{i_{1}},$ (83)

where

$v_{i_{1},\ldots,i_{k}}(b_{1},c_{1},\ldots,b_{s},c_{s})=\sum_{1\leq
j_{1}\leq\cdots\leq j_{k}\leq s}\frac{b_{j_{1}}\cdots
b_{j_{k}}}{\sigma(j_{1},\ldots,j_{k})}\,c_{j_{1}}^{i_{1}-1}\cdots
c_{j_{k}}^{i_{k}-1}.$ (84)

In this way, a consistent splitting method is at least of order $r$ if and
only if

$\sum_{1\leq j_{1}\leq\cdots\leq j_{k}\leq s}\frac{b_{j_{1}}\cdots
b_{j_{k}}}{\sigma(j_{1},\ldots,j_{k})}\,c_{j_{1}}^{i_{1}-1}\cdots
c_{j_{k}}^{i_{k}-1}=\frac{1}{(i_{1}+\cdots+i_{k})\cdots(i_{1}+i_{2})i_{1}}$
(85)

holds for each multi-index $(i_{1},\ldots,i_{k})$ such that
$i_{1}+\cdots+i_{k}\leq r$. For illustration, in Table 2.3 we give explicit
conditions (85) corresponding to multi-indices with up to three indices.

However, such order conditions are not independent. The situation is very
similar to that of the previous subsection: instead of series indexed by the
set of words in the alphabet $\\{1,2\\}$, now they are indexed by the set of
words in the alphabet $\mathbb{N}=\\{1,2,3,4,\ldots\\}$. To distinguish the
words of both sets, we will keep referring to the words in the alphabet
$\mathbb{N}$ as multi-indices, and will write them as $(i_{1},\ldots,i_{k})$
instead of $i_{1}\cdots i_{k}$. Analogously to the previous subsection, the
corresponding coefficients of the series expansions of both
$\mathrm{e}^{-hF_{1}}\Psi(h)$ and
$\mathrm{e}^{-hF_{1}}\mathrm{e}^{h(F_{1}+F_{2})}$ satisfy the shuffle
relations

$v_{i_{1},\ldots,i_{n}}v_{i_{n+1},\ldots,i_{n+m}}=\sum_{\sigma\in\mathrm{Sh}(n,m)}v_{i_{\sigma(1)},\dots,i_{\sigma(n+m)}}.$
(86)

This can be seen by showing that such coefficients can be written in both
cases as iterated integrals. From the discussion above, we already know that,
for each multi-index $(i_{1},\ldots,i_{k})$, the coefficients

$\displaystyle\frac{1}{(i_{1}+\cdots+i_{k})\cdots(i_{1}+i_{2})i_{1}}=\int_{0}^{1}\int_{0}^{\tau_{k}}\cdots\int_{0}^{\tau_{2}}\tau_{1}^{i_{1}-1}\cdots\tau_{k}^{i_{k}-1}d\tau_{1}\cdots
d\tau_{k}$

of the series expansion of $\mathrm{e}^{-hF_{1}}\mathrm{e}^{h(F_{1}+F_{2})}$
are indeed iterated integrals.

We next show that the coefficients (84) of the series expansion (83) can also
be defined as iterated integrals.

We begin by showing that the infinite series expansion
(LABEL:eq:Psi(h)_expansion) of $\mathrm{e}^{-hF_{1}}\,\Psi(h)$ can be
represented as a truncated series plus a remainder. For that purpose, we
consider (73) with $A(\tau)$ defined as follows:

$A(\tau)=\left\\{\begin{matrix}b_{i}C(hc_{i})&\mbox{ if }&\tau\in[i-1,i],\
i\in\\{1,\ldots,s\\},\\\ 0&\mbox{ if }&\tau>s.\end{matrix}\right.$ (87)

In that case, if $\tau\in[i,i+1]$ with $i<s$,

$Y(\tau)=\mathrm{e}^{(\tau-i)hb_{i+1}C(hc_{i+1})}\mathrm{e}^{hb_{i}C(hc_{i})}\cdots\mathrm{e}^{hb_{1}C(hc_{1})}$

and in particular, $\mathrm{e}^{-hF_{1}}\,\Psi(h)=Y(s)$. Since the solution
$Y(\tau)$ of (73) admits the representation (76) with remainder (77), then

$\begin{split}\mathrm{e}^{-hF_{1}}\Psi(h)&=I+\sum_{k=1}^{n}h^{k}\int_{0}^{s}\int_{0}^{\tau_{k}}\cdots\int_{0}^{\tau_{2}}A(\tau_{k})A(\tau_{k-1})\cdots
A(\tau_{1})d\tau_{1}\cdots d\tau_{k}\\\
&+h^{n+1}\mathcal{R}_{n+1}(s,h),\end{split}$

where for each $n$,

$\mathcal{R}_{n}(s,h)=\int_{0}^{s}\int_{0}^{\tau_{n}}\cdots\int_{0}^{\tau_{2}}A(\tau_{n})A(\tau_{n-1})\cdots
A(\tau_{1})Y(\tau_{1})d\tau_{1}\cdots d\tau_{n}.$

By comparison with (LABEL:eq:Psi(h)_expansion), we conclude that

$\displaystyle\sum_{1\leq j_{1}\leq\cdots\leq j_{k}\leq
s}\frac{b_{j_{1}}\cdots b_{j_{k}}}{\sigma(j_{1},\ldots
j_{k})}\,C(c_{j_{k}}h)\cdots C(c_{j_{1}}h)$
$\displaystyle=\int_{0}^{s}\int_{0}^{\tau_{k}}\cdots\int_{0}^{\tau_{2}}A(\tau_{k})A(\tau_{k-1})\cdots
A(\tau_{1})d\tau_{1}\cdots d\tau_{k}$

Now, $A(\tau)$ can be written as

$A(\tau)=d_{1}(\tau)C_{1}+hd_{2}(\tau)C_{2}+h^{2}d_{3}(\tau)C_{3}+\cdots$

where $d_{j}(\tau)=b_{i}c_{i}^{j-1}$ if $\tau\in[i,i+1]$ with $i<s$.
Proceeding as in previous section, one obtains

$v_{i_{1},\ldots,i_{k}}(b_{1},c_{1},\ldots,b_{s},c_{s})=\int_{0}^{s}\int_{0}^{\tau_{1}}\cdots\int_{0}^{\tau_{k-1}}d_{i_{1}}(\tau_{1})\cdots
d_{i_{k}}(\tau_{k})d\tau_{k}\cdots d\tau_{1},$ (88)

which implies that the shuffle relations (86) hold for the polynomials
$v_{i_{1},\ldots,i_{k}}$.

As in the previous subsection, a set of independent conditions that imply the
order conditions (85) can be obtained by considering (85) for each Lyndon
multi-index $(i_{1},\ldots,i_{k})$ such that $1<i_{1}+\cdots+i_{k}\leq r$.
Here, we exclude the multi-index $(1)$ as in that case (85) coincides with the
second equality in the consistency condition (69). For instance, the subset of
Lyndon multi-indices $(i_{1},\ldots,i_{k})$ such that
$1<i_{1}+\cdots+i_{k}\leq 5$ is

$\\{(2),(3),(4),(5),(1,2),(1,3),(1,4),(2,3),(1,1,2),(1,1,3),(1,2,2),(1,1,1,2)\\}.$

Multi-index |  Condition
---|---
$(i)$ | $\displaystyle\sum_{j=1}^{s}b_{j}\,c_{j}^{i-1}=\frac{1}{i}$
$(i_{1},i_{2})$ | $\displaystyle\frac{1}{2}\sum_{j=1}^{s}b_{j}^{2}c_{j}^{i_{1}+i_{2}-2}+\sum_{1\leq j_{1}<j_{2}\leq s}b_{j_{1}}b_{j_{2}}c_{j_{1}}^{i_{1}-1}c_{j_{2}}^{i_{2}-1}=\frac{1}{(i_{1}+i_{2})i_{1}}$
$(i_{1},i_{2},i_{3})$ | $\displaystyle\frac{1}{6}\sum_{j=1}^{s}b_{j}^{3}c_{j}^{i_{1}+i_{2}+i_{3}-3}+\frac{1}{2}\sum_{1\leq j_{1}<j_{3}\leq s}b_{j_{1}}^{2}b_{j_{3}}c_{j_{1}}^{i_{1}+i_{2}-2}c_{j_{3}}^{i_{3}-1}$
| $\displaystyle+\frac{1}{2}\sum_{1\leq j_{1}<j_{2}\leq
s}b_{j_{1}}b_{j_{2}}^{2}c_{j_{1}}^{i_{1}-1}c_{j_{2}}^{i_{2}+i_{3}-2}$
| $\displaystyle+\sum_{1\leq j_{1}<j_{2}<j_{3}\leq
s}b_{j_{1}}b_{j_{2}}b_{j_{3}}c_{j_{1}}^{i_{1}-1}c_{j_{2}}^{i_{2}-1}c_{j_{2}}^{i_{3}-1}=\frac{1}{(i_{1}+i_{2}+i_{3})(i_{1}+i_{2})i_{1}}$
Table 2.3: Conditions (85) corresponding to multi-indices with up to three
indices.

For time-symmetric splitting methods, a set of independent order conditions
will be obtained by considering (85) restricted to Lyndon multi-indices
$(i_{1},\ldots,i_{k})$ with odd weight $i_{1}+\cdots+i_{k}$. For instance, the
subset of Lyndon multi-indices $(i_{1},\ldots,i_{k})$ such that
$1<i_{1}+\cdots+i_{k}\leq 5$ with odd weight $i_{1}+\cdots+i_{k}$ is

$\\{(3),(5),(1,2),(1,4),(2,3),(1,1,3),(1,2,2),(1,1,1,2)\\}.$

Notice that the treatment carried out in this subsection may also be formally
applied when $F_{1}$ is an unbounded operator. In that case, however, one has
to get rigorous estimates of the remainders to prove stability and convergence
of the corresponding schemes, as done in Thalhammer (2008), for example.

###  2.5 Order conditions IV: composition methods with Lyndon multi-indices

We now turn our attention to compositions (48) of a basic second-order time-
symmetric scheme $S_{h}^{[2]}$ with appropriate coefficients
$\gamma_{1},\ldots,\gamma_{s}$ chosen to achieve higher orders. Of course, a
set of conditions that guarantee that the scheme (48) attains a given order
can be obtained by rewriting that composition in terms of basic maps, such as
(66), and using the characterization of the order of the splitting method (66)
described in Subsection 2.4. The corresponding parameters $a_{j},b_{j}$ can be
obtained in terms of $\gamma_{1},\ldots,\gamma_{s}$ as follows:
$a_{1}=\gamma_{1}/2$, and

$a_{j+1}=\frac{\gamma_{j}+\gamma_{j+1}}{2},\qquad\quad
b_{j}=\gamma_{j},\qquad\mbox{ for }\,j=1,\ldots,s$

with $\gamma_{s+1}=0$. However, the resulting polynomial equations for any
given order $r$, once written in terms of the coefficients
$\gamma_{1},\ldots,\gamma_{s}$, are no longer independent. An alternative
formulation of the order of (48) in terms of explicit independent algebraic
equations in the coefficients $\gamma_{1},\ldots,\gamma_{s}$ will be presented
in Paragraph 2.5.2 below. This characterization is based on the treatment of
the more general composition (62), which is treated next.

####  2.5.1 Order conditions of compositions of a basic method and its adjoint

The composition (62) is at least of order $r$ if
$\Psi(h)-\exp(h(F_{1}+F_{2}))=\mathcal{O}(h^{r+1})$, where $\Psi(h)$ is the
associated Lie transformation (63).

To get the series expansion of $\Psi(h)$, we first consider the expansion in
powers of $h$ of the Lie transformation
$\mathcal{X}(h)=I+hX_{1}+h^{2}X_{2}+\cdots$ associated to the basic integrator
$\chi_{h}$. If $\chi_{h}$ is the Lie–Trotter scheme,
$\chi_{h}=\varphi_{h}^{[1]}\circ\varphi_{h}^{[2]}$, then

$\mathcal{X}(h)=\mathrm{e}^{hF_{2}}\mathrm{e}^{hF_{1}}=I+h(F_{1}+F_{2})+h^{2}(\frac{1}{2}F_{1}^{2}+F_{2}F_{1}+\frac{1}{2}F_{2}^{2})+\cdots.$

To deal with the most general problem, however, from now on we only assume
that $\chi_{h}$ is a smooth consistent integrator, so that
$X_{1}=F_{1}+F_{2}$, and each $X_{n}$ can be defined so that for each smooth
function $g$, $X_{n}g$ is a new smooth function given by (56).

Let us consider $\Psi_{0}(h)=I$, and for each $j\geq 1$,

$\begin{split}\Psi_{2j-1}(h)&=\mathcal{X}(-\alpha_{1}h)^{-1}\mathcal{X}(\alpha_{2}h)\cdots\mathcal{X}(-\alpha_{2j-1}h)^{-1}\\\
\Psi_{2j}(h)&=\mathcal{X}(-\alpha_{1}h)^{-1}\mathcal{X}(\alpha_{2}h)\cdots\mathcal{X}(-\alpha_{2j-1}h)^{-1}\mathcal{X}(\alpha_{2j}h),\end{split}$
(89)

so that in particular, $\Psi(h)=\Psi_{2s}(h)$. Notice that

$\displaystyle\mathcal{X}(-h)^{-1}$ $\displaystyle=I+\sum_{k\geq
1}(-1)^{k}(-hX_{1}+h^{2}X_{2}-h^{3}X_{3}+\cdots)^{k}$
$\displaystyle=I+hX_{1}+h^{2}(X_{1}^{2}-X_{2})+h^{3}(X_{1}^{3}-X_{1}X_{2}-X_{2}X_{1}+X_{3})+\cdots,$

which implies that for each $k\geq 1$, there exist polynomials
$w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{k})$ on the coefficients
$\alpha_{1},\ldots,\alpha_{k}$ such that

$\Psi_{k}(h)=I+\sum_{n\geq 1}h^{n}\sum_{m\geq
1}\sum_{i_{1}+\cdots+i_{m}=n}w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{k})\,X_{i_{1}}\cdots
X_{i_{m}}.$ (90)

We next determine such polynomial coefficients
$w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{k})$ recursively from the
relations $\Psi_{1}(h)=\mathcal{X}(-\alpha_{2j-1}h)^{-1}$ and

$\Psi_{2j-1}(h)=\Psi_{2j-2}(h)\mathcal{X}(-\alpha_{2j-1}h)^{-1},\quad\quad\Psi_{2j}(h)=\Psi_{2j-1}(h)\mathcal{X}(\alpha_{2j}h),$

or equivalently, $\Psi_{1}(h)\mathcal{X}(-\alpha_{1}h)=I$ and

$\Psi_{2j-1}(h)\mathcal{X}(-\alpha_{2j-1}h)=\Psi_{2j-2}(h),\quad\quad\Psi_{2j}(h)=\Psi_{2j-1}(h)\mathcal{X}(\alpha_{2j}h).$
(91)

Specifically, for arbitrary coefficients $w_{i_{1},\ldots,i_{m}}$ and
$\lambda$ one has

$\displaystyle\left(I+\sum_{n\geq 1}h^{n}\sum_{m\geq
1}\sum_{i_{1}+\cdots+i_{m}=n}w_{i_{1},\ldots,i_{m}}\,X_{i_{1}}\cdots
X_{i_{m}}\right)\mathcal{X}(\lambda h)=I+\sum_{i\geq
1}h^{i}(w_{i}+\lambda^{i})X_{n}$ $\displaystyle\qquad\qquad+\sum_{n\geq
1}h^{n}\sum_{m\geq
2}\sum_{i_{1}+\cdots+i_{m}=n}(w_{i_{1},\ldots,i_{m}}+\lambda^{i_{m}}w_{i_{1},\ldots,i_{m-1}})\,X_{i_{1}}\cdots
X_{i_{m}},$

so that, taking this expression into account, (91) leads to the following
identities

$\displaystyle w_{i}(\alpha_{1})$ $\displaystyle=-(-\alpha_{1})^{i},$
$\displaystyle w_{i}(\alpha_{1},\ldots,\alpha_{2j-1})$
$\displaystyle=w_{i}(\alpha_{1},\ldots,\alpha_{2j-2})-(-\alpha_{2j-1})^{i},$
$\displaystyle w_{i}(\alpha_{1},\ldots,\alpha_{2j})$
$\displaystyle=w_{i}(\alpha_{1},\ldots,\alpha_{2j-1})+\alpha_{2j}^{i},$
$\displaystyle w_{i_{1},\ldots,i_{m}}(\alpha_{1})$
$\displaystyle=-(-\alpha_{1})^{i_{m}}w_{i_{1},\ldots,i_{m-1}}(\alpha_{1}),$
$\displaystyle w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{2j-1})$
$\displaystyle=w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{2j-2})$
$\displaystyle-(-\alpha_{2j-1})^{i_{m}}w_{i_{1},\ldots,i_{m-1}}(\alpha_{1},\ldots,\alpha_{2j-1}),$
$\displaystyle w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{2j})$
$\displaystyle=w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{2j-1})$
$\displaystyle+\alpha_{2j}^{i_{m}}w_{i_{1},\ldots,i_{m-1}}(\alpha_{1},\ldots,\alpha_{2j-1}).$

Clearly, (90) holds for the coefficients
$w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{k})$ determined by the
relations above. Equivalently, the functions $w_{i_{1},\ldots,i_{m}}$ can be
defined as

$\begin{split}w_{i}(\alpha_{1},\ldots,\alpha_{2\ell})&=\sum_{j=1}^{s}\big{(}\alpha_{2j}^{i}-(-\alpha_{2j-1})^{i}\big{)},\\\
w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{2\ell})&=\sum_{j=1}^{s}\big{(}\alpha_{2j}^{i_{m}}-(-\alpha_{2j-1})^{i_{m}}\big{)}w_{i_{1},\ldots,i_{m-1}}(\alpha_{1},\ldots,\alpha_{2j-1}).\end{split}$
(92)

Notice that

$w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{k},0,\ldots,0)=w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{k}),$

as expected from (89), (90), and $\mathcal{X}(0)=I$.

Comparing the series expansion of $\Psi(h)=\Psi_{2s}(h)$ with

$\exp(h\,(F_{1}+F_{2}))=\exp(h\,X_{1})=I+\sum_{n\geq
1}\frac{h^{n}}{n!}X_{1}^{n},$

we finally conclude that the scheme (62) is of order at least $r$ if and only
if

$w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{2s})=\left\\{\begin{array}[]{ll}\frac{1}{m!}&\mbox{
if}\quad(i_{1},\ldots,i_{m})=\stackrel{{\scriptstyle
m}}{{\overbrace{(1,\ldots,1)}}}\\\ \\\ 0&\mbox{ otherwise},\end{array}\right.$
(93)

for each multi-index $(i_{1},\ldots,i_{k})$ such that $1\leq
i_{1}+\cdots+i_{k}\leq r$.

Furthermore, the order conditions (93) are not all independent. For instance,
it is straightforward to check from (92) that
$w_{i_{1},i_{2}}+w_{i_{2},i_{1}}+w_{i_{1}+i_{2}}=w_{i_{1}}w_{i_{2}}$ for
arbitrary indices $i_{1},i_{2}$. In particular, $2w_{1,1}+w_{2}=w_{1}^{2}$,
which implies that if the order conditions (93) for the multi-indices $(1)$
and $(2)$ are satisfied, then the condition for the multi-index $(1,1)$ is
automatically fulfilled. Actually, such dependences are similar to the shuffle
relations (79) that hold for the coefficients $v_{i_{1},\ldots,i_{m}}$
considered in Subsection 2.4. Indeed,

$w_{i_{1},\ldots,i_{n}}w_{i_{n+1},\ldots,i_{n+m}}=\sum_{\sigma\in\mathrm{Sh}(n,m)}w_{i_{\sigma(1)},\dots,i_{\sigma(n+m)}}+\cdots,$
(94)

where $\cdots$ refers to sums of products of coefficients corresponding to
multi-indices with $m-1$ or fewer indices. In fact, as shown in Chartier and
Murua (2009), the dependences (94) are directly related to the quasi-shuffle
product $*$ on the linear span of multi-indices introduced in Hoffman (2000),
and due to such dependences, it is enough to consider (93) for Lyndon multi-
indices $(i_{1},\ldots,i_{m})$ such that $i_{1}+\cdots+i_{k}\leq r$. That is,
the scheme (62) is of order at least $r$ if $\alpha_{1}+\cdots+\alpha_{2s}=1$
and for each Lyndon multi-index $(i_{1},\ldots,i_{k})$ such that
$1<i_{1}+\cdots+i_{k}\leq r$,

$w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{2s})=0.$ (95)

In particular, a method of order 3 must satisfy, besides consistency, the
conditions $w_{2}=w_{3}=w_{1,2}=0$.

This provides an alternative characterization of the order conditions of
general splitting schemes (64) in terms of the coefficients $\alpha_{j}$
obtained from the method parameters $a_{j},b_{j}$ from (65). Most importantly,
this also allows one to characterize the order of the scheme (48) obtained by
composing Strang maps. This will be presented next in Section 2.5.2.

####  2.5.2 Explicit characterization of the order conditions of scheme (48)

In Murua and Sanz-Serna (1999), such an explicit characterization is obtained
in terms of a set of polynomials indexed by certain sets of rooted trees
decorated by the set of odd positive integers. We now describe a related
formulation in terms of polynomials indexed by the set of Lyndon multi-indices
with odd indices based on the formalism developed in the previous subsection
for the more general composition (62).

We begin by considering (48) as the particular case of (62) with

$\alpha_{2j-1}=\alpha_{2j}=\gamma_{j}/2\quad\mbox{for}\quad j=1,\ldots,s.$
(96)

For each multi-index $(i_{1},\ldots,i_{k})$, we define the function
$u_{i_{1},\ldots,i_{m}}$ on the set of finite sequences
$(\gamma_{1},\ldots,\gamma_{s})$ of real numbers as follows:

$u_{i_{1},\ldots,i_{m}}(\gamma_{1},\ldots,\gamma_{s})=2^{i_{1}+\cdots+i_{m}-1}w_{i_{1},\ldots,i_{m}}(\alpha_{1},\ldots,\alpha_{2s})$
(97)

with $\alpha_{2j-1}=\alpha_{2j}=\gamma_{j}/2$. Clearly, scheme (48) is of
order at least $r$ if and only if $\gamma_{1}+\cdots+\gamma_{s}=1$ and

$u_{i_{1},\ldots,i_{m}}(\gamma_{1},\ldots,\gamma_{s})=0$ (98)

for each Lyndon multi-index $(i_{1},\ldots,i_{k})$ such that
$1<i_{1}+\cdots+i_{k}\leq r$. However, by definition,
$u_{i_{1},\ldots,i_{m}}(\gamma_{1},\ldots,\gamma_{s})\equiv 0$ if $i_{m}$ is
even. Moreover, for any multi-index $(i_{1},\ldots,i_{k})$ with some even
index, (98) holds provided that it holds for every Lyndon multi-index with
fewer indices.

Therefore, scheme (48) is of order at least $r$ if and only if
$\gamma_{1}+\cdots+\gamma_{s}=1$ and (98) for each Lyndon multi-index
$(i_{1},\ldots,i_{k})$ with odd indices such that $1<i_{1}+\cdots+i_{k}\leq
r$. For instance, the set of Lyndon multi-indices $(i_{1},\ldots,i_{k})$ of
odd indices such that $1<i_{1}+\cdots+i_{k}\leq 7$ is

$\\{(3),(5),(7),(1,3),(1,5),(1,1,3),(1,1,5),(1,3,3),(1,1,1,3),(1,1,1,1,3)\\}.$

The resulting number of order conditions, denoted as $m_{n}$, is gathered in
Table 2.1.

For multi-indices $(i_{1},\ldots,i_{m})$ with odd indices, the functions
$u_{i_{1},\ldots,i_{m}}$ can be written more explicitly as follows:

$\displaystyle u_{i}(\gamma_{1},\ldots,\gamma_{s})$
$\displaystyle=\sum_{j=1}^{s}\gamma_{j}^{i},$ (99) $\displaystyle
u_{i_{1},i_{2}}(\gamma_{1},\ldots,\gamma_{s})$
$\displaystyle=\sum_{j_{2}=1}^{s}\gamma_{j_{2}}^{i_{2}}\sum_{j_{1}=1}^{j_{2}*}\gamma_{j_{1}}^{i_{1}},$
(100) $\displaystyle u_{i_{1},i_{2},i_{3}}(\gamma_{1},\ldots,\gamma_{s})$
$\displaystyle=\sum_{j_{3}=1}^{s}\gamma_{j_{3}}^{i_{3}}\sum_{j_{2}=1}^{j_{3}*}\gamma_{j_{2}}^{i_{2}}\sum_{j_{1}=1}^{j_{2}*}\gamma_{j_{1}}^{i_{1}},$
(101) $\displaystyle
u_{i_{1},i_{2},i_{3},i_{4}}(\gamma_{1},\ldots,\gamma_{s})$
$\displaystyle=\sum_{j_{4}=1}^{s}\gamma_{j_{4}}^{i_{4}}\sum_{j_{3}=1}^{j_{4}*}\gamma_{j_{3}}^{i_{3}}\sum_{j_{2}=1}^{j_{3}*}\gamma_{j_{2}}^{i_{2}}\sum_{j_{1}=1}^{j_{2}*}\gamma_{j_{1}}^{i_{1}},$
(102)

and so on. Here, as in Murua and Sanz-Serna (1999), we have used the notation

$\sum_{j=1}^{k*}A_{j}=\frac{A_{k}}{2}+\sum_{j=1}^{k-1}A_{j}$

for $A_{1},\ldots,A_{k}\in\mathbb{R}$.

For time-symmetric integration schemes (48), a set of independent conditions
for even order $r$ will be obtained by considering (98) restricted to Lyndon
multi-indices $(i_{1},\ldots,i_{k})$ of odd indices and odd weight
$i_{1}+\cdots+i_{k}<r$. For instance, for order $8$, one only needs to
consider the following subset of Lyndon multi-indices:

$\\{(3),(5),(7),(1,1,3),(1,1,5),(1,3,3),(1,1,1,1,3)\\}.$

The above characterization of the order of the scheme (48) is also true in the
more general case where the Strang map $S^{[2]}_{h}$ is replaced by an
arbitrary time-symmetric integrator of order $2\ell$. In that case, only
Lyndon multi-indices with indices from the set