# High order integrators obtained by linear combinations of symmetric-
conjugate compositions

F. Casas1, A. Escorihuela-Tomàs2
1 Departament de Matemàtiques and IMAC, Universitat Jaume I, 12071-Castellón,
Spain
Email<EMAIL_ADDRESS>
2 Departament de Matemàtiques, Universitat Jaume I, 12071-Castellón, Spain
Email<EMAIL_ADDRESS>

###### Abstract

A new family of methods involving complex coefficients for the numerical
integration of differential equations is presented and analyzed. They are
constructed as linear combinations of symmetric-conjugate compositions
obtained from a basic time-symmetric integrator of order $2n$ ($n\geq 1$). The
new integrators are of order $2(n+k)$, $k=1,2,\ldots$, and preserve time-
symmetry up to order $4n+3$ when applied to differential equations with real
vector fields. If in addition the system is Hamiltonian and the basic scheme
is symplectic, then they also preserve symplecticity up to order $4n+3$. We
show that these integrators are well suited for a parallel implementation,
thus improving their efficiency. Methods up to order 10 based on a 4th-order
integrator are built and tested in comparison with other standard procedures
to increase the order of a basic scheme.

AMS numbers: 65L05, 65P10, 37M15

Keywords: Composition methods, symmetric-conjugate compositions, complex
coefficients. preservation of properties, parabolic equations

##  1 Introduction

Composition methods constitute a standard tool to construct high-order
numerical integrators for the initial value problem

$\dot{x}=f(x),\qquad x(t_{0})=x_{0}\in\mathbb{R}^{d},$ (1.1)

in particular when the vector field $f$ possesses some qualitative property
whose preservation by numerical approximations is deemed relevant [6, 16]. Let
$\mathcal{S}_{h}^{[2n]}$ denote a $2n$-th order method, so that
$\mathcal{S}_{h}^{[2n]}(x_{0})=\varphi_{h}(x_{0})+\mathcal{O}(h^{2n+1})$,
where $x(h)=\varphi_{h}(x_{0})$ is the exact solution of Eq. (1.1) for a time
step $h$. Then, if the coefficients $\alpha_{1},\alpha_{2},\ldots,\alpha_{s}$
satisfy some algebraic conditions, the composition of the basic scheme with
step sizes $\alpha_{1}h,\alpha_{2}h,\ldots,\alpha_{s}h$, i.e.,

$\psi_{h}=\mathcal{S}_{\alpha_{1}h}^{[2n]}\circ\mathcal{S}_{\alpha_{2}h}^{[2n]}\circ\cdots\circ\mathcal{S}_{\alpha_{s-1}h}^{[2n]}\circ\mathcal{S}_{\alpha_{s}h}^{[2n]}$
(1.2)

is a new method of higher order $2n+m$ [9]. If in particular $f$ is
Hamiltonian and $\mathcal{S}_{h}^{[2n]}$ is symplectic, then the composition
method (1.2) is also symplectic [16]. In general, any geometric property the
basic method has in common with the exact solution is still shared by the
higher-order scheme (1.2) if this property is preserved by composition [18].
Moreover, suppose $\mathcal{S}_{h}^{[2n]}$ is time-symmetric, namely, it
satisfies

$\mathcal{S}_{h}^{[2n]}\circ\mathcal{S}_{-h}^{[2n]}=\mathrm{id},$

where $\mathrm{id}$ is the identity map, for any $h$. Then, method (1.2) is
also time-symmetric if the composition is left-right palindromic, i.e.,
$\alpha_{s+1-j}=\alpha_{j}$, $j=1,2,\ldots$.

A well known class of composition methods is obtained by applying the triple-
jump procedure [21, 24]:

$\mathcal{S}_{h}^{[2n+2]}=\mathcal{S}_{\alpha_{1}h}^{[2n]}\circ\mathcal{S}_{\alpha_{2}h}^{[2n]}\circ\mathcal{S}_{\alpha_{1}h}^{[2n]},$
(1.3)

with

$\alpha_{1}=\frac{1}{2-2^{1/(2n+1)}},\qquad\alpha_{2}=1-2\alpha_{1},$ (1.4)

is a new method of order $2n+2$. The same technique can be applied again to
$\mathcal{S}_{h}^{[2n+2]}$, so that one can construct recursively time-
symmetric methods of any order $2n+2k$, $k=1,2,\ldots$.

When constructing high-order composition methods, real coefficients
$\alpha_{1},\ldots,\alpha_{s}$ are not the only option, however. In fact, the
unavoidable existence of negative $\alpha_{j}$ in (1.2) when the order is
higher than two [5, 15, 20, 22] typically imposes stability restrictions on
the step size. This occurs in particular when Eq. (1.1) is the outcome of a
parabolic differential equation discretized in space. In that case,
considering complex coefficients with positive real part is also a valid
alternative [12, 17]. Even for problems where the presence of some
$\alpha_{j}<0$ is not particularly troublesome, composition methods with
complex coefficients have also been proposed and analyzed from the
preservation of properties viewpoint [10, 7, 13].

In the particular case of the triple-jump composition (1.3), in addition to
the real solution (1.4), the complex one with the smallest phase is

$\alpha_{1}=\frac{\mathrm{e}^{i\pi/(2n+1)}}{2^{1/(2n+1)}-2\mathrm{e}^{i\pi/(2n+1)}},\qquad\alpha_{2}=1-2\alpha_{1},$
(1.5)

and the resulting method has in fact smaller truncation errors than its real
counterpart (1.4). If the basic scheme is time-symmetric and of order 2, then
time-symmetric methods up to order 14 with coefficients having positive real
part are possible by applying this technique [8].

The order can be raised by one instead with the simplest composition [3, 22]

$\psi_{h}^{[2n+1]}=\mathcal{S}_{\alpha_{1}h}^{[2n]}\circ\mathcal{S}_{\alpha_{2}h}^{[2n]}$
(1.6)

if

$\alpha_{1}=\bar{\alpha}_{2}=\frac{1}{2}+\frac{i}{2}\,\frac{\sin\frac{2\ell+1}{2n+1}\pi}{1+\cos\frac{2\ell+1}{2n+1}\pi}\qquad\mbox{
for }\qquad-n\leq\ell\leq n-1.$

The choice $\ell=0$ gives the solution with the smallest phase, which we
denote by $\gamma^{[2n]}$:

$\alpha_{1}=\gamma^{[2n]}:=\frac{1}{2}+\frac{i}{2}\,\frac{\sin\frac{\pi}{2n+1}}{1+\cos\frac{\pi}{2n+1}},\qquad
n=1,2,\ldots$ (1.7)

When the vector field $f$ in (1.1) is real, then
$x_{1}=\psi_{h}^{[2n+1]}(x_{0})$ is complex, and so it is quite natural to
project $x_{1}$ on the real axis and proceed to the next step only with
$\Re(x_{1})$. This is equivalent of course to integrating with the scheme

$R_{h}^{(1)}=\frac{1}{2}\left(\psi_{h}^{[2n+1]}+\overline{\psi}_{h}^{[2n+1]}\right).$
(1.8)

Method (1.8) is not time-symmetric, even when $\mathcal{S}_{h}^{[2n]}$ is.
Nevertheless, it has been shown in [11] that $R_{h}^{(1)}$ is _pseudo-
symmetric_ of order $4n+3$, in the sense that

$R_{h}^{(1)}\circ R_{-h}^{(1)}=\mathrm{id}+\mathcal{O}(h^{4n+4})$

if the vector field $f$ in (1.1) is real. If in addition $f$ is Hamiltonian
and $\mathcal{S}_{h}^{[2n]}$ is symplectic, then $R_{h}^{(1)}$ is also
_pseudo-symplectic_ of order $4n+3$. In other words, projecting
$\psi_{h}^{[2n+1]}$ at each integration step leads to a numerical method that
preserves geometric properties of the exact solution up to an order that is
much higher than the order of the method itself. Pseudo-symplectic integrators
have been previously considered in the literature, both in the context of
Runge–Kutta [2] and polynomial extrapolation methods [4, 14].

Moreover, as shown in [11], $R_{h}^{(1)}$ can be taken as the basis of the
recursion

$R_{h}^{(k)}=\frac{1}{2}\left(R_{\gamma^{[2k]}h}^{(k-1)}\circ
R_{\bar{\gamma}^{[2k]}h}^{(k-1)}+R_{\bar{\gamma}^{[2k]}h}^{(k-1)}\circ
R_{\gamma^{[2k]}h}^{(k-1)}\right),\qquad\qquad k=2,3\ldots,$ (1.9)

producing methods of order $2(n+k)$, also pseudo-symmetric of order $4n+3$.
Here the coefficients $\gamma^{[2k]}$ are given by Eq. (1.7). For future
reference, we call (1.9) $R$-methods.

Scheme (1.6) is a particular example of a _symmetric-conjugate_ composition.
These are composition methods of the form

$\psi_{h}=\mathcal{S}_{\alpha_{1}h}^{[2n]}\circ\mathcal{S}_{\alpha_{2}h}^{[2n]}\circ\cdots\circ\mathcal{S}_{\bar{\alpha}_{2}h}^{[2n]}\circ\mathcal{S}_{\bar{\alpha}_{1}h}^{[2n]},$
(1.10)

i.e., compositions (1.2) with $\alpha_{j}\in\mathbb{C}$ and

$\bar{\alpha}_{s+1-j}=\alpha_{j},\qquad j=1,2,\ldots.$

Methods of this class, as shown in [7], possess remarkable preservation
properties when considering its real part,

$\Re(\psi_{h})=\frac{1}{2}\left(\psi_{h}+\overline{\psi}_{h}\right).$

In particular, if one takes a time-symmetric 2nd-order scheme as the basic
method and the coefficients $\alpha_{1},\alpha_{2},\ldots$ are chosen in such
a way that $\psi_{h}$ is of order $2n-1$, then $\Re(\psi_{h})$ is of order
$2n$ and pseudo-symmetric of order $4n-1$ when the vector field $f$ in (1.1)
is real. If in addition $f$ is a (real) Hamiltonian vector field and
${\mathcal{S}}_{h}^{[2]}$ is a symplectic integrator, then $\Re(\psi_{h})$ is
pseudo-symplectic of order $4n-1$.

Since taking the real part of a symmetric-conjugate method is just a very
special linear combination, it is quite natural to ask what happens when one
considers a more general linear combination of symmetric-conjugate
compositions and their complex-conjugate, $\psi_{h}^{(j)}$,
$\overline{\psi}_{h}^{(j)}$: is it possible to construct new methods of higher
order whereas still preserving time-symmetry (and symplecticity) up to the
order prescribed by the composition $\psi_{h}^{(j)}$? If yes, how the new
methods are built? Addressing these questions is precisely the subject of the
present paper. In doing so, we present a new family of schemes of increasingly
higher order well adapted for implementation in a parallel environment,
requiring less computational effort than the $R$-methods (1.9) but with the
same qualitative properties.

If we denote for simplicity the symmetric-conjugate composition (1.10) by its
sequence of coefficients,

$\psi_{h}^{(j)}=(\alpha_{1},\alpha_{2},\ldots,\alpha_{s-1},\alpha_{s}),$

with $\bar{\alpha}_{s+1-j}=\alpha_{j}$, these new schemes have the basic
structure

$T_{h}^{(k)}=\frac{1}{2^{k}}\sum_{j=1}^{2^{k-1}}\Big{(}(\alpha_{j_{2^{k}}},\ldots,\alpha_{j_{1}})+\mbox{c.c.}\Big{)}$
(1.11)

and are of order $2(n+k)\leq 4n+3$ and pseudo-symmetric of order $4n+3$. We
designate them as $T$-methods.

##  2 Construction of the family of $T$-methods

In this section we construct the new family of integrators $T_{h}^{(k)}$ and
show explicitly that they are of order $2(n+k)$ and pseudo-symmetric of order
$4n+3$ for $k=1,2,3$. The same procedure can be formally extended to any
$k>3$. The analysis is based on the Lie formalism applied to the series of
differential operators associated to the integrators.

###  2.1 Series of differential operators

As is well known, given a time-symmetric integrator $\mathcal{S}_{h}^{[2n]}$
of order $2n\geq 2$ one can associate a series of linear operators
$\exp(Y(h))$ so that

$g(\mathcal{S}_{h}^{[2n]}(x))=\exp(Y(h))[g](x)$

for all functions $g$ [9], with

$Y(h)=hY_{1}+h^{2n+1}Y_{2n+1}+h^{2n+3}Y_{2n+3}+\cdots.$

Here $Y_{k}$ are certain operators depending on the particular method and, for
consistency, $Y_{1}=F$, where $F$ is the Lie derivative associated with $f$:

$F=\sum_{i\geq 1}\,f_{i}(x)\,\frac{\partial}{\partial x_{i}}.$ (2.1)

The composition (1.2) then has the associated series

$\Psi(h)=\exp(Y(h\alpha_{s}))\,\exp(Y(h\alpha_{s-1}))\,\cdots\,\exp(Y(h\alpha_{2}))\,\exp(Y(h\alpha_{1})),$
(2.2)

which can be formally written as $\Psi(h)=\exp(V(h))$ by repeated application
of the Baker–Campbell–Hausdorff formula, with

$V(h)=hF+h^{2n+1}V_{2n+1}+h^{2n+2}V_{2n+2}+\cdots.$

Here $V_{2n+1},V_{2n+2},\ldots$ are linear combinations of Lie brackets
involving the operators $Y_{1},Y_{2n+1},Y_{2n+3},\ldots$ [18]. In the
particular case of a symmetric-conjugate composition (1.10), terms $V_{2k}$ in
$V(h)$ of even powers in $h$ are pure imaginary, whereas terms $V_{2k+1}$ are
real [7].

For a consistent symmetric-conjugate composition (1.10), i.e., verifying

$\sum_{j=1}^{s}\alpha_{j}=\alpha_{1}+\alpha_{2}+\cdots+\bar{\alpha}_{2}+\bar{\alpha}_{1}=1,$
(2.3)

we get explicitly

$V(h)=hE_{1,1}+h^{2n}\sum_{j\geq
0}h^{2j+1}\sum_{k=1}^{\ell_{2j+1}}\mu_{2n+2j+1,k}E_{2j+1,k}+i\,h^{2n}\sum_{j\geq
1}h^{2j}\sum_{k=1}^{\ell_{2j}}\sigma_{2n+2j,k}E_{2j,k},$ (2.4)

where $\mu_{n,k}$, $\sigma_{n,k}$ are homogeneous real polynomials of degree
$n$ in the coefficients $\alpha_{l}$, $l=1,\ldots,s$, and $E_{n,k}$ are
elements $Y_{j}$ and independent Lie brackets involving these operators. In
particular

$\mu_{2n+2j+1,1}=\sum_{l=1}^{s}\alpha_{l}^{2(n+j)+1},\quad\quad j\geq 0$

and

$E_{1,1}=Y_{1},\quad E_{2n+2n+2j+1,1}=Y_{2n+2j+1},\quad
E_{2n+2j,1}=[E_{1,1},E_{2n+(2j-1),1}],\quad j=1,2,\ldots$

###  2.2 Linear combinations of symmetric-conjugate compositions

Let us now consider the linear combination

$\phi_{h}=\frac{1}{2k}\sum_{j=1}^{k}\left(\psi_{h}^{(j)}+\overline{\psi}_{h}^{(j)}\right),$
(2.5)

where each $\psi_{h}^{(j)}$ is a consistent symmetric-conjugate composition of
the form (1.10) with different coefficients $\alpha_{k}^{(j)}$. Then, clearly,
$\phi_{h}$ has

$\Phi(h)\equiv\frac{1}{2k}\sum_{j=1}^{k}\left(\Psi^{(j)}(h)+\overline{\Psi}^{(j)}(h)\right)=\frac{1}{2k}\sum_{j=1}^{k}\left(\mathrm{e}^{V_{j}(h)}+\mathrm{e}^{\overline{V}_{j}(h)}\right)$
(2.6)

as the associated series of operators, where each $V_{j}(h)$ is of the form
(2.4). Now, by following the same approach as in [11], we express $\Phi(h)$ as

$\Phi(h)=\frac{1}{2k}\mathrm{e}^{\frac{h}{2}F}\,\sum_{j=1}^{k}\left(\mathrm{e}^{W_{j}(h)}+\mathrm{e}^{\overline{W}_{j}(h)}\right)\,\mathrm{e}^{\frac{h}{2}F},$

where

$\displaystyle
W_{j}(h)=h^{2n+1}\mu_{2n+1,1}^{(j)}E_{2n+1,1}+i\,h^{2n+2}\sigma_{2n+2,1}^{(j)}E_{2n+2,1}$
(2.7)
$\displaystyle\qquad+h^{2n+3}\left(\mu_{2n+3,1}^{(j)}E_{2n+3,1}+\Big{(}\mu_{2n+3,2}^{(j)}+\frac{1}{24}\mu_{2n+1,1}^{(j)}\Big{)}\,E_{2n+3,2}\right)$
$\displaystyle\qquad+i\,h^{2n+4}\left(\sigma_{2n+4,1}^{(j)}E_{2n+4,1}+\Big{(}\sigma_{2n+4,2}^{(j)}+\frac{1}{24}\sigma_{2n+2,1}^{(j)}\Big{)}\,E_{2n+4,2}\right)$
$\displaystyle\qquad+h^{2n+5}\left(\mu_{2n+5,1}^{(j)}E_{2n+5,1}+\Big{(}\mu_{2n+5,2}^{(j)}+\frac{1}{24}\mu_{2n+3,1}^{(j)}\Big{)}\,E_{2n+5,2}\,+\right.$
$\displaystyle\qquad\qquad\quad\left.\Big{(}\mu_{2n+5,3}^{(j)}+\frac{1}{24}\mu_{2n+3,2}^{(j)}+\frac{1}{1920}\mu_{2n+1,1}^{(j)}\Big{)}\,E_{2n+5,3}\right)$
$\displaystyle\qquad+i\,h^{2n+6}\left(\sigma_{2n+6,1}^{(j)}E_{2n+6,1}+\Big{(}\sigma_{2n+6,2}^{(j)}+\frac{1}{24}\sigma_{2n+4,1}^{(j)}\Big{)}\,E_{2n+6,2}+\right.$
$\displaystyle\qquad\qquad\quad\left.\Big{(}\sigma_{2n+6,3}^{(j)}+\frac{1}{24}\sigma_{2n+4,2}^{(j)}+\frac{1}{1920}\sigma_{2n+2,1}^{(j)}\Big{)}\,E_{2n+6,3}\right)$
$\displaystyle\qquad+\mathcal{O}(h^{2n+7}).$

Here

$\displaystyle E_{2n+3,2}=[E_{1,1},E_{2n+2,1}],\quad
E_{2n+4,2}=[E_{1,1},E_{2n+3,2}],\quad E_{2n+5,2}=[E_{1,1},E_{2n+4,1}],$
$\displaystyle E_{2n+6,2}=[E_{1,1},E_{2n+5,1}],\quad
E_{2n+5,3}=[E_{1,1},E_{2n+4,2}],\quad E_{2n+6,2}=[E_{1,1},E_{2n+5,2}].$

This is done by applying the symmetric Baker–Campbell–Hausdorff formula to
each product
$\mathrm{e}^{-\frac{h}{2}F}\,\mathrm{e}^{V_{j}(h)}\,\mathrm{e}^{-\frac{h}{2}F}$.
From (2.7), a straightforward calculation shows that

$\displaystyle(W_{j}+\overline{W}_{j})^{2}=4h^{4n+2}(\mu_{2n+1,1}^{(j)})^{2}\,E_{2n+1,1}^{2}+\mathcal{O}(h^{4n+4})$
$\displaystyle
W_{j}^{2}+\overline{W}_{j}^{2}=2h^{4n+2}(\mu_{2n+1,1}^{(j)})^{2}\,E_{2n+1,1}^{2}+\mathcal{O}(h^{4n+4}).$

Therefore,

$\frac{1}{2}\left(\mathrm{e}^{W_{j}}+\mathrm{e}^{\overline{W}_{j}}\right)-\mathrm{e}^{\frac{1}{2}(W_{j}+\overline{W}_{j})}=\frac{1}{4}(W_{j}^{2}+\overline{W}_{j}^{2})-\frac{1}{8}(W_{j}+\overline{W}_{j})^{2}+\cdots=\mathcal{O}(h^{4n+4})$

and $\Phi(h)$ can also be written as

$\Phi(h)=\frac{1}{k}\,\sum_{j=1}^{k}\mathrm{e}^{\frac{h}{2}F}\,\mathrm{e}^{\frac{1}{2}(W_{j}(h)+\overline{W}_{j}(h))}\,\mathrm{e}^{\frac{h}{2}F}+\mathcal{O}(h^{4n+4}).$

In consequence, each term in $\phi_{h}$ is time-symmetric up to terms
$h^{4n+3}$, with independence of the polynomials $\mu_{k,l}^{(j)}$,
$\sigma_{k,l}^{(j)}$, since the sum $W_{j}(h)+\overline{W}_{j}(h)$ only
contains odd powers of $h$.

On the other hand, one has

$\displaystyle\frac{1}{2k}\sum_{j=1}^{k}\left(\mathrm{e}^{W_{j}}+\mathrm{e}^{\overline{W}_{j}}\right)-\exp\left(\frac{1}{2k}\sum_{j=1}^{k}(W_{j}+\overline{W}_{j})\right)$
$\displaystyle=\frac{1}{4k}\sum_{j=1}^{k}(W_{j}^{2}+\overline{W}_{j}^{2})-\frac{1}{8k^{2}}\left(\sum_{j=1}^{k}(W_{j}+\overline{W}_{j})\right)^{2}+\cdots$
$\displaystyle=h^{4n+2}\,\frac{1}{2k}\left(\sum_{j=1}^{k}(\mu_{2n+1}^{(j)})^{2}-\frac{1}{k}\Big{(}\sum_{j=1}^{k}\mu_{2n+1,1}^{(j)}\Big{)}^{2}\right)\,E_{2n+1,1}^{2}+\mathcal{O}(h^{4n+4}),$

so that it is also true that

$\Phi(h)=\exp\left(\frac{h}{2}F\right)\exp\left(\frac{1}{2k}\sum_{j=1}^{k}\big{(}W_{j}(h)+\overline{W}_{j}(h)\big{)}\right)\exp\left(\frac{h}{2}F\right)+\mathcal{O}(h^{4n+2}).$
(2.8)

###  2.3 Order conditions

It is thus possible to obtain the order conditions for the method $\phi_{h}$
in (2.5) by analyzing just the exponent of the central term in (2.8). From
(2.7) it follows that

$\displaystyle\frac{1}{2}\sum_{j=1}^{k}\left(W_{j}(h)+\overline{W}_{j}(h)\right)=h^{2n+1}\,c_{2n+1,1}\,E_{2n+1,1}+h^{2n+3}c_{2n+3,1}E_{2n+3,1}$
$\displaystyle\qquad+h^{2n+3}\left(c_{2n+3,2}+\frac{1}{24}c_{2n+1,1}\right)\,E_{2n+3,2}+h^{2n+5}\,c_{2n+5,1}\,E_{2n+5,1}$
$\displaystyle\qquad+h^{2n+5}\left(c_{2n+5,2}+\frac{1}{24}c_{2n+3,1}\right)\,E_{2n+5,2}$
$\displaystyle\qquad+h^{2n+5}\left(c_{2n+5,3}+\frac{1}{24}c_{2n+3,2}+\frac{1}{1920}c_{2n+1,1}\right)\,E_{2n+5,3}+\mathcal{O}(h^{2n+7}),$

with

$\displaystyle c_{2n+1,1}=\sum_{j=1}^{k}\mu_{2n+1,1}^{(j)},\qquad
c_{2n+3,1}=\sum_{j=1}^{k}\mu_{2n+3,1}^{(j)},\qquad
c_{2n+3,2}=\sum_{j=1}^{k}\mu_{2n+3,2}^{(j)}$ $\displaystyle
c_{2n+5,1}=\sum_{j=1}^{k}\mu_{2n+5,1}^{(j)},\qquad
c_{2n+5,2}=\sum_{j=1}^{k}\mu_{2n+5,2}^{(j)},\qquad
c_{2n+5,3}=\sum_{j=1}^{k}\mu_{2n+5,3}^{(j)}$

In consequence, for consistent compositions $\psi_{h}^{(j)}$, $j=1,\ldots,k$,
the conditions to be satisfied so that $\phi_{h}$ is a method of order $r$ are
the following:

  * •

$r=2n+2$: $\,c_{2n+1,1}=0$

  * •

$r=2n+4$: $\,c_{2n+1,1}=c_{2n+3,1}=c_{2n+3,2}=0$

  * •

$r=2n+6$:
$\,c_{2n+1,1}=c_{2n+3,1}=c_{2n+3,2}=c_{2n+5,1}=c_{2n+5,2}=c_{2n+5,3}=0$

###  2.4 New schemes

Once identified the relevant order conditions, our next goal is to solve these
equations with the minimum number of basic schemes in the compositions
$\psi_{h}^{(j)}$ and the minimum value of $k$ in the linear combination (2.5).

#### Order $r=2n+2$.

One needs to solve two equations to get a method $\phi_{h}$ of order $2n+2$:
consistency and $c_{2n+1,1}=0$. These can be satisfied by taking $k=1$ and the
simplest composition
$\psi_{h}=\mathcal{S}_{\alpha_{1}h}^{[2n]}\circ\mathcal{S}_{\bar{\alpha}_{1}h}^{[2n]}$,
in which case one has

$\alpha_{1}+\bar{\alpha}_{1}=1,\qquad\alpha_{1}^{2n+1}+\bar{\alpha}_{1}^{2n+1}=0.$

In other words, we recover the composition (1.6) and the $R$-method (1.8). Our
first $T$-method (1.11) is thus

$T_{h}^{(1)}=\frac{1}{2}\Big{(}(\gamma^{[2n]},\bar{\gamma}^{[2n]})+(\bar{\gamma}^{[2n]},\gamma^{[2n]})\Big{)}$
(2.9)

or in more detail

$T_{h}^{(1)}=\frac{1}{2}\Big{(}{\mathcal{S}}_{\gamma^{[2n]}h}^{[2n]}\circ{\mathcal{S}}_{\bar{\gamma}^{[2n]}h}^{[2n]}+{\mathcal{S}}_{\bar{\gamma}^{[2n]}h}^{[2n]}\circ{\mathcal{S}}_{\gamma^{[2n]}h}^{[2n]}\Big{)}.$

#### Order $r=2n+4$.

Now we have to solve 3 order conditions in addition to consistency for the
compositions $\psi_{h}^{(j)}$ involved. As before, one could take in principle
$k=1$. In that case, the minimum number of basic maps in $\psi_{h}^{(1)}$ is
4, just to have enough parameters to satisfy the order conditions. It turns
out, however, that there are no solutions with the required symmetry
$\alpha_{4}=\bar{\alpha}_{1}$, $\alpha_{3}=\bar{\alpha}_{2}$. In fact, if we
take

$\psi_{h}^{(1)}=(\bar{\alpha}_{1},\bar{\alpha}_{2},\alpha_{2},\alpha_{1}),\qquad\mbox{
with
}\qquad\alpha_{1}=\bar{\gamma}^{[2n+4]}\bar{\gamma}^{[2n+2]},\quad\alpha_{2}=\bar{\gamma}^{[2n+4]}\gamma^{[2n]},$

then $\mu_{2n+1,1}^{(1)}=\mu_{2n+3,1}^{(1)}=0$, but $\mu_{2n+3,2}^{(1)}\neq
0$. On the other hand, if we take

$\psi_{h}^{(2)}=(\bar{\alpha}_{2},\bar{\alpha}_{1},\alpha_{1},\alpha_{2})$

with the _same_ values of $\alpha_{1}$, $\alpha_{2}$ as before, then
$\mu_{2n+3,2}^{(2)}=-\mu_{2n+3,2}^{(1)}$, whereas still verifying that
$\mu_{2n+1,1}^{(2)}=\mu_{2n+3,1}^{(2)}=0$. In consequence, by combining both
compositions,

$\phi_{h}=\frac{1}{4}\left(\psi_{h}^{(1)}+\overline{\psi}_{h}^{(1)}+\psi_{h}^{(2)}+\overline{\psi}_{h}^{(2)}\right),$

one gets a method of order $2n+4$ and pseudo-symmetric of order $4n+3$. This
corresponds to our second $T$-method, which reads explicitly

$\displaystyle T_{h}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{1}{4}\Big{(}(\gamma^{[2n+2]}\gamma^{[2n]},\gamma^{[2n+2]}\bar{\gamma}^{[2n]},\bar{\gamma}^{[2n+2]}\gamma^{[2n]},\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]})$
$\displaystyle\quad+(\gamma^{[2n+2]}\bar{\gamma}^{[2n]},\gamma^{[2n+2]}\gamma^{[2n]},\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]},\bar{\gamma}^{[2n+2]}\gamma^{[2n]})$
$\displaystyle\quad+(\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]},\bar{\gamma}^{[2n+2]}\gamma^{[2n]},\gamma^{[2n+2]}\bar{\gamma}^{[2n]},\gamma^{[2n+2]}\gamma^{[2n]})$
$\displaystyle\quad+(\bar{\gamma}^{[2n+2]}\gamma^{[2n]},\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]},\gamma^{[2n+2]}\gamma^{[2n]},\gamma^{[2n+2]}\bar{\gamma}^{[2n]})\Big{)}.$

Again, the coefficients $\gamma^{[2m]}$ are given by Eq. (1.7).

#### Order $r=2n+6$.

A total of 7 equations (including consistency) have to be solved in this case,
so that we take a symmetric-conjugate composition involving $s=8$ basic maps,

$\psi_{h}^{(1)}=(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\bar{\alpha}_{4},\bar{\alpha}_{3},\bar{\alpha}_{2},\bar{\alpha}_{1}).$

With the choice

$\displaystyle\alpha_{1}=\gamma^{[2n+4]}\gamma^{[2n+2]}\gamma^{[2n]},\qquad\quad\alpha_{2}=\gamma^{[2n+4]}\gamma^{[2n+2]}\bar{\gamma}^{[2n]},$
$\displaystyle\alpha_{3}=\gamma^{[2n+4]}\bar{\gamma}^{[2n+2]}\gamma^{[2n]},\qquad\quad\alpha_{4}=\gamma^{[2n+4]}\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]}$

it turns out that conditions $c_{2n+1}=c_{2n+3,1}=c_{2n+5,1}=0$ are
automatically satisfied. By following the same approach as before, we permute
the position of the coefficients and take the composition

$\psi_{h}^{(2)}=(\alpha_{2},\alpha_{1},\alpha_{4},\alpha_{3},\bar{\alpha}_{3},\bar{\alpha}_{4},\bar{\alpha}_{1},\bar{\alpha}_{2}).$

Then, one has $\mu_{2n+3,2}^{(2)}=-\mu_{2n+3,2}^{(1)}$, so that
$\psi_{h}^{(1)}+\psi_{h}^{(2)}$ leads to a method of order $2n+4$. More
composition have to be incorporated, however, in order to verify conditions
$c_{2n+5,2}=0$ and $c_{2n+5,3}=0$. The former is accomplished by both sums
$\psi_{h}^{(1)}+\psi_{h}^{(4)}$ and $\psi_{h}^{(2)}+\psi_{h}^{(3)}$, where

$\displaystyle\psi_{h}^{(3)}=(\alpha_{3},\alpha_{4},\alpha_{1},\alpha_{2},\bar{\alpha}_{2},\bar{\alpha}_{1},\bar{\alpha}_{4},\bar{\alpha}_{3})$
$\displaystyle\psi_{h}^{(4)}=(\alpha_{4},\alpha_{3},\alpha_{2},\alpha_{1},\bar{\alpha}_{1},\bar{\alpha}_{2},\bar{\alpha}_{3},\bar{\alpha}_{4}),$

but the later is satisfied only by adding up the four compositions. In
summary, the linear combination

$\frac{1}{4}\Re(\psi_{h}^{(1)}+\psi_{h}^{(2)}+\psi_{h}^{(3)}+\psi_{h}^{(4)})$

leads to a method of order $2n+6$, denoted as $T_{h}^{(3)}$. More explicitly,

$T_{h}^{(3)}=\frac{1}{8}\Big{(}\psi_{h}^{(1)}+\psi_{h}^{(2)}+\psi_{h}^{(3)}+\psi_{h}^{(4)}+\overline{\psi}_{h}^{(1)}+\overline{\psi}_{h}^{(2)}+\overline{\psi}_{h}^{(3)}+\overline{\psi}_{h}^{(4)}\Big{)}.$
(2.11)

The same procedure can be carried out in general, although more order
conditions (and consequently more compositions involving more basic maps) have
to be dealt with. This class of methods can be represented in a convenient way
as follows. If we introduce the matrix of coefficients

$\Gamma_{2n}:=\frac{1}{2}\left(\begin{array}[]{cc}\gamma^{[2n]}&\bar{\gamma}^{[2n]}\\\
\bar{\gamma}^{[2n]}&\gamma^{[2n]}\end{array}\right)$

then, according with the previous results, method $T_{h}^{(1)}$ (of order
$2n+2$) can be represented by $\Gamma_{2n}$,

$T_{h}^{(1)}\leadsto\Gamma_{2n},$

whereas $T_{h}^{(2)}$ (of order $2n+4$) can be associated with the matrix

$\Gamma_{2n+2}\otimes\Gamma_{2n}=\frac{1}{4}\left(\begin{array}[]{cccc}\gamma^{[2n+2]}\gamma^{[2n]}&\gamma^{[2n+2]}\bar{\gamma}^{[2n]}&\bar{\gamma}^{[2n+2]}\gamma^{[2n]}&\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]}\\\
\gamma^{[2n+2]}\bar{\gamma}^{[2n]}&\gamma^{[2n+2]}\gamma^{[2n]}&\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]}&\bar{\gamma}^{[2n+2]}\gamma^{[2n]}\\\
\bar{\gamma}^{[2n+2]}\gamma^{[2n]}&\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]}&\gamma^{[2n+2]}\gamma^{[2n]}&\gamma^{[2n+2]}\bar{\gamma}^{[2n]}\\\
\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]}&\bar{\gamma}^{[2n+2]}\gamma^{[2n]}&\gamma^{[2n+2]}\bar{\gamma}^{[2n]}&\gamma^{[2n+2]}\gamma^{[2n]}\end{array}\right),$

in the sense that each file of $\Gamma_{2n+2}$ corresponds to a particular
symmetric-conjugate composition entering into the formulation of
$T_{h}^{(2)}$. We can write analogously

$T_{h}^{(2)}\leadsto\Gamma_{2n+2}\otimes\Gamma_{2n},$

and moreover

$T_{h}^{(3)}\leadsto\Gamma_{2n+4}\otimes(\Gamma_{2n+2}\otimes\Gamma_{2n}).$

In general, the coefficients in the $T$-method of order $r=2n+2k$ are
distributed according with the pattern

$T_{h}^{(k)}\leadsto\Gamma_{2(n+k-1)}\otimes(\Gamma_{2(n+k-2)}\otimes\cdots\otimes(\Gamma_{2n+2}\otimes\Gamma_{2n})\cdots).$

##  3 Numerical examples

We illustrate next the behavior of some of the previously constructed
$T$-methods on a pair of numerical examples. The first one (the 2-dimensional
Kepler problem) allows one to check preservation properties, whereas the
second (a simple diffusion equation) is used as a test of their relative
performance. In all cases we take as basic scheme $\mathcal{S}_{h}^{[2n]}$ the
4th-order ($n=2$) time-symmetric splitting method

${\mathcal{S}}_{h}^{[4]}=\varphi_{b_{1}h}^{[b]}\circ\,\varphi_{a_{1}h}^{[a]}\circ\,\varphi_{b_{2}h}^{[b]}\circ\,\varphi_{a_{2}h}^{[a]}\circ\,\varphi_{b_{3}h}^{[b]}\circ\,\varphi_{a_{2}h}^{[a]}\circ\,\varphi_{b_{2}h}^{[b]}\circ\,\varphi_{a_{1}h}^{[a]}\circ\,\varphi_{b_{1}h}^{[b]}$
(3.1)

with coefficients

$\displaystyle b_{1}$ $\displaystyle=$ $\displaystyle
0.060078275263542357774-0.060314841253378523039\,i,$ (3.2) $\displaystyle
a_{1}$ $\displaystyle=$ $\displaystyle 0.18596881959910913140,$ $\displaystyle
b_{2}$ $\displaystyle=$ $\displaystyle
0.27021183913361078161+0.15290393229116195895\,i,$ $\displaystyle a_{2}$
$\displaystyle=$ $\displaystyle 0.31403118040089086860,$ $\displaystyle b_{3}$
$\displaystyle=$ $\displaystyle
0.33941977120569372122-0.18517818207556687181\,i,$

previously considered in [8]. This integrator is intended for Eq. (1.1) when
$f$ can be decomposed as $f(x)=f_{a}(x)+f_{b}(x)$ in such a way that each sub-
problem

$\dot{x}=f_{a}(x),\qquad\qquad\dot{x}=f_{b}(x),$

with $x(0)=x_{0}$, has solution $x(t)=\varphi_{t}^{[a]}(x_{0})$, and
$x(t)=\varphi_{t}^{[b]}(x_{0})$, respectively.

The implementation of all the integrators has been done in Python 3.7 running
on Debian GNU/Linux 10 and the operations with complex arithmetics have been
coded using the complex class of the numpy library.

#### Kepler problem.

The Hamiltonian function for the planar two-body problem reads

$H(q,p)=T(p)+V(q)=\frac{1}{2}p^{T}p-\mu\frac{1}{r}.$ (3.3)

Here $q=(q_{1},q_{2})$, $p=(p_{1},p_{2})$, $r=\|q\|$, $\mu=GM$, $G$ is the
gravitational constant and $M$ is the sum of the masses of the two bodies. The
corresponding equations of motion are then

$\dot{q}_{i}=\frac{\partial H}{\partial
p_{i}}=p_{i},\qquad\qquad\dot{p}_{i}=-\frac{\partial H}{\partial
q_{i}}=-\mu\frac{q_{i}}{r^{3}},\qquad i=1,2.$

Taking $\mu=1$ and initial conditions

$q_{1}(0)=1-e,\quad q_{2}(0)=0,\quad p_{1}(0)=0,\quad
p_{2}(0)=\sqrt{\frac{1+e}{1-e}},$ (3.4)

the resulting trajectory is an ellipse of eccentricity $0\leq e<1$. In this
case $\varphi_{h}^{[a]}$ (respectively, $\varphi_{h}^{[b]}$) corresponds to
the exact solution obtained by integrating the kinetic energy $T(p)$ (resp.,
potential energy $V(q)$) in (3.3).

We take $e=0.6$, integrate until the final time $t_{f}=20\pi$ with the basic
splitting method ${\mathcal{S}}_{h}^{[4]}$ given by (3.1) and schemes
$T_{h}^{(k)}$, with $k=1,2,3$ for several time steps and then we compute the
average error in energy along the integration interval. Figure 1 (left) shows
this error as a function of the number of evaluations of the basic scheme
${\mathcal{S}}_{h}^{[4]}$. The diagram clearly exhibits the order of
convergence of each method: order 4 for ${\mathcal{S}}_{h}^{[4]}$, and orders
6, 8 and 10 for $T_{h}^{(1)}$, $T_{h}^{(2)}$ and $T_{h}^{(3)}$, respectively.

Figure 1: Left: Relative error in energy vs. number of evaluations of the
basic scheme $\mathcal{S}^{[4]}$ (blue) for the different $T$-methods:
$T^{(1)}$ (orange), $T^{(2)}$ (green) and $T^{(3)}$ (red), in the interval
$t\in[0,20\pi]$. Right: Evolution of this error along the integration when
$t\in[0,2000\pi]$. In this case the step size is chosen so that all schemes
involve the same number of evaluations of the basic method.

In the right panel we show the long-time behavior of the error in energy for
each method when the step size is chosen so that all of them involve the same
computational cost. We see that the error in energy is almost constant for
$t\leq 2000\pi$, as is the case for symplectic integrators. In other words,
the lack of symplecticity at order $h^{12}$ has no effect in this integration
interval. In addition, the scheme $T_{h}^{(3)}$ provides the smaller error.

#### A linear parabolic equation.

Our second example concerns the linear equation in one-dimension

$\frac{\partial}{\partial t}u(x,t)=\frac{\partial^{2}}{\partial
x^{2}}u(x,t)+V(x)u(x,t),\qquad u(x,0)=\sin(2\pi x),$ (3.5)

with periodic boundary conditions in the space domain $[0,1]$. We take
$V(x)=8+4\sin(2\pi x)$ and partition the interval $[0,1]$ into $N$ parts of
length $\Delta x=1/N$, so that the vector
$U=(U_{0},\ldots,U_{N-1})^{T}\in\mathbb{R}^{N}$ is formed, with
$U_{j}=u(x_{j},t)$ and $x_{j}=j/N$, $j=0,1,\ldots,N-1$. If a Fourier spectral
collocation method is used, we end up with the $N$-dimensional linear ODE

$\frac{dU}{dt}=A\,U+B\,U,$ (3.6)

where $B=\mathrm{diag}(V(x_{0}),\ldots,V(x_{N-1}))$ and $A$ is a (full)
differentiation matrix related with the second derivative $\partial_{xx}$. The
splitting here corresponds to solving separately the systems $\dot{U}=A\,U$
and $\dot{U}=B\,U$. Notice that, since $B$ is diagonal, then

$(\mathrm{e}^{hB}U)_{j}=\mathrm{e}^{hV(x_{j})}U_{j}$

and only requires the computation of $N$ multiplications. On the other hand,
$A\,U=\mathcal{F}^{-1}D_{A}\mathcal{F}\,U$, where $\mathcal{F}$ and
$\mathcal{F}^{-1}$ are the forward and backward discrete Fourier transform,
and $D_{A}$ is again diagonal [23]. In consequence,

$\mathrm{e}^{hA}\,U=\mathcal{F}^{-1}\mathrm{e}^{hD_{A}}\mathcal{F}\,U,$

requiring $\mathcal{O}(N\log N)$ operations when the transformation
$\mathcal{F}$ (and its inverse) is computed with the fast Fourier transform
(FFT) algorithm.

Figure 2: Error vs. number of evaluations of the basic method (3.1) obtained
by schemes $T_{h}^{(1)}$ (orange), $T_{h}^{(2)}$ (green) and $T_{h}^{(3)}$
(red). The blue line corresponds to ${\mathcal{S}}_{h}^{[4]}$.

We take $N=128$ and integrate until $t_{f}=1$, where we compute the relative
error $\|U-U_{ex}\|/\|U_{ex}\|$ with each method $T_{h}^{(k)}$, $k=1,2,3$, in
addition to the basic scheme (3.1). The ‘exact’ solution $U_{ex}$ is taken as
the output of the 8th-order composition method P8S15 of [8]. The corresponding
efficiency diagram is shown in Figure 2, where the same notation is used for
the curves depicted. Here also the higher degree integrators provide the best
efficiency.

##  4 Discussion

###  4.1 $T$-methods and $R$-methods

Methods $T_{h}^{(k)}$ have indeed close similarities with the compositions
$R_{h}^{(k)}$ (1.9) previously analyzed in [11]: not only their starting point
is the same (the basic time-symmetric method ${\mathcal{S}}_{h}^{[2n]}$), but
one has in addition $T_{h}^{(1)}=R_{h}^{(1)}$ and also the same coefficients
$\gamma^{[2m]}$ defined in (1.7) enter into their formulation. Finally, they
have the same preservation properties. There is, however, a fundamental
difference: whereas $T$-methods are linear combinations of symmetric-conjugate
compositions only, this is not the case of $R$-methods, and in fact schemes
$R_{h}^{(k)}$ involve a much larger number of compositions. This can be
clearly seen by writing explicitly the expression of $R_{h}^{(2)}$:

$\displaystyle R_{h}^{(2)}$ $\displaystyle=$
$\displaystyle\frac{1}{8}\Big{(}(\gamma^{[2n+2]}\gamma^{[2n]},\gamma^{[2n+2]}\bar{\gamma}^{[2n]},\bar{\gamma}^{[2n+2]}\gamma^{[2n]},\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]})$
$\displaystyle\quad+(\gamma^{[2n+2]}\bar{\gamma}^{[2n]},\gamma^{[2n+2]}\gamma^{[2n]},\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]},\bar{\gamma}^{[2n+2]}\gamma^{[2n]})$
$\displaystyle\quad+(\gamma^{[2n+2]}\gamma^{[2n]},\gamma^{[2n+2]}\bar{\gamma}^{[2n]},\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]},\bar{\gamma}^{[2n+2]}\gamma^{[2n]})$
$\displaystyle\quad+(\gamma^{[2n+2]}\bar{\gamma}^{[2n]},\gamma^{[2n+2]}\gamma^{[2n]},\bar{\gamma}^{[2n+2]}\gamma^{[2n]},\bar{\gamma}^{[2n+2]}\bar{\gamma}^{[2n]})$
$\displaystyle\;\;+\mathrm{c.c.}\Big{)},$

whereas $R_{h}^{(3)}$ is the sum of 64 compositions containing 8 basic schemes
with weights $\gamma^{[2n+4]}\gamma^{[2n+2]}\gamma^{[2n]}$,
$\bar{\gamma}^{[2n+4]}\gamma^{[2n+2]}\gamma^{[2n]}$, etc. plus their complex
conjugate divided by 128. In general, $R_{h}^{(k)}$ involves the sum of
$2^{2^{k}-2}$ compositions of $2^{k}$ appropriately weighted basic schemes:

$R_{h}^{(k)}=\frac{1}{2^{2^{k}-1}}\sum_{j=1}^{2^{2^{k}-2}}\Big{(}(\alpha_{j_{2^{k}}},\ldots,\alpha_{j_{1}})+\mbox{c.c.}\Big{)},$

where $\alpha_{j_{i}}$ are products of the $k$ coefficients
$\gamma^{[2n]},\ldots,\gamma^{[2(n+k-1)]}$ and their complex conjugate. This
should be compared with the $T$-methods: in general, $T_{h}^{(k)}$ involves
the sum of $2^{k-1}$ compositions of $2^{k}$ basic schemes. In either case,
the computation of the complex conjugate part can be avoided just by taking
the real part, with no extra evaluations of ${\mathcal{S}}_{h}^{[2n]}$.

$k$ |  $R$ (explicit) |  $R$ (recursive) |  $T$ (explicit)
---|---|---|---
$1$ | $2$ | $2$ | $2$
$2$ | $16$ | $8$ | $8$
$3$ | $512$ | $32$ | $32$
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
$m$ | $2^{m}\cdot 2^{2^{m}-2}$ | $2^{m}\cdot 2^{m-1}$ | $2^{m}\cdot 2^{m-1}$
Table 1: Number of basic maps ${\mathcal{S}}_{h}^{[2n]}$ necessary to compute
when formulating $R$\- and $T$-methods explicitly or recursively (in the case
of $R$-methods).

These numbers are collected in Table 1, when schemes $R_{h}^{(k)}$ (second
column) and $T_{h}^{(k)}$ (last column) are formulated explicitly. Of course,
a recursive implementation of $R$-methods by applying the procedure (1.9)
turns out to be more efficient. In that case the required computational
effort, measured as the number of basic schemes, is shown in the third column
of the table. Again, in this case we only have to compute the real part in the
last iteration.

In view of the number of basic maps required by the recursive implementation
of $R$-methods and the explicit formulation (1.11) of $T$-methods, it is
natural to ask what are the advantages (if any) of the later schemes with
respect to the former ones. In this respect, one should take into account that
both explicit formulations (1.11) and (4.1) are directly amenable to
parallelization, whereas this is less obvious for the recursion (1.9).

If one has a computer with, say, $2\ell$ threads, it is easy to estimate the
effective number of evaluations of ${\mathcal{S}}_{h}^{[2n]}$ both for $R$\-
and $T$-methods. Thus, for $R_{h}^{(k)}$ one has:

  * •

if $\ell\leq 2^{k}-2$ then the number of evaluations is $2^{k}$;

  * •

if $\ell>2^{k}-2$ then the number of evaluations is $2^{k}\cdot
2^{2^{k}-2-\ell}$,

whereas this number is considerably reduced for schemes $T_{h}^{(k)}$:

  * •

if $\ell\leq k-1$ then the number of evaluations is $2^{k}$;

  * •

if $\ell>k-1$ then the number of evaluations is $2^{k}\cdot 2^{k-1-\ell}$.

In Table 2 we collect these numbers for the first values of $k$ in the
particular case of $2^{2}=4$ and $2^{5}=32$ threads. We see that the
implementation of the explicit expression of the $R$-methods is more
advantageous than the recursive procedure already with a relatively small
number of threads, and that, in any case, $T$-methods require less
computational effort.

| $4$ threads | $32$ threads
---|---|---
$k$ | $R$ | $T$ | $R$ | $T$
$1$ | $2$ | $2$ | $2$ | $2$
$2$ | $4$ | $4$ | $4$ | $4$
$3$ | $128$ | $8$ | $16$ | $8$
Table 2: Effective number of evaluations of the basic map
${\mathcal{S}}_{h}^{[2n]}$ when the corresponding $R$\- and $T$-method is
implemented in parallel with 4 and 32 threads.

To better illustrate this issue, we next compare the efficiency of the
different methods when implemented on a computer able to execute 4 threads
without loss of performance. The corresponding results are displayed in Figure
3 for the Kepler problem (left) and the linear parabolic equation (3.5)
(right). The gain in efficiency of the new schemes is clearly visible.

Figure 3: Efficiency diagram for $R$\- and $T$-methods in a computer accepting
4 threads without loss of performance. Left: Kepler problem. Right: linear
parabolic equation.

Even in the case when one could run the schemes on a machine such that the
effective number of evaluations of both $R_{h}^{(k)}$ and $T_{h}^{(k)}$ is the
same, i.e., $2^{k}$ in both cases, the latter turn out to be more efficient.
This is clearly visible in Figure 4, obtained again by applying the previous
schemes to the Kepler problem (left) and the linear parabolic equation
(right).

Figure 4: Efficiency diagram for $R$\- and $T$-methods when the effective
number of evaluations of the basic scheme is the same ($2^{k}$ in both cases).
Left: Kepler problem. Right: linear parabolic equation.

Finally, it is also illustrative to compare the efficiency of the new
$T$-methods with the standard triple-jump procedure, Eqs. (1.3)-(1.4), both
applied to the same basic scheme (3.1). Thus, in Figure 5 we depict the
results achieved by projecting ${\mathcal{S}}_{h}^{[6]}$,
${\mathcal{S}}_{h}^{[8]}$, and ${\mathcal{S}}_{h}^{[10]}$ at each step,
together with $T_{h}^{(k)}$, $k=1,2,3$ for the Kepler problem with the same
parameters and final time $t_{f}=20\pi$. Here the effective number of
evaluations of the basic scheme has been taken as $2^{k}$ for $T$-methods and
$3^{k}$ for triple-jump. Not surprisingly, the new schemes turn out to be much
more efficient.

Figure 5: Error in energy vs. number of basic maps ${\mathcal{S}}_{h}^{[4]}$
for $T$-methods in comparison with schemes obtained by triple-jump for the
Kepler problem.

###  4.2 Concluding remarks

The standard triple-jump procedure is a popular technique that allows one to
construct numerical integrators for differential equations of arbitrarily high
order by composition of a basic integrator of low order. It has nevertheless
certain limitations: the number of basic maps grows rapidly with the order,
and the main error terms are quite large in comparison with other specially
built integrators. Moreover, they involve some negative coefficients when the
order $r\geq 3$, so that the resulting schemes cannot be used in particular
when the initial value problem (1.1) results from the space discretization of
a parabolic partial differential equation involving the Laplace operator. In
this context it is quite natural to explore whether it is still possible using
the triple-jump technique (1.3), but with the complex coefficients furnished
by (1.5) as long as their real part is positive. It has been established that
this is indeed the case, although once again they require an exceedingly large
number of basic methods. For this reason, other alternatives for constructing
high-order composition methods have also been proposed [8, 12, 17]. Among
them, the class of schemes (1.8) possess some special features: starting from
a time-symmetric basic scheme ${\mathcal{S}}_{h}^{[2n]}$ of order $2n$, it is
possible to construct recursively methods of order $2n+2k$, $k=1,2,\ldots$
that are still time-symmetric up to order $4n+3$. Moreover, if the
differential equation in (1.1) has some qualitative properties (such as
symplecticity or volume preservation) then these properties are still shared
by the numerical solution up to order $4n+3$ [11].

Methods (1.9) are based on the simple symmetric-conjugate composition (1.6).
As shown in [7], symmetric-conjugate composition methods still possess
remarkable preservation properties when projected on the real axis at each
integration step, and so it makes sense to consider more general linear
combinations of methods within this class. The corresponding analysis has been
carried out here, and as a result we have built a new class of schemes that
essentially have the same preservation properties as methods (1.9), but
requiring a much reduced computational cost. In addition, these methods are
particularly well suited for their parallel implementation. The examples
included show a significant improvement in efficiency with respect to schemes
(1.9) and those obtained by applying the triple-jump procedure.

### Acknowledgements

This work has been funded by Ministerio de Ciencia e Innovación (Spain)
through project PID2019-104927GB-C21 (AEI/FEDER, UE) and by Universitat Jaume
I (UJI-B2019-17). A.E.-T. has been additionally supported by the predoctoral
contract BES-2017-079697 (Spain).

## References

  * [1] V. Arnold, Mathematical Methods of Classical Mechanics, Springer-Verlag, Second ed., 1989.
  * [2] A. Aubry and P. Chartier, Pseudo-symplectic Runge–Kutta methods, BIT Num. Math., 38 (1998), pp. 439–461.
  * [3] A. Bandrauk and H. Shen, Improved exponential split operator method for solving the time-dependent Schrödinger equation, Chem. Phys. Lett., 176 (1991), pp. 428–432.
  * [4] S. Blanes, F. Casas, and J. Ros, Extrapolation of symplectic integrators, Celest. Mech. & Dyn. Astr., 75 (1999), pp. 149–161.
  * [5] S. Blanes and F. Casas, On the necessity of negative coefficients for operator splitting schemes of order higher than two, Appl. Numer. Math., 54 (2005), pp. 23–37.
  * [6] S. Blanes and F. Casas, A Concise Introduction to Geometric Numerical Integration, CRC Press, 2016.
  * [7] S. Blanes, F. Casas, P. Chartier, and A. Escorihuela-Tomàs, On symmetric-conjugate composition methods in the numerical integration of differential equations, Tech. Rep. 2101.04100, arXiv, 2021.
  * [8] S. Blanes, F. Casas, P. Chartier, and A. Murua, Optimized high-order splitting methods for some classes of parabolic equations, Math. Comput., 82 (2013), pp. 1559–1576.
  * [9] S. Blanes, F. Casas, and A. Murua, Splitting and composition methods in the numerical integration of differential equations, Bol. Soc. Esp. Mat. Apl., 45 (2008), pp. 89–145.
  * [10] S. Blanes, F. Casas, and A. Murua, Splitting methods with complex coefficients, Bol. Soc. Esp. Mat. Apl., 50 (2010), pp. 47–61.
  * [11] F. Casas, P. Chartier, A. Escorihuela-Tomàs, and Y. Zhang, Compositions of pseudo-symmetric integrators with complex coefficients for the numerical integration of differential equations, J. Comput. Appl. Math., 381 (2021), p. 113006.
  * [12] F. Castella, P. Chartier, S. Descombes, and G. Vilmart, Splitting methods with complex times for parabolic equations, BIT Numer. Math., 49 (2009), pp. 487–508.
  * [13] J. Chambers, Symplectic integrators with complex time steps, Astron. J., 126 (2003), pp. 1119–1126.
  * [14] R. Chan and A. Murua, Extrapolation of symplectic methods for Hamiltonian problems, Appl. Numer. Math., 34 (2000), pp. 189–205.
  * [15] D. Goldman and T. Kaper, $n$th-order operator splitting schemes and nonreversible systems, SIAM J. Numer. Anal., 33 (1996), pp. 349–367.
  * [16] E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations, Springer-Verlag, Second ed., 2006.
  * [17] E. Hansen and A. Ostermann, High order splitting methods for analytic semigroups exist, BIT Numer. Math., 49 (2009), pp. 527–542.
  * [18] R. McLachlan and R. Quispel, Splitting methods, Acta Numerica, 11 (2002), pp. 341–434.
  * [19] J. Sanz-Serna and M. Calvo, Numerical Hamiltonian Problems, Chapman & Hall, 1994.
  * [20] Q. Sheng, Solving linear partial differential equations by exponential splitting, IMA J. Numer. Anal., 9 (1989), pp. 199–212.
  * [21] M. Suzuki, Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations, Phys. Lett. A, 146 (1990), pp. 319–323.
  * [22] M. Suzuki, General theory of fractal path integrals with applications to many-body theories and statistical physics, J. Math. Phys., 32 (1991), pp. 400–407.
  * [23] L.N. Trefethen, Spectral Methods in MATLAB, SIAM, 2000.
  * [24] H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A, 150 (1990), pp. 262–268.