Splitting and composition methods with embedded error estimators

TL;DR

This paper introduces new local error estimators for splitting and composition methods that enable efficient adaptive step sizing with minimal additional computational cost, demonstrated through numerical examples.

Contribution

The authors develop novel local error estimators based on linear combinations of intermediate stages, enhancing adaptive step size control in splitting and composition methods.

Findings

Estimators are computationally inexpensive to evaluate.

Adaptive step sizing improves integration efficiency.

Numerical tests confirm estimator effectiveness.

Abstract

We propose new local error estimators for splitting and composition methods. They are based on the construction of lower order schemes obtained at each step as a linear combination of the intermediate stages of the integrator, so that the additional computational cost required for their evaluation is almost insignificant. These estimators can be subsequently used to adapt the step size along the integration. Numerical examples show the efficiency of the procedure.

Tables

Table 1. Table 1: Number of order conditions, in additional to the trivial one for w 0 subscript 𝑤 0 w_{0} , required by a linear combination of intermediate outputs to achieve order ℓ ℓ \ell for symmetric compositions of 2nd-order symmetric schemes ( 𝒮 ​ 𝒮 𝒮 𝒮 \mathcal{SS} ), compositions of a first order method with its adjoint ( 28 ) (method-adjoint) and a splitting method ( 36 ) (splitting).

Order $\ell$	1	2	3	4	5	6
$𝒮𝒮$	1	2	4	7	12	20
method-adjoint	1	3	7	15	31	63
splitting	2	6	14	30	62	126

Equations

$$\dot{x}=f(x),x_0=x(0)\in {\mathbb{R}}^D$$

$$L_{f}G(x)=\sum _{i=1}^d{f}_i(x)\frac{\partial G}{\partial x_i},$$

$$L_f=\sum _{i=1}^d{f}_i\frac{\partial }{\partial x_i}.$$

$${\Phi }^t[G](x)=(G\circ {\phi }_t)(x).$$

$$G({\phi }_t(x_0))=\sum _{k\ge 0}\frac{t^k}{k!}(L_f^{k}G)(x_0)\equiv \exp (t{L}_f)[G](x_0),$$

$${\Phi }^t[G](x)=\exp (t{L}_f)[G](x)\equiv \exp (tF)[G](x),$$

$${\phi }_t(x_0)=\exp (tF)[\mathrm{Id}](x_0).$$

$$X(h)=I+\sum _{n\ge 1}{h}^n{X}_n,$$

$$X_n[G](x)={\frac{1}{n!}\frac{d^n}{d{h}^n}|}_{h=0}G({\chi }_h(x)),$$

$$(G\circ {\chi }_h)(x)=X(h)[G](x)=\exp (Y(h))[G](x),$$

$$Y(h)=\sum _{n\ge 1}{h}^n{Y}_n=\log (X(h)).$$

$$Y_1=F,\text{mbox}and{Y}_n=0\text{mbox}for2\le n\le r.$$

$$\exp(Y(h))=\exp\big{(}hF+h^{2}Y_{2}+h^{3}Y_{3}+\mathcal{O}(h^{4})\big{)},$$

$$(G\circ {\chi }_h^{\ast })(x)=\exp (-Y(-h))[G](x)$$

$$\exp(-Y(-h))=\exp\big{(}hF-h^{2}Y_{2}+h^{3}Y_{3}+\mathcal{O}(h^{4})\big{)}.$$

$${\psi }_h={\mathcal{S}}_{h{\alpha }_s}^{[2]}\circ \cdots \circ {\mathcal{S}}_{h{\alpha }_2}^{[2]}\circ {\mathcal{S}}_{h{\alpha }_1}^{[2]}.$$

$${\alpha }_j={\alpha }_{s-j+1},\text{mbox}forallj$$

$$\Psi (h)=\exp (Y(h{\alpha }_1))\exp (Y(h{\alpha }_2))\cdots \exp (Y(h{\alpha }_2))\exp (Y(h{\alpha }_1)),$$

$$\exp(Y(h\alpha_{k}))=\exp\big{(}h\alpha_{k}F+h^{3}\alpha_{k}^{3}Y_{3}+h^{5}\alpha_{k}^{5}Y_{5}+\mathcal{O}(h^{7})\big{)}$$

$$\Psi (h)=\exp (hF+\mathcal{O}(h^{r+1})),$$

$$\sum _{j=1}^s{\alpha }_j=1,\sum _{j=1}^s{\alpha }_j^3=0$$

$$\sum _{j=1}^s{\alpha }_j^5=0,\sum _{j=1}^s{\alpha }_j^3\left({\left(\sum _{\ell =1}^{j-1}{\alpha }_{\ell }\right)}^2+{\alpha }_j\sum _{\ell =1}^{j-1}{\alpha }_{\ell }\right)=0.$$

$$x_{n,k}={\mathcal{S}}_{h{\alpha }_k}^{[2]}\circ \cdots \circ {\mathcal{S}}_{h{\alpha }_1}^{[2]}{x}_n,k=1,\dots ,s-1,$$

$${\tilde{x}}_{n+1}=\sum _{k=0}^{s-1}{w}_k{x}_{n,k}$$

$$w_s=1,w_k=0,k=0,1,\dots ,s-1$$

$$w_{0}I+\sum _{k=1}^{s-1}{w}_k\prod _{i=1}^k\exp (Y(h{\alpha }_i))=\exp (hF)+\mathcal{O}(h^{\ell +1}),$$

$$\exp (hF)=I+hF+\frac{h^2}{2}{F}^2+\frac{h^3}{3!}{F}^3+\frac{h^4}{4!}{F}^4+\mathcal{O}(h^5)$$

$$w_{0}I+\sum _{k=1}^{s-1}{w}_k\prod _{i=1}^k\exp (Y(h{\alpha }_i))=f_{0}I+h{f}_{1}F+\frac{h^2}{2}{f}_2{F}^2$$

$$+h^3\left(\frac{1}{3!}{f}_{3,1}{F}^3+f_{3,2}{Y}_3\right)+h^4\left(\frac{1}{4!}{f}_{4,1}{F}^4+\frac{1}{2}{f}_{4,2}F{Y}_3+\frac{1}{2}{f}_{4,3}{Y}_{3}F\right)+\mathcal{O}(h^5),$$

$$f_0\equiv w_0+\sum _{k=1}^{s-1}{w}_k=1$$

$$f_1\equiv \sum _{k=1}^{s-1}{w}_k\sum _{j=1}^k{\alpha }_j=1$$

$$\displaystyle f_{2}\equiv\sum_{k=1}^{s-1}w_{k}\,\Big{(}\sum_{j=1}^{k}\alpha_{j}\Big{)}^{2}=1$$

Full text

1 Introduction

Splitting and composition methods are of particular interest in the numerical integration of differential equations when the vector field is separable into solvable parts or when a low order basic method is known, and the goal is to construct higher order schemes by composing the basic method with fractional time steps [25, 26].

Although integrators of this class have a long history in numerical mathematics and have been applied, sometimes with different names, in many different contexts (partial differential equations [32], quantum statistical mechanics [34], chemical physics [16, 18], molecular dynamics [36], celestial mechanics [11, 23], etc.), it has been with the advent of the so-called Geometric Numerical Integration that the interest in splitting and composition has revived and new and very efficient schemes have been designed in the simulation of physical systems. The goal in Geometric Numerical Integration is to construct schemes in such a way that the numerical approximation shares with the exact solution many of its relevant qualitative (very often, geometrical) properties, such as symplecticity, unitarity, orthogonality, etc. [5, 19]. If the basic method possesses (some of) these geometric properties, so do the schemes obtained by composing them. In addition, when they are used with a constant time step, they show a more favorable error growth behavior than standard integrators, especially in long term integrations. Symplectic integration schemes for Hamiltonian dynamical systems constitute a classical example of geometric numerical integrators [30].

Even in problems where no qualitative properties have to be preserved and/or only short time integrations are required, splitting and composition methods have shown to be an excellent option (see e.g. [17] and references therein), even when compared with other standard integrators.

As is well known, some of the most popular and efficient standard schemes are embedded methods: the numerical procedure contains, besides the numerical approximation $x_{n}$, a second approximation $\widetilde{x}_{n}$ (usually of a lower order) obtained from intermediate outputs, so that the difference is used as an estimate of the local error for the less precise result and can subsequently be used for step size control [20]. Well known examples in this area are the class of high order embedded Runge–Kutta methods constructed by Verner [37] (implemented as the DVERK code) and Prince & Dormand [29], giving rise to the code DOP853 [20].

Since splitting and composition methods also provide intermediate outputs when computing the numerical approximation at every step, it seems then natural to analyze whether these intermediate outputs can also be used along the same lines as standard embedded methods to endow the schemes with a step size control. We will see that this is indeed the case as long as the splitting scheme involves a sufficiently large number of stages and, furthermore, we will show how to construct explicitly the lower order approximation $\widetilde{x}_{n}$ from these intermediate outputs at virtually cost free.

It is important to remark that, whereas splitting and composition methods implemented with a constant step size are specially well suited in geometric numerical integration for long time integrations, this is not the case of the variable step size schemes constructed by applying the strategy proposed here [10]. In any case, the second approximation $\widetilde{x}_{n}$ is only used to estimate the local error and this is not propagated along the integration interval.

Of course, the idea of endowing splitting methods with a local error estimator is not new. We can mention in particular references [13, 14], where a embedded splitting method is constructed for the second-order Strang splitting for stiff evolutionary partial differential equations, and [2, 3, 35], where a controller splitting method of order $r+1$ is selected and then an integrator of order $r$ is constructed for which a maximal number of compositions coincide with those of the controller. The methods thus built are then applied for the numerical solution of nonlinear parabolic problems with periodic boundary conditions.

By contrast, the approach we follow here allows one, given a splitting or composition method of order $r$, to construct a second, lower order approximation as a linear combination of the outputs generated at the intermediate stages. This is essentially similar to the procedure presented in [8] for computing cheap approximations to the optimal postprocessor in composition methods with processing, and can be done virtually cost-free. The lower order methods thus designed can be used to endow some of the most popular splitting and composition schemes with a reliable and easy-to-evaluate error estimator [9, 12, 28]

The plan of the paper is the following. In section 2 we briefly summarize the mathematical formalism to be used in the subsequent analysis. Then, in section 3 we proceed to obtain estimators for symmetric compositions of second order basic schemes and of a first order method with its adjoint, whereas an analogous treatment is discussed in section 4. The relationship between composition and splitting methods, together with their respective estimators, is treated in section 5. The new estimators are illustrated in section 6 in comparison with other well established techniques. Finally, section 7 contains some concluding remarks.

2 Flows and Lie derivatives

The analysis of splitting and composition methods can be conveniently carried out with the formalism of Lie derivatives. In that case both the exact flow and the numerical flow corresponding to an integrator, as well as compositions of this integrator, can be associated to the exponential (or products of exponentials) of operators, just as in the linear case, so that the order conditions can be obtained by applying the familiar Baker–Campbell–Hausdorff formula.

To be more specific, given the initial value problem

[TABLE]

with $f:\mathbb{R}^{D}\longrightarrow\mathbb{R}^{D}$ and flow $\varphi_{t}$, we can associate with $f$ the first order differential operator (the Lie derivative) $L_{f}$, whose action on differentiable functions $G:\mathbb{R}^{d}\longrightarrow\mathbb{R}$ is (see [1, Chap. 8])

[TABLE]

so that formally

[TABLE]

Moreover, one can also introduce an operator $\Phi^{t}$ acting on functions $G$ as [27]

[TABLE]

Then, the Taylor series of $G(\varphi_{t}(x_{0}))$ at $t=0$ is given by [19, 5]

[TABLE]

and so

[TABLE]

where, for the sake of simplicity in the notation, we write $F\equiv L_{f}$. If we replace $G$ in (5) by the identity map $\mathrm{Id}(x)=x$, we get for the exact solution of (1)

[TABLE]

In the same way as for the exact flow $\varphi_{t}$, we can associate to each numerical integrator for a time step $h$, $\chi_{h}:\mathbb{R}^{d}\longrightarrow\mathbb{R}^{d}$, the operator

[TABLE]

where $I$ denotes the identity operator and each $X_{n}$ acts on smooth functions $G$ as

[TABLE]

so that $X(h)[G](x)=(G\circ\chi_{h})(x)$. It is then possible to write $X(h)$ formally as the exponential of another operator $Y(h)$,

[TABLE]

where

[TABLE]

Clearly, the integrator $\chi_{h}$ is of order $r$ if $\exp(Y(h))=\exp(hF)$ up to terms $h^{r}$, or equivalently, if

[TABLE]

Thus, in particular, if $r=1$, then

[TABLE]

whereas for its adjoint method $\chi^{*}_{h}\equiv\chi_{-h}^{-1}$, one has analogously

[TABLE]

with

[TABLE]

A second-order method ${\mathcal{S}}^{[2]}_{h}$ is (time-)symmetric if and only if $({\mathcal{S}}^{[2]}_{h})^{*}={\mathcal{S}}^{[2]}_{h}$, or equivalently, if its corresponding operator has the form $Y(h)=hF+h^{3}Y_{3}+h^{5}Y_{5}+\mathcal{O}(h^{7})$.

3 Estimators for composition methods

3.1 Composition of symmetric second order methods

Suppose now that, starting with a basic symmetric second order integrator ${\mathcal{S}}^{[2]}_{h}$, we form the composition

[TABLE]

If the coefficients $\alpha_{1},\ldots,\alpha_{s}$ satisfy some requirements (the order conditions), then $\psi_{h}$ provides an approximation of order $r$ to the exact solution. The number of order conditions is considerably reduced for symmetric compositions, i.e.,

[TABLE]

in (11). In that case its associated series of differential operators reads

[TABLE]

where

[TABLE]

is the operator associated with ${\mathcal{S}}^{[2]}_{h}$. By requiring that

[TABLE]

one gets the order conditions to be satisfied by the coefficients $\alpha_{1},\ldots,\alpha_{s}$ in the composition (11). Up to order $r=6$ these conditions read explicitly

[TABLE]

Notice that, when computing the numerical approximation $x_{n+1}=\psi_{h}(x_{n})\approx x(t_{n+1})=x(t_{n}+h)$ with (11), the procedure also provides $s-1$ intermediate outputs in addition to $x_{n}$, i.e.,

[TABLE]

and the question we pose is whether one can obtain another approximation $\widetilde{x}_{n+1}$ of $x(t_{n+1})$ by a linear combination

[TABLE]

of these intermediate values $x_{n,k}$, with $x_{n,0}=x_{n}$. It turns out that this is indeed possible, but the highest order of approximation that can be achieved in this way depends on the number of intermediate stages $s$. The procedure is similar to the technique used in [7, 8] to construct cheap postprocessors for composition methods with processing.

One should note that $w_{s}$ is not included in the linear combination (15). Otherwise, only the trivial solution

[TABLE]

is obtained.

Our goal is then to find coefficients $w_{k}$ so that, given a number of stages $s$, the linear combination (15) is an approximation to $x(t_{n+1})$ of order $\ell$, or equivalently,

[TABLE]

where $\exp(Y(h\alpha_{k}))$ is given by (13) and $\ell$ is as large as possible. Since a linear combination of exponential operator is not, in general, a exponential operator, the conditions to be satisfied by $w_{k}$ can be derived by expanding both terms in (16) in powers of $h$ and equating their respective coefficients. Thus, in particular, up to order $\ell=4$, one has explicitly

[TABLE]

and

[TABLE]

whence the following system of linear equations results:

[TABLE]

Notice that the first equation is trivially solved in $w_{0}$, so to achieve an approximation $\widetilde{x}_{n+1}$ of order 4, we have to verify 7 linear equations. More generally, the total number of equations (in addition to the trivial one) required to achieve a given order $\ell$ is collected in Table 1 for orders $\ell=1,\ldots 6$. Strictly speaking, this number is the sum of the dimensions $m_{k}$, $k\geq 1$, of the subspaces $\mathcal{A}_{k}$ of the universal enveloping algebra $\mathcal{A}=\bigoplus_{k\geq 0}\mathcal{A}_{k}$ associated to the graded Lie algebra of operators corresponding to the composition method, with $\mathcal{A}_{0}=\mbox{span}(I)$ [8].

Next we analyze in detail the construction of numerical schemes of orders 3, 4 and 5 within this approach to be used as error estimators for symmetric compositions of the form (11).

Third-order estimators.

Only the first five equations in (17) have to be satisfied to get order three. This can be achieved if the composition (11) has at least $s=5$. For $s=5$, when the symmetry of the coefficients (12) (i.e., $\alpha_{5}=\alpha_{1}$, $\alpha_{4}=\alpha_{2}$) and the order conditions of a 4th-order composition (i.e., equations in the first line of (14)) are taken into account, then the unique solution of the system is given by

[TABLE]

so that $w_{0}=-1$. A popular (and efficient) 4th-order composition method within this class is the one devised by Suzuki [33], with coefficients

[TABLE]

so that its third-order estimator reads

[TABLE]

Another widely used 4th-order method involving $s=7$ stages is due to McLachlan [24], with coefficients

[TABLE]

Its corresponding estimator now involves a free parameter, which can be taken to be $w_{3}$, and reads

[TABLE]

Here

[TABLE]

with $g_{1}=\alpha_{1}$, $g_{i}=g_{i-1}+\alpha_{i}$, $i=2,3$.

The same strategy can also be applied to the popular 4th-order 3-stage Yoshida’s method [38]

[TABLE]

with

[TABLE]

which is known to lead to small errors when complex coefficients are taken [6]. Since only three intermediate outputs per step are available, one needs at least two steps of it as if it were one single method, i.e., one can take as integrator the composition

[TABLE]

In this case the corresponding estimator reads

[TABLE]

with

[TABLE]

We can adopt the terminology of embedded Runge–Kutta methods [20] and denote the previous compositions with their respective estimators as methods of order 4(3).

Compositions of order 6(4).

To get linear combinations (15) of order four one has to solve the whole set of equations (17). Although in principle this would require $s=8$ stages, it turns out that if the underlying time-symmetric composition (11) satisfies the order conditions up to order 6 given by (14) with the minimum number of stages ($s=7$), one gets a unique solution of the form

[TABLE]

where $w_{i}$ can be expressed analytically in terms of the $\alpha_{i}$ coefficients of the composition. For the particular method found by Yoshida [38], with coefficients

[TABLE]

one has

[TABLE]

The same strategy can be applied of course if 6th-order compositions with more stages are considered. For instance, we have found an estimator within this class for the symmetric method proposed by Kahan & Li [21], with $s=9$ stages.

Compositions of order 6(5) and 8(5).

A system of 13 linear equations has to be solved for getting an estimator of order five. Although not all of them are independent when the time-symmetry and the order conditions for the underlying composition are introduced, at least $s=11$ stages are necessary. Starting from the 6th-order symmetric composition obtained by Sofroniou & Spaletta [31] with coefficients

[TABLE]

there is just one set of coefficients satisfying all the order conditions. The resulting method of order 6(5) is of the form

[TABLE]

with

[TABLE]

The same strategy can be applied to compositions (11) of order 8. A well known example within this class is the symmetric method proposed by Kahan & Li [21] with $s=17$ and coefficients

[TABLE]

the estimator reads

[TABLE]

with

[TABLE]

The DOP853 algorithm based on a 12-stage RK8(6) method by Dormand & Prince (announced but not published in [15]), where the embedded 6th-order method is replaced by a pair of embedded methods of order five and three by Hairer & Wanner [20]), is one of the most efficient schemes within this framework. In comparison, the previous composition method involves more stages, but on the other hand does not require to keep up to 12 vectors in memory.

As a matter of fact, we can apply the same strategy to the 8th-order composition method considered here and construct a second estimator of order 3 to avoid any possible over-estimation of the error. One possible 3th-order estimator is given by

[TABLE]

with $w_{1}$, $w_{7}$ verifying

[TABLE]

where $g_{1}=\alpha_{1},\ g_{7}=\alpha_{1}+\cdots+\alpha_{7}$, i.e.

[TABLE]

We then have two error estimators for the scheme (11) with coefficients (25),

[TABLE]

Applying now the same strategy as in [20], we consider

[TABLE]

as an error estimator that behaves asymptotically like the global error of the method.

Notice that we can obtain error estimators for other composition schemes in a similar way. For example, at order eight one can find in the literature methods with up to 21 stages [19, 21, 31], and their relative performance depend on the particular problem to solve as well as on the symmetric second order scheme used as the basic scheme for the composition.

3.2 Composition of a first order method with its adjoint

Higher order methods can also be obtained by composing a first order basic method $\chi_{h}$ and its adjoint $\chi^{*}_{h}=\chi_{-h}^{-1}$,

[TABLE]

with appropriately chosen real coefficients $(\alpha_{1},\ldots,\alpha_{2s})$. The associated series of differential operators is of the form

[TABLE]

where $Y(h\alpha_{k})=h\alpha_{k}F+h^{2}\alpha_{k}^{2}Y_{2}+h^{3}\alpha_{k}^{3}Y_{3}+\mathcal{O}(h^{4})$. Again, by requiring that $\Psi(h)=\exp(hF+\mathcal{O}(h^{r+1}))$, one gets the order conditions to be satisfied by the coefficients to achieve order $r$. These order conditions are considerably simplified if $\alpha_{2s-j+1}=\alpha_{j}$ for all $j$. In that case the composition (28) is time-symmetric.

As with symmetric compositions of symmetric second order schemes, here we can also take a linear combination

[TABLE]

of intermediate outputs

[TABLE]

to produce an approximation of order $\ell<r$ to be used as an error estimator for the composition (28). The coefficients $w_{k}$ can be determined by requiring that

[TABLE]

By expanding the product of exponentials we get the number of conditions the $w_{k}$ have to satisfy at a given order in a similar way as with compositions of 2nd-order symmetric methods. This number is collected in Table 1.

In particular, 8 linear equations are required to get a 3rd-order approximation in this way. Since several efficient 4th-order methods of this class with up to 6 stages (or 12 intermediate outputs) are available in the literature, it is in principle possible to get third order estimators for them (even with free parameters for optimization). As an illustration, for the symmetric 4th-order method (28) with $s=6$ and coefficients

[TABLE]

we propose the linear combination (30) with $w_{0}=-1$ and

[TABLE]

4 Estimators for splitting methods

If $f$ in equation (1) can be split as $f=\sum_{i=1}^{m}f^{[i]}$ for certain functions $f^{[i]}:\mathbb{R}^{D}\longrightarrow\mathbb{R}^{D}$, in such a way that the equations

[TABLE]

can be integrated exactly, with solutions $x(h)=\varphi_{h}^{[i]}(x_{0})$ at $t=h$, then the basic first-order method in the composition (28) can be taken simply as

[TABLE]

whereas its adjoint is just the reversed composition

[TABLE]

For $m=2$, i.e., when $f(x)$ is decomposed in just two pieces,

[TABLE]

one could also consider a time-symmetric composition

[TABLE]

with appropriately chosen coefficients $a_{i}$, $b_{i}$ verifying

[TABLE]

to achieve a prescribed order. Here it is also possible to take advantage of the intermediate outputs to construct a lower order approximation which may be used as an error estimator for the integrator (36). In this case it has the form

[TABLE]

with

[TABLE]

As before, the analysis can be carried out with the associated series of differential operators, which in this case reads

[TABLE]

where $A$ and $B$ denote the Lie derivatives corresponding to $f^{[1]}$ and $f^{[2]}$, respectively:

[TABLE]

Analogously, the conditions to be satisfied by the $w_{i}$ are determined by expanding the exponentials in

[TABLE]

The number to achieve a given order is collected in Table 1 (last line).

Now a system of 15 equations have to be satisfied by the coefficients $w_{i}$ in the linear combination (37) to achieve order 3. As in the preceding cases, we can take several efficient splitting methods of the form (36) involving enough intermediate steps and construct estimators for them. In particular, for the 4th-order symmetric splitting scheme designed by Blanes & Moan [9], with 12 intermediate outputs

[TABLE]

and coefficients

[TABLE]

we propose the linear combination

[TABLE]

solving all order conditions with

[TABLE]

Another particularly efficient 4th-order splitting method designed for systems of the form

[TABLE]

when written as a first order system

[TABLE]

corresponds to the composition (38) with

[TABLE]

In this case the estimator has also the form (39) with

[TABLE]

This splitting method, as well as the error estimator, can also be used to integrate in time the Schrödinger equation

[TABLE]

where $m$ is the reduced mass, $\Delta$ is the Laplacian operator and $V(x)$ is the potential. After spatial discretisation one has to solve a linear system of ODEs

[TABLE]

where $A$ corresponds to the spatial discretization of the kinetic part and $B$ to the potential part. Here $B$ is a diagonal matrix in the coordinates space, whereas $A$ is diagonal in the momentum space, so fast Fourier transform (FFT) algorithms ${\cal F}$ can be used to compute the action of a $A$ on a vector, $Au=F^{-1}D_{A}Fu$, with $D_{A}$ a diagonal matrix.

5 Connection between splitting and composition

Splitting and composition methods for system $\dot{x}=f^{[1]}(x)+f^{[2]}(x)$ are closely connected. On the one hand, if ${\mathcal{S}}^{[2]}_{h}=\varphi_{h/2}^{[2]}\circ\varphi_{h}^{[1]}\circ\varphi_{h/2}^{[2]}$ or ${\mathcal{S}}^{[2]}_{h}=\varphi_{h/2}^{[1]}\circ\varphi_{h}^{[2]}\circ\varphi_{h/2}^{[1]}$, then the composition scheme (11) can be written as (36), although the opposite is not true in general. On the other hand, if $\chi_{h}=\varphi^{[2]}_{h}\circ\varphi^{[1]}_{h}$, then $\chi_{h}^{*}=\varphi^{[1]}_{h}\circ\varphi^{[2]}_{h}$ and the composition (28) reads

[TABLE]

Since $\varphi^{[i]}_{h}$ ($i=1,2$) are exact flows, then they verify444This property is not satisfied, in general, if the exact flows are replaced by numerical approximations. $\varphi^{[i]}_{\beta h}\circ\varphi^{[i]}_{\delta h}=\varphi^{[i]}_{(\beta+\delta)h},$ and the method can be rewritten as the splitting scheme

[TABLE]

if $b_{1}=\alpha_{1}$ and

[TABLE]

(with $\alpha_{2s+1}=0$). Conversely, any integrator of the form (45) with $\sum_{i=1}^{s}a_{i}=\sum_{i=1}^{s+1}b_{i}$ can be expressed in the form (28) with $\chi_{h}=\varphi^{[2]}_{h}\circ\varphi^{[1]}_{h}$ and

[TABLE]

with $\alpha_{0}=0$ for consistency. Nevertheless, the intermediate outputs are different in each implementation as well as the number of order conditions for the estimators. In general this number grows faster with the order for splitting methods. Moreover, implementing the splitting scheme $\psi_{h}$ as a composition method is in general more costly because explicitly obtaining the intermediate values requires the computation of additional basic flows. In more detail, suppose we write (45) as a composition:

[TABLE]

Then, for the first intermediate output we have

[TABLE]

However, whereas obviously $\varphi^{[1]}_{(a_{1}-b_{1})h}\circ\varphi^{[1]}_{b_{1}h}=\varphi^{[1]}_{a_{1}h}$, the computational cost of computing $z=\varphi^{[1]}_{b_{1}h}(y)$ and then $\varphi^{[1]}_{(a_{1}-b_{1})h}(z)$ can be in many cases up to twice more costly than directly evaluating $\varphi^{[1]}_{a_{1}h}(y)$.

For example, taking this composition for solving the Schrödinger equation requires the computation of $s$ additional inverse FFTs with respect to the same scheme written as a splitting method. Similarly, taking a composition with the symmetric second order scheme ${\mathcal{S}}^{[2]}_{h}=\varphi_{h/2}^{[2]}\circ\varphi_{h}^{[1]}\circ\varphi_{h/2}^{[2]}$ requires the same number of FFTs as the corresponding splitting composition, but taking instead ${\mathcal{S}}^{[2]}_{h}=\varphi_{h/2}^{[1]}\circ\varphi_{h}^{[2]}\circ\varphi_{h/2}^{[1]}$ as the basic scheme, requires $s$ additional inverse FFTs for the intermediate outputs because $\varphi_{h}^{[1]}$ carries the costly part of the scheme.

A noteworthy exception is the case in which $f^{[1]}$ and $f^{[2]}$ originate from a partitioned ordinary differential equation of the form

[TABLE]

The system can then be written as

[TABLE]

and

[TABLE]

where the same evaluation $g(p_{n})$ is used in both cases.

The algorithm corresponding to the splitting method (45) for the step $(q_{0},p_{0})\mapsto(q_{1},p_{1})$ reads

[TABLE]

so that it can be seen as an explicit partitioned Runge–Kutta method. On the other hand, the composition (28) with (44) leads to the algorithm

[TABLE]

requiring exactly the same evaluations of $f$ and $g$. If in addition $g(p)=p$ (i.e., if we are solving the second order differential equation $\ddot{q}=f(q)$), then the estimator for (47) takes the form (for appropriate choices of the parameters $w_{i}$)

[TABLE]

in a similar way as for embedded Runge–Kutta–Nyström methods. In any case, other choices of $\delta_{i},\gamma_{i}$ can also lead to estimators associated to a given $s$-stage composition scheme [10], and that can not be obtained by taking intermediate outputs.

6 Numerical examples

In this section we analyze the accuracy and reliability of the estimators presented in this work in comparison with other well established schemes for a simple example. Specifically, the methods (and notation) we consider are the following:

• •

$\mathcal{RKN}_{6}43$: The 6-stage 4th-order splitting method (38) for systems of the form (41) with the 3rd-order estimator (43).

• •

$\mathcal{PRK}_{6}43$: The 6-stage 4th-order splitting (38), with the 3rd-order estimator (39) and coefficients given by (40).

• •

$\mathcal{S}_{6}43$: The 6-stage 4th-order method-adjoint symmetric composition (28) with coefficients (31) and 3rd-order estimator (32).

• •

$\mathcal{SS}_{5}43$: The 5-stage 4th-order symmetric composition (11) with coefficients (19) and 3rd-order estimator (20).

• •

$\mathcal{SS}_{11}65$: The 11-stage 6th-order symmetric composition (11) with coefficients (23) and 5th-order estimator (24).

• •

$\mathcal{SS}_{17}853$: The 17-stage 8th-order symmetric composition (11) with coefficients (25) with the 5th- and 3rd-order estimators (26) and (27).

These are compared with:

• •

$\mathrm{eRKN}_{4}43$: the non-symmetric 4-stage 4th-order Runge–Kutta–Nyström (RKN) method with a 3rd-order estimator presented in [10]. This method has an error estimator that is only valid for equations of the form (41), so that it cannot be used in particular for the Schrödinger equation.

• •

$\mathrm{ePRK}_{5}43$: The 5-stage 4th-order splitting method given by the composition

[TABLE]

with the symmetry $b_{6-i}=a_{i},\ i=1,2,\ldots,5$, and the 3rd-order estimator given by a similar composition sharing the first stages555The idea to consider estimators using a second composition sharing some of the stages was first proposed in [22]

[TABLE]

with $\widetilde{a}_{i},\widetilde{b}_{i},\ i=3,4,5$ chosen appropriately. The estimator requires three new evaluations. We take in particular the scheme666The corresponding coefficients are available at http://www.asc.tuwien.ac.at/~winfried/splitting Emb 4/3 AK p, in which case $\widetilde{a}_{3}=0$, so that only two new evaluations are required and the overall cost is taken as 7 evaluations per step.

• •

RK6(5): the well known 8-stage Verner’s method of order 6(5) (see Table 5.4 in [20], page 181) that is implemented in the routine DVERK.

• •

DOP853: the 12-stage embedded Runge–Kutta method of order 8(5) by Dormand & Prince [15] and improved as the routine DOP853 in [20].

Specifically, we consider as a test bench the two-dimensional Kepler problem with Hamiltonian

[TABLE]

Here $q=(q_{1},q_{2}),p=(p_{1},p_{2})$, $\mu=GM$, $G$ is the gravitational constant and $M$ is the sum of the masses of the two bodies. Taking $\mu=1$ and initial conditions

[TABLE]

if $0\leq e<1$, then the total energy is $H=H_{0}=-1/2$, the solution is periodic with period $2\pi$, and the trajectory is an ellipse of eccentricity $e$.

The performance of an embedded Runge–Kuta method depends on the performance of the high order method used to propagate the solution, but also on the accuracy of the lower order one as well as how the error estimator approaches the true error of the high order method. Some times the error estimator is much larger than the true error and the algorithm uses smaller time steps than necessary to reach a given accuracy. Some other times, however, this error can be considerably smaller than the true error (usually due to cancellations because the methods share internal stages) and the algorithm takes longer time steps than required which lead to undesirable large errors.

In this example we integrate with a constant time step and compute the maximum true error

[TABLE]

and the maximum error estimator

[TABLE]

An efficient method should give ${\cal E}_{2}\sim{\cal E}_{1}$, while being both as small as possible at a given computational cost.

The integration is carried out in the time interval $t\in[0,20]$ with a constant time step, and this integration is repeated for different values of the time step and for several values of the eccentricity, in particular for $e=\frac{1}{5},\ \frac{2}{5},\ \frac{3}{5},\ \frac{4}{5}$. This is done first for RK6(5) (or DVERK subroutine) and the composition scheme $\mathcal{SS}_{11}65$.

Figure 1 shows in double logarithmic scale the error ${\cal E}_{1}$ (thin lines) and the estimate ${\cal E}_{2}$ (thick lines) versus the computational cost measured as the number of force evaluations. Dashed lines are obtained with RK6(5), whereas solid lines correspond to $\mathcal{SS}_{11}65$.

We notice from the figure that the composition method is not only more accurate at the same cost (even for such a short time integration) but also the error estimator is much closer to the true error. The error estimator of DVERK is very optimistic: ${\cal E}_{2}$ is much smaller that ${\cal E}_{1}$, especially when the eccentricity takes large values (and thus adjusting the step size is increasingly relevant). The reason lies in the fact that both $x_{n}$ and $\widetilde{x}_{n}$ are computed using very similar procedures, since they share the intermediate stages. This is not the case for the error estimators proposed here, and thus the error ${\cal E}_{2}$ is reasonably close to the true error of the method, even when the coefficients for this specific method are not particularly small.

Next the same numerical experiment is carried out again, but this time with DOP853 and the composition scheme $\mathcal{SS}_{17}853$. Figure 2 shows the results obtained.

We observe that, for this example, the symplectic composition method is as efficient as the 8th-order RK method even for such a short time integration. In addition, our error estimator for the composition method is closer to the true error providing a better error estimator and as a result allowing to choose more appropriate time steps.

Next we compare the results achieved by methods of order 4(3) that are valid for general splitting methods and symmetric-symmetric compositions. This is shown in Figure 3 for eccentricity $e=1/2$ in eq. (53): $\mathcal{PRK}_{6}43$ (dashed lines); $\mathrm{ePRK}_{5}43$ (dot-dashed lines); and $\mathcal{SS}_{5}43$ (solid lines). We observe that the embedded scheme $\mathrm{ePRK}_{5}43$ provides an exceedingly optimistic error estimator as well as a lower performance due to its higher cost per step.

Finally, Figure 4 shows the same results as Figure 3 for the RKN methods of order 4(3) and the composition method-adjoint obtained from the coefficients of the 6-stage RKN method and the relation (5). It provides the same results for the 4th-order method, but different outputs for the estimator. Specifically, we collect the results obtained with $\mathcal{S}_{6}43$ (dashed lines), $\mathrm{eRKN}_{4}43$ (dot-dashed lines), and $\mathcal{RKN}_{6}43$ (solid lines). We observe that the scheme $\mathrm{eRKN}_{4}43$ provides an optimistic error estimator as well as a lower performance.

7 Concluding remarks

In this work we have proposed a procedure to estimate the local error of splitting and composition methods based on the construction of a second lower order integrator by linear combinations of the intermediate outputs of the original scheme. The difference can then be combined with standard strategies of automatic step size control [20] to use the original splitting and composition methods with adaptive step size along the integration. In contrast with other approaches, the proposed strategy does not increase the computational cost of the overall scheme and provides a reliable estimate of the error, so that it can be safely used in problems where keeping the step size constant is not of paramount importance, such as it is the case in certain partial differential equations of evolution. In any event, in that case one should use a very precise discretization in space to guarantee that the main source of error originates when integrating in time.

We should remark in particular the good properties exhibited by the estimator constructed for the 17-stage 8th-order composition scheme (11) with coefficients (25) in comparison with the well known routine DOP853. Taking into account that even more efficient composition methods involving 19 and 21 stages do exist within this class, we conclude that these can constitute a worthwhile alternative for integrating problems when high accuracy is required.

The error estimator proposed here coupled with a variable step size strategy could be most useful for the application of splitting methods for solving the Schrödinger eigenvalue problem with the imaginary time propagation technique, in order to reduce the overall computational cost, as illustrated e.g. in [4].

Although only several representative schemes have been considered, it is clear that the same strategy can be applied to any other splitting and composition method. In particular, we can also construct estimators for the high-order methods with complex coefficients collected in [6] and schemes involving double commutators, such as those presented in [12, 23, 28], as long as they involve a sufficiently large number of intermediate stages to form the required linear combinations.

Acknowledgements

Part of this work was developed during a research stay at the Wolfgang Pauli Institute Vienna; the authors are grateful to the director Norbert Mauser and the staff members for their support and hospitality. The first two authors acknowledge funding by the Ministerio de Economía y Competitividad (Spain) through project MTM2016-77660-P (AEI/FEDER, UE).