Runge–Kutta methods determined from extended phase space methods for Hamiltonian systems

Many commonly used numerical methods for the time integration of differential equations can be expanded in B-series which elucidate their geometric and numerical properties [2, 10]. However, symplectic integrators with a B-series are implicit, the implicit midpoint rule being a central example [5]. Explicit symplectic integrators exist for some systems, such as separable classical mechanical systems [13]. To avoid this restriction, Pihajoki [15] introduced extended phase space methods: a new Hamiltonian, defined on the product of two copies of the original phase space, is constructed, that is amenable to explicit symplectic integration in the extended phase space.

In place of the Hamiltonian system $X_{H}$ associated with the Hamiltonian $H$ and canonical symplectic form $\omega$,

($(q,p)\in\mathbb{R}^{2d}$), Pihajoki considered the extended system

($(q,x,p,y)\in\mathbb{R}^{4d}$) with initial condition

such that the solution obeys $q(t)=x(t)$ and $p(t)=y(t)$ for all $t$ and $(q(t),p(t))$ satisfies the original system (1) with initial condition $(q(0),p(0))=(q_{0},p_{0})$. The extended system is Hamiltonian with extended Hamiltonian $\hat{H}=\hat{H}_{A}+\hat{B}_{B}$, $\hat{H}_{A}=H(x,p)$, $\hat{H}_{B}=H(q,y)$ and symplectic form $\hat{\omega}:=dq\wedge dp+dx\wedge dy$.

As $x$ and $p$ (resp. $q$ and $y$) commute, the flow $\exp\big{(}tX_{\hat{H}_{A}}\big{)}$ of Hamilton’s equations for $\hat{H}_{A}$ is given explicitly by Euler’s method:

(and analogously for $\hat{H}_{B}$). The integrator

is therefore explicit, second order, and preserves the extended symplectic form $\hat{\omega}$.

Two key issues, however, immediately arise: the ‘duplicate’ point $(x,y)$ may move away from the ‘base’ point $(q,p)$; and it is not clear how symplecticity in the extended phase space is advantageous in the original phase space. Tao [16] addressed the first point by adding a coupling term $\frac{1}{2}\alpha(\|x-q\|^{2}+\|y-p\|^{2})$ to the extended Hamiltonian, finding that this could suppress the growth of $(q-x,p-y)$. Other authors have addressed the second point by projecting the solution to the original phase space in different ways. Two of these, the symmetric projection method of Ohsawa [14] and the midpoint projection method of Luo et al. [7], are the subject of this paper.

The midpoint projection method, considered in Section 2, is shown to be equivalent to an explicit Runge–Kutta method that is pseudosymplectic (that is, approximately symplectic) and pseudosymmetric up to surprisingly high order – order 5 for the leapfrog-based method of classical order 2, and order 9 for the methods of classical order 4. We suggest that these properties account for the methods’ good performance in astrophysical applications.

The symmetric projection methods, considered in Section 3, are shown to be equivalent to monoimplicit symplectic Runge–Kutta methods, revealing their affine equivariance and generality.

Luo et al. [7] composed such an extended phase space integrator with the midpoint projection¹¹1Called the midpoint permutation in [7].

to yield an explicit integrator on the original phase space. These methods were called ‘symplectic-like’ ‘because they, like standard implicit symplectic integrators, show no drift in the energy error’ [7]. This lack of energy drift has been observed in many astrophysical simulations with nonseparable Hamiltonians without explanation, and the method has become quite popular [6, 18]. The order can be increased using composition methods [13]. The following result accounts for the greatly reduced energy drift.

Pseudosymplecticity was introduced by Aubry and Chartier [1], who derive various explicit pseudosymplectic Runge–Kutta methods; the second order method above appears there as a member of a 1-parameter family of 3-stage methods of pseudosymplectic order 4, it being the unique member of that family that is pseudosymplectic of order 5.

The numerical examples in the astrophysics literature do not show energy drift. However, the potential drift effect is rather small and we have confirmed numerically in planar Hamiltonian systems that the energy does drift proportionally to $h^{5}t$ resp. $h^{9}t$ for the 2nd resp. 4th order methods above. Symmetry also affects time integration and can moderate energy drift [12], so it is possible that pseudosymmetry is having some positive effect as well. The 3-stage 4th order method is known to have large error constants and it is possible that the 5-stage method given here, or some other explicit pseudosymplectic method, may be have advantages in these applications. On the other hand, energy behaviour is not the only manifestation of symplecticity and the choice of a symplectic vs a pseudosymplectic integrator may depend on the application.

Ohsawa [14] defined the extended phase space integrator with symmetric projection as follows. Let

Given an extended phase space integrator $\Phi_{h}\colon\mathbb{R}^{4d}\to\mathbb{R}^{4d}$ and initial condition $z_{0}=(q_{0},p_{0})\in\mathbb{R}^{2d}$, the integrator is the map $z_{0}\mapsto z_{1}$ defined by

1. 1.

   $\zeta_{0}:=(q_{0},q_{0},p_{0},p_{0})$

2. 2.

   Find $\mu\in\mathbb{R}^{2d}$ such that $\Phi_{h}(\zeta_{0}+D^{\top}\mu)+D^{\top}\mu\in\mathcal{N}$

3. 3.

   $\hat{\zeta}_{0}:=\zeta_{0}+D^{\top}\mu$

4. 4.

   $\hat{\zeta}_{1}:=\Phi_{h}(\hat{\zeta}_{0})$

5. 5.

   $\zeta_{1}=(q_{1},q_{1},p_{1},p_{1}):=\hat{\zeta}_{1}+D^{\top}\mu$

6. 6.

   $z_{1}=(q_{1},p_{1})$

The method is well-defined, symmetric, and symplectic [14]. Jayawardana and Ohsawa [8] further show that the method preserves arbitrary quadratic invariants. Together these results give a strong indication that the method may be a symplectic B-series method. We show below that this is true and that it is in fact equivalent to a monoimplicit Runge–Kutta method, a class introduced by Cash [3] that have only a single implicit stage. (The Simpson–AVF method $z_{1}=z_{0}+\frac{1}{6}h(f(z_{0})+4f((z_{0}+z_{1})/2)+f(z_{1}))$ is an example [4].) In a monoimplicit Runge–Kutta method the $A$ matrix of the Butcher tableau takes the form of a rank one matrix plus a matrix of zero spectral radius.

The extended phase space integrator with symmetric projection is also an instance of the class of extended Runge–Kutta methods that we now define.

Recall the central result of Munthe-Kaas et al. [9]: Let $\Phi=\{\Phi_{n}\}_{n\in\mathbb{N}}$ be an integration method, defined for all vector fields on all dimensions $n$. Then $\Phi$ is a B-series method if and only if the property of affine equivariance is fulfilled: if $a(x):=Ax+b$ is an affine map from $\mathbb{R}^{m}$ to $\mathbb{R}^{n}$, $f$ a vector field on $\mathbb{R}^{m}$, and $g$ a vector field on $\mathbb{R}^{n}$ such that $g(Ax+b)=Af(x)$, then $a\circ\Phi_{m}(f)=\Phi_{n}(f)\circ a$.

The results here have been established by direct calculation, but the methods appear quite natural and intrinsically defined. One could seek direct geometric proofs, not relying on B-series, of the pseudosymplecticity and pseudosymmetry orders for the midpoint projection method when the extended method is an arbitrary symplectic integrator. Likewise, for the symmetric projection method, the diagonal $\mathcal{N}$ is a symplectic subspace of the extended symplectic space, suggesting an approach using the symplectic geometry of constraints [11]. The question mentioned earlier, of the best pseudosymplectic method to use in applications, and any potential issues arising from nonsymplecticity, should be examined. Finally, we suggest determining the entire set of monoimplicit symplectic Runge–Kutta methods and their relative merits.