Composite symmetric general linear methods (COSY-GLMs) for the long-time integration of reversible Hamiltonian systems

In choosing a numerical method for the long-time integration of reversible Hamiltonian systems one must take into consideration several key factors: order of the method, ability to preserve invariants of the system, and efficiency of the computation. In this paper, 6th-order composite symmetric general linear methods (COSY-GLMs) are constructed using a generalisation of the composition theory associated with Runge–Kutta methods (RKMs). A novel aspect of this approach involves a nonlinear transformation which is used to convert the GLM to a canonical form in which its starting and finishing methods are trivial. Numerical experiments include efficiency comparisons to symmetric diagonally-implicit RKMs, where it is shown that COSY-GLMs of the same order typically require half the number of function evaluations, as well as long-time computations of both separable and non-separable Hamiltonian systems which demonstrate the preservation properties of the new methods.

Hamiltonian system, the class of symmetric second-order multistep methods have been identified as excellent candidates for integration, as they are explicit, high-order, and preserve invariants of the system over long-times [14,Ch. XV], [18]. However, for non-separable Hamiltonian systems, these methods are generally not applicable. Thus, one would usually consider a symmetric (or possibly symplectic) Runge-Kutta method (RKM) as an alternative. It is well-known that structure-preserving RKMs are necessarily implicit, which inevitably limits their long-time efficiency. The impact of implicitness can be mitigated to some extent by considering the class of symmetric diagonally-implicit RKMs (DIRKs). Here, since the methods are essentially compositions of the implicit midpoint rule (IMR) [19], the computational cost amounts to the number of its stages multiplied by the cost of IMR.
Alternatively, one could consider selecting a method from the class of symmetric general linear methods (GLMs). In recent work, it has been shown that these methods can preserve quadratic and Hamiltonian invariants over long-times without suffering from the parasitic instability associated with multistep methods [5][6][7]. Furthermore, these methods can be designed such that they consist of a mixture of implicit and explicit stage equations, suggesting they have the potential to outperform symmetric DIRKs, in terms of cost for a given order.
The choice of method will also depend on its order and computational efficiency. At present, symmetric GLMs are mostly limited to 4th-order [6,17], with the exception of the method presented in [7] that is of 6th-order. Beyond this, the construction of higher-order methods is difficult. In this paper, we present a construction process that is based on the theory of composition for one-step methods (OSMs), such as RKMs, which in the past has been applied to generate composite symmetric (COSY) methods of arbitrarily high order [10,12,16,20,21]. In particular, we generalise the following OSM-composition formula (see e.g. [14,Ch. II.4]) to GLMs: where φ h denotes the OSM and φ * h := φ −1 −h denotes its adjoint. In addition, we also consider the GLM generalisation of the following composition which is frequently used with symmetric methods of even order p ∈ N. Notable examples include the triple jump composition: (1.3) and the : (1.4) The main challenge with constructing a COSY-GLM is understanding how its inputs change between each method evaluation. For OSMs, this is not an issue as there is only one input which approximates the exact solution. However, GLMs possess multiple inputs, each of which takes the form of a B-series. For the special case where these are in Nordsieck form, an appropriate re-scaling of each input is sufficient to guarantee an order increase [17,Ch. 5]. In this paper, we address the issue of general inputs by introducing a nonlinear transformation that puts the GLM into a canonical form. Here, canonical is taken to mean that the resulting method has trivial starting and finishing methods given by its preconsistency vectors u and w H . As a result, its inputs are essentially a Kronecker product of the preconsistency vector u and an approximation to the exact solution. Thus, the method behaves in this respect as if it were a OSM and the theory of composition can be straightforwardly applied to generate high-order COSY-GLMs.
This paper is organised as follows: In Sect. 2 we give an introduction to GLMs. In Sect. 3 we introduce a transformation that puts a GLM into a canonical form. In Sect. 4 we use the canonical transformation to develop composition formulae for GLMs. In Sect. 5, we test our GLM-composition formula by constructing 6th-order COSY-GLMs. Efficiency comparisons are then made between these methods and symmetric DIRKs of the same order. Finally, several long-time Hamiltonian experiments are performed to demonstrate the preservation properties of the methods.

General linear methods
Throughout the paper, we assume numerical methods are applied to the following autonomous, initial value problem (IVP) where f : X → X , and the solution y : R → X is expressed in terms of the flow map ϕ t : X → X and initial data y 0 such that Of particular interest will be Hamiltonian IVPs (see e.g. [ where H : X → R, is the Hamiltonian and X = R 2d , d ∈ N.

Method definition
An r ∈ N-input GLM, for fixed time-step h ∈ R, is written as the map M h : X r → X r acting on inputs y [n] ∈ X r at step n ∈ N 0 such that where ⊗ denotes a Kronecker product, I X is the identity matrix defined on X , s ∈ N denotes the number of stages and its coefficient matrices are denoted by For compactness, we will refer to a GLM using its tableau: A U B V .

Convergence, consistency and stability
The convergence of a GLM is guaranteed if the method is both consistent and stable.

Starting and finishing methods
Inputs to a GLM generally take the form of B-series. Consequently, a starting method is required to generate the starting vector y [0] and a finishing method is required to obtain approximations to the solution y(nh). [3,4] A starting method is defined as the map S h : X → X r , where

Definition 2.2
A finishing method is defined as the map F h : X r → X such that Both S h and F h can be viewed as GLMs with tableaux respectively given by where u, w are as in the definition for preconsistency, and it is assumed that . This is desirable as the finishing method need not perform any additional function evaluations to approximate y(nh).
The numerical method in its entirety may now be expressed as the composition where C(y 0 ) = 0 is a constant vector depending on the method and various derivatives of f evaluated at y 0 .
It is worth mentioning that starting and finishing methods have found application in a variety of areas outside of GLMs. For example, Butcher [2] uses them to achieve an effective order increase in RKMs. Chan and Gorgey [8] apply passive symmetrizers to the Gauss methods, where the symmetrizers are essentially finishing methods. Also, the Gragg smoother used in extrapolation codes (see e.g. [15,Ch. II.9]) is a finishing method that eliminates the leading order parasitic term in the Leapfrog method.

Composition
The composition of two GLMs, M 2 h • M 1 h , can be computed using the following tableau: ⎡

Equivalence
Definition 2.4 [6] Consider a pair of GLMs with coefficient matrices . Then, the two are said to be (T, P)-equivalent if there exists an invertible matrix T ∈ C r ×r and an s × s permutation matrix P such that their coefficient matrices satisfy Permutation-based equivalence arises from the fact that F(Y ) = P P −1 F(Y ) = P F(P −1 Y ). Transformation-based equivalence arises from studying the numerical method as a whole: Notice that under the transformation T , both starting and finishing methods also undergo a transformation, with their tableaux now reading as For a GLM, the direct analogue of this definition is quite restrictive. To see this, consider the tableau of the GLM-adjoint method M * h := M −1 −h : ) to obtain the inverse method. Then, reverse the sign of h to obtain the adjoint method.
Here, we note that the inverse and adjoint methods exist if and only if V −1 exists. Furthermore, unless V is an involution, then a GLM cannot be symmetric according to the OSM definition. However, it is possible that a GLM is equivalent to its adjoint. Definition 2.6 (GLM symmetry [6]) A GLM is said to be (L , P)-symmetric if there exists an r × r involution matrix L, and an s × s symmetric permutation matrix P such that it is (L , P)-equivalent to its adjoint, i.e. if its tableau satisfies In terms of maps, the symmetry condition may be equivalently expressed as An important observation that should be noted is that the adjoint method requires a different set of starting and finishing methods, namely, S −h (y 0 ) and F −h (y [n] ). This choice ensures that the order of the method is preserved, as is summarised by the following lemma.

Symmetric starting and finishing methods
Since symmetry is defined in terms of equivalence to the adjoint method, it follows that there must exist an alternative set of starting and finishing methods that will not affect the order of the method [6].
Definition 2.7 [6] Consider an (L , P)-symmetric GLM with starting and finishing methods, S h and F h , described by the tableaux (2.2). Then, S h and F h are said to be (L , P S )-symmetric if there exists a symmetric permutation matrix P S such that In terms of maps, condition (2.8) can be written as

A canonical form for GLMs
In this section we consider equivalent methods that arise from a nonlinear change of coordinates. As we have seen in Sect. 2.5, any change of coordinates will also transform the corresponding starting and finishing methods. It is shown below that there exists a nonlinear transformation that can convert the starting and finishing methods to a trivial form, i.e. they are given by the preconsistency vectors u and w H .

Definition 3.1
A pth-order GLM is said to be canonical if its starting and finishing methods are given by its preconsistency vectors u and w H .
Canonical methods have the important property that their inputs are independent of h. Thus, we can compose multiple canonical methods of different time-steps provided only the preconsistency vectors agree.

Theorem 3.1 Every pth-order GLM M h , with starting and finishing methods, S h and F h , determined by the tableaux (2.2), is equivalent to a pth-order canonical GLM defined by the composition
Now, consider a nonlinear transformation of the numerical method as a whole, i.e.
Note that the corresponding starting and finishing methods of C h are given by where the tableaux for S h and F h are given by (2.2). Observe that the composition T h • u yields a tableau of the form where we have used U F u = 1 S from (2.3). This agrees with the tableau for S h and thus it follows that S C h (y 0 ) = u ⊗ y 0 . Also observe that w H T −1 h yields a tableau of the form Thus, (C h , u) is of order p, and by Definition 3.1 it follows that C h is a canonical method.
Remark 3.1 In the above theorem, a different notion of equivalence is used than that was introduced in Definition 2.4, i.e. w.r.t. a nonlinear transformation T h . However, note that if T h = T 0 = T for some T ∈ X r ×r then equivalence is defined in the usual sense.
The tableau for the corresponding canonical method of a GLM may be obtained using the tableau composition formula (2.5): Preservation of symmetry In general, performing a nonlinear change of coordinates runs the risk of destroying certain properties of the underlying GLM. In particular, symmetry is not preserved unless the starting and finishing methods are also symmetric.
Proof Since S h and F h are symmetric, this implies that their coefficient matrices satisfy condition (2.8), i.e.

Upon substitution into the tableaux for T h and T
However, Thus, the canonical method C h = LC * h • L is symmetric, as required.

Preservation of parasitism-free behaviour:
In recent work [5,6], it has been shown that structure-preserving GLMs can be designed such that they are free from the influence of parasitism over intervals of length O(h −2 ) [11]. In particular, if their coefficient matrices and preconsistency vectors satisfy the following condition [6]: Corollary 3.2 Suppose that M h satisfies the parasitism-free condition (3.3). Then, C h is also parasitism-free.
Proof From (3.2), we take the expressions for the B, U matrices of C h and insert into the LHS of (3.3) to find where we have used V u = u, w T V = w T from the definition of the preconsistency vectors, and w T BU u = 0 since M h satisfies (3.3). Thus, C h is also parasitism-free.
Example 3.1 It has been shown by Gragg [13] that the Leapfrog method, Here, we observe that the second and third stage equations are equivalent, which implies there is a redundancy in the representation of the method. By removing one of these redundant stages (i.e. combining the second and third columns together, then removing the third row), we obtain the irreducible representation given by the tableau It can now be verified using the symmetry conditions (2.7), that the canonical method is (L C , P C )-symmetric where In addition, we observe that the starting and finishing methods, u and w H , are trivially symmetric, i.e. L C u = u and w H L C = w H . Thus, the Leapfrog method written out in full is given by where w H , C h and u are all symmetric with respect to L C . Finally, we may now apply Theorem 14 of [6], which states that a symmetric GLM of even order p, with symmetric starting and finishing methods, yields a global error expansion in even powers of h.

Composition of canonical methods
Consider a canonical GLM C h with an invertible matrix V . Since inputs to canonical methods are given by their preconsistency vector u we can consider a straightforward generalisation of composition (1.1) to GLMs: It can be shown that the conditions on α 1 , . . . , α k , β 1 , . . . , β k required for an order increase agree with those used for the composition of OSMs.

Theorem 4.1 Suppose the pair (C h , u) is of order p ∈ N and its V -matrix is invertible. Then, the pair (C A h , u) is at least of order p + 1 provided
For canonical methods, S h = S −h = u. Thus, for each j ∈ {1, . . . , k}, we have Composing these expressions, we find that Recursively applying the above result to each C * β j h •C α j h in the order of j = 1, . . . , k we find Thus, if (4.2) and (4.3) are satisfied, it follows that the pair (C A h , u) is at least of order p + 1.
To obtain an adjoint-free composition, i.e. a GLM-generalisation of (1.2), we set β j = 0 for j = 1, . . . , k in (4.1). This choice replaces each C * β j h by V −1 to give Notice here that the final left-hand multiplication by V −1 will not affect the order of the method as V u = u. In other words, if Thus, we define the adjoint-free composition of canonical GLMs as

Corollary 4.1 Let the assumptions of Theorem 4.1 hold. Then, the pair (C B h , u) is at least of order p + 1 provided
Proof This follows from Theorem 4.1 with β 1 = · · · = β s = 0, noting that V u = u. Preservation of symmetry Suppose now that the canonical method is (L , P)symmetric. Without loss of generality, we restrict our attention to compositions of the form (4.4) as, by definition, a symmetric method is similar to its adjoint.

Corollary 4.2 Let C h be an (L , P)-symmetric, canonical GLM. Then, composition
Proof Taking the adjoint of C B h , we find By assumption, α j = α k− j+1 , for j = 1, . . . , k. Thus, this becomes Since C h is symmetric, we have that C h = LC * h • L and LV L = V −1 . Therefore, and the method is symmetric as required.
The above result coupled with the necessity of even order for symmetric methods (cf. Theorem 14 of [6]) implies that the composite method will achieve an increase of two orders, i.e. p → p + 2. Furthermore, this composition can be repeatedly applied to generate canonical COSY-GLMs of arbitrarily high order. Preservation of parasitism-free behaviour Consider now the case that the canonical method is also parasitism-free.

Corollary 4.3 Let C h satisfy the parasitism-free condition (3.3). Then, composition (4.4) is also parasitism-free.
Proof Repeatedly applying the tableau composition formula of (2.5), we find that the corresponding B, U matrices of (4.4) are respectively given by Inserting these into the LHS of (3.3), we find where we have applied w T BU u = 0 since C h is parasitism-free. Thus, composition (4.4) is also parasitism-free.

Composition of non-canonical methods
Consider now a composition method based on an invertible GLM with arbitrary inputs. The corresponding composition formulae and results all extend straightforwardly from those given in the previous section after making the substitution In particular, the general form of (4.1) is written as where R A (a, b) := T ah • T * bh , and for (4.4) this is bh . In addition, the starting and finishing methods are given by where S h and F h are the starting and finishing methods of the base GLM M h .
Example 4.1 From composition (4.7), we can obtain the GLM version of the triple jump: where, for M h of even order p, α 1 and α 2 are given in (1.3). Similarly, we can obtain the GLM version of the Suzuki 5-jump: where α 1 and α 2 are given in (1.4).

Stage reductions:
Consider the straightforward implementation of a COSY-GLM (e.g. in its canonical form), possibly using multiple iterations, such that the resulting method attains order p D ∈ N. Then, the total number of stages that require evaluation are given by · (s + 2s) when only using (4.8), · (s + 2s) when only using (4.9), wheres denotes the number of stages in the starting method. Improvements to the implementation can be made by identifying and removing any redundant stages prior to integration. For example, the nonlinear map R B (a, b) which is performed between method evaluations has the GLM tableau: Suppose A S is an s × s matrix, then this tableau suggests that a total of 2 s stage equations must be solved for each R B (a, b)-evaluation. However, if we choose U F = 1 S w H 1 , where w 1 is the left eigenvector of V corresponding to eigenvalue ζ = 1, then in the case a = b (see Suzuki composition (4.9)) we find a reduction to s-many stages occurs, i.e. the tableau for R B (a, a) actually reads As reductions of this type are both method and composition dependent, we suggest that each (distinct) nonlinear map R B (a, b) is implemented as an individual GLM, with redundant stages removed. Then, compositions such as the  would be performed in the fashion α 2 ) are each distinct, and irreducible GLMs.

Numerical experiments
In following set of experiments, we consider compositions of methods under the triple jump (4.8) and the . To indicate which specific composition is used we adopt the naming convention T.method for the triple jump and S.method for Suzuki 5-jump. For example, a triple jump of a method named GLM4A would be referred to as T.GLM4A.
The numerical experiments we consider are as follows: Firstly, we computationally verify that the proposed COSY-formulae yield an appropriate order increase. Here, we consider compositions of 4th to 6th-order methods (see also [17,Ch. 7] for higher order). Secondly, we perform efficiency comparisons in terms of accuracy versus function evaluations. Here, accuracy is measured by either the trajectory error or the maximum absolute deviation in the Hamiltonian. Finally, we investigate long-time Hamiltonian preservation of the composition methods on several reversible Hamiltonian systems. In particular, we consider (P2) Bead on a wire (see [1]): For y = [p, q] T , ( p 10 , p 20 , q 10 , q 20 ) = 0, 1 + e 1 − e , 1 − e, 0 , e = 0.6.

Symmetric GLMs
The following GLMs are both 4th-order, (L , P)-symmetric and satsify the parasitismfree condition ( In order to attain 4th-order, both GLMs must approximate the starting vector , which can be achieved using the following (L , P S )-symmetric starting method: where the choice U F = 1 S e T 1 has been made such that T −1 h is also explicit.

Symmetric RKMs
For comparison, we have chosen methods from the class of symmetric DIRKs which are closest, in the sense of structure, to GLM4A and GLM4B, i.e. their stage equations can be solved sequentially.
Remark 5.1 Higher-order Gauss and Lobatto methods are not considered as their stage equations are typically solved 'all-at-once', i.e. in the product space X s , in conjunction with a sophisticated iteration scheme (see e.g. [9]).
The Butcher tableaux for two 4th-order symmetric DIRKs are given as follows: Note that both DIRKs are formed by compositions of the implicit midpoint rule. DIRK43 is based on the triple jump with α 1 , α 2 given by (1.3) and DIRK45 is a Suzuki 5-jump with α 1 , α 2 given by (1.4).

Implementation
We have ensured the update procedure is consistent across all methods. Specifically, the stage equations are solved sequentially using fixed point iteration. The termination criteria are given by where Δ k denotes the difference between iterates k and k + 1. Also, redundant stages in the COSY-GLMs have been identified and removed prior to each integration.

Order confirmation
Consider problem (P3) which is known to be 2π -periodic (see e.g. [ The order of a method is then approximately given by the gradient of the corresponding plot of log(ε h ) − log(h). In Fig. 1, we demonstrate that compositions of GLM4A and GLM4B by both the triple jump and Suzuki 5-jump can be applied to yield methods of order p = 6.

Efficiency comparisons
Next, we investigate the computational efficiency of the COSY-GLMs. In particular, we repeat the experiment of Sect. 5.2 and record the total number of function evaluations made during the integration. Comparing this quantity to the ε h error provides us with an estimate on computational cost versus accuracy. The results from this experiment are shown in Fig. 2. For a fixed accuracy, it can be seen on average that COSY-GLMs require between 1.66 and 2.5 times fewer function evaluations than the DIRK methods. In addition to ε h error versus function evaluations, we also consider the maximum deviation of the absolute Hamiltonian error, i.e. max n |H (y n ) − H (y 0 )|, n ∈ 1, . . . , T h , versus function evaluations. These results are given in Fig. 3. Similar to before, for a fixed accuracy it can been seen that COSY-GLMs require approximately half the number of function evaluations when compared against the DIRK methods.

Long-time hamiltonian preservation
Lastly, we investigate the long-time preservation of the Hamiltonian for problems (P1)-(P3) using a COSY-GLM. As the results on computational efficiency indicate that compositions of GLM4A and GLM4B perform similarly, we restrict our attention to compositions of only GLM4B as this requires fewer function evaluations for a fixed time-step. The results given in Figs. 4, 5 and 6 show that COSY-GLMs approximately preserve the Hamiltonian over long-times. Also, Fig. 6 shows that quadratic invariants such as angular momentum, i.e. L : R 4 → R, are also approximately preserved. Furthermore, we observe that parasitic instability has yet to manifest itself despite the use of coarse time-steps.

Conclusion
A composition technique for generating composite symmetric general linear methods (COSY-GLMs) of arbitrarily high order has been developed and then applied to create new methods of order 6. The process involves using a canonical transformation that alters the starting and finishing methods of a GLM to be in terms of only the preconsistency vectors u and w. The new methods have been shown to be suitable for the long-time integration of reversible Hamiltonian systems that are either separable or non-separable. In particular, numerical experiments have been performed which show that the methods approximately preserve the Hamiltonian and other invariants for the duration of the computation. Furthermore, these methods have been found to be computationally more efficient than symmetric diagonally-implicit Runge-Kutta methods (DIRKs) of the same order, typically requiring half the number of total function evaluations for a fixed accuracy.