Symmetric second derivative integration methods

Accepted Manuscript Symmetric second derivative integration methods M. Hosseini Nasab, A. Abdi, G. Hojjati

PII: DOI: Reference:

S0377-0427(17)30447-8 http://dx.doi.org/10.1016/j.cam.2017.09.016 CAM 11296

To appear in:

Journal of Computational and Applied Mathematics

Received date : 4 May 2017 Revised date : 6 August 2017 Please cite this article as: M. Hosseini Nasab, A. Abdi, G. Hojjati, Symmetric second derivative integration methods, Journal of Computational and Applied Mathematics (2017), http://dx.doi.org/10.1016/j.cam.2017.09.016 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript Click here to view linked References

Symmetric second derivative integration methods M. Hosseini Nasaba , A. Abdia,∗, G. Hojjatia a Faculty

of Mathematical Sciences, University of Tabriz, Tabriz, Iran

Abstract The purpose of this paper is the derivation of symmetric second derivative general linear methods (SGLMs) for general time-reversible differential equations. For this purpose, we find algebraic conditions on the coefficients matrices of the methods which are sufficient for the method to be symmetric. Some symmetric methods of order 4 are constructed and implemented on the famous reversible problems. Numerical experiments show that the constructed symmetric SGLMs approximately conserve the invariants of motion over long time intervals for reversible Hamiltonian systems. Keywords: Hamiltonian systems, General linear methods, Second derivative methods, Symmetric methods, G-symplecticity, Reversible problems

1. Introduction The construction of numerical schemes for solving ordinary differential equations (ODEs) which preserve some qualitative geometrical properties of the flow of the differential equation is a relatively new area of numerical analysis. The symmetric property has been recognized as an important factor for the good long-time behaviour of a numerical method applied to time-reversible differential equations, as it is the case of conservative mechanical systems including Hamiltonian systems. Symmetric methods play a central role in the geometric integration of differential equations which a theoretical explanation of the good long-time behaviour of these methods is given in [20]. The differential equation y′ = f (y) is called ρ -reversible if f (ρ y) = −ρ f (y) for all y,

(1)

with ρ as an invertible linear transformation in the phase space [20]. This implies that the flow of the system, denoted by ϕt , satisfies ρ ◦ ϕt = ϕt−1 ◦ ρ . A Hamiltonian system with the Hamiltonian function satisfying H(p, q) = H(−p, q) is an example of reversible problems with the reflection ρ (p, q) = (−p, q). ∗ Corresponding

author Email addresses: [email protected] (M. Hosseini Nasab), [email protected] (A. Abdi), [email protected] (G. Hojjati)

Preprint submitted to Elsevier

August 6, 2017

To solve a reversible differential equation, we need to search for numerical methods which produce solutions with the same behaviour as the exact ones. A generic one-step integrator yn+1 = Φh (yn ) is symmetric or time-reversible if it satisfies the condition Φh = Φ−1 −h . It is proved in [21, 28] that for all ρ -reversible problems, the numerical flow of a Runge–Kutta method will be also ρ -reversible if and only if it is symmetric. It is shown in [20, 32] that the explicit Runge–Kutta methods cannot be symmetric. However, it is possible to get explicit symmetric integrators with partitioned Runge–Kutta methods for separable problems originating from a system of second order differential equations. Also, within the implicit symmetric Runge–Kutta methods, diagonally implicit Runge–Kutta methods (DIRKs) formed by compositions of the implicit midpoint method are very popular and efficient [23, 27, 29, 33]. The early study of the symmetric multistep methods started by Quinlan and Tremaine [26], where new symmetric methods were successfully tested for simulations of the outer solar system. Multistep methods also suffer from parasitic solutions. However, bounds for parasitic solution components are obtained by Hairer and Lubich in [19, 20]. The properties of symmetric multistep methods were further investigated in [11, 15]. A definition of symmetry for general linear methods (GLMs) is formulated in [20, 7]. A symmetric GLM is defined in a similar fashion to, but not as simple as, a Runge–Kutta method. Sufficient algebraic conditions for a GLM with the coefficients matrices A, U, B, and V to be symmetric are obtained by using the involution matrix L to rearrange the input and output vector such that L2 = I and LV −1 L = V , and a stage reversing permutation matrix P, i.e. P2 = I, (cf. [20, 7]). According to [7], a GLM is symmetric with respect to the involution and permutation matrices L and P if PAP = UV −1 B − A,

UL = PUV −1 ,

BP = V LB,

V L = LV −1 .

Some characterizations of symmetry and its connections with G-symplecticity for the GLMs are given in [7]. Recently, partitioned GLMs with the G-symplecticity property and a special type of symmetry, known as interchange symmetry, for separable Hamiltonian problems were introduced [6]. The advantage of such methods is that for separable Hamiltonian problems the overall numerical scheme is fully explicit. To construct methods with high orders and desirable stability properties, some efficient second derivative methods were introduced [12, 16, 22, 30]. Second derivative general linear methods (SGLMs) as an extension of GLMs to cover second derivative methods were introduced by Butcher and Hojjati in [10] and investigated more in [3, 4, 2, 5, 17, 24]. An additional advantage of SGLMs over GLMs is that these methods have more free parameters which can be used to achieved better properties like higher order and wider region of absolute stability while they usually have lower computational cost, see [5, 1, 24]. G-symplectic SGLMs which preserve a generalization of quadratic invariants along the long-time integration are given in [25]. The spirit of this paper is that of deriving symmetric methods belonging to the family of SGLMs that are well suited for the integration reversible Hamiltonian system. So, we seek conditions of symmetry by making some relations between the coefficients matrices of the methods. The paper is organized as follows. In Section 2, after a short 2

review of SGLMs, a concept of symmetry for SGLMs has defined by using coincidence between the numerical method and its adjoint. Examples of symmetric methods of order p = 4 are constructed in Section 3. These methods have diagonally implicit stage matrices, and one of these methods is also G-symplectic. Numerical experiments are shown in Section 4 and compared with those of an existing symmetric GLM and a symmetric Runge–Kutta method. The numerical results of Section 4 show that the symmetric SGLMs approximately conserve the Hamiltonian of reversible problems over long time intervals. 2. Symmetric SGLMs It is our intention to construct symmetric SGLMs. We begin by reviewing the fundamental theory of SGLMs [3]. 2.1. Introduction to SGLMs SGLMs are multistage and multivalue numerical methods for solving the initial value problem ( ′ y (x) = f (y(x)), x ∈ [x0 , x], (2) y(x0 ) = y0 , where f : Rd → Rd . An SGLM for the numerical solution (2) with fixed stepsize h is defined by  r s s [n−1] [n] [n] [n]  2  a g(Y ) + a f (Y ) + h Y = h  i j i j ∑ ui j y j , i = 1, 2, · · · , s, ∑ ∑ j j  i j=1

j=1

j=1

s s r  [n] [n] [n−1] [n] 2  b f (Y ) + h y = h vi j y j , b g(Y ) +  i j i j ∑ ∑ ∑ j j  i j=1

j=1

j=1

(3)

i = 1, 2, · · · , r,

[n]

where n = 1, 2, · · · , N, h = (x − x0 )/N. Here Yi , i = 1, 2, . . . , s, are approximations [n] to the solution of (2) in the internal points xn−1 + ci h and the vectors f (Yi ) and [n] g(Yi ), i = 1, 2, . . . , s are the stage first and second derivative values, where g(·) = [n] f ′ (·) f (·). The output values yi , i = 1, 2, . . . , r are approximations of order p to the linear combinations of scaled derivatives of y(xn ), i.e., [n]

yi =

p

∑ αik hk y(k) (xn ) + O(h p+1),

k=0

for some real numbers αik , k = 0, 1, . . . , p (cf. [3]). The vector y[n] denotes the computed quantities at step number n, which also represents the incoming values for step number n + 1. These methods are characterized by the abscissa vector c = [c1 c2 . . . cs ]T , the coefficient matrices A = [ai j ] ∈ Rs×s ,

B = [bi j ] ∈ Rr×s ,

A = [ai j ] ∈ Rs×s ,

B = [bi j ] ∈ Rr×s , 3

U = [ui j ] ∈ Rs×r , V = [vi j ] ∈ Rr×r ,

(4)

and the four integers: p, the order of the method; q, the stage order; r, the number of external stages; and s, the number of internal stages. They can be represented by a partitioned (s + r) × (2s + r) matrix of the form   A A U . B B V An SGLM with coefficient matrices (A, A,U, B, B,V ) is: (i) pre-consistent if there exists a vector u ∈ Rr such that Vu = u, Uu = e, where e = [1, 1, . . . , 1]T ∈ Rs .

(ii) consistent if it is pre-consistent and there exists a vector υ ∈ Rr such that Be + V υ = u + υ. (iii) zero-stable if V is power bounded, i.e. supn≥0 kV n k < ∞.

Convergence of an SGLM is guaranteed when we ensure that it is both consistent and zero-stable [3]. 2.2. Adjoint and Symmetric methods In this subsection, we look for some characterizations of symmetric SGLMs. The key for understanding symmetry is the concept of the adjoint method. The adjoint method Φ∗h of a method Φh is the inverse of the original method with reversed timestep −h, i.e. Φ∗h := Φ−1 −h . A numerical integration scheme Φh is symmetric if it is equal to its adjoint Φ∗h , i.e., Φh = Φ∗h . The adjoint method of an SGLM is also an SGLM which undoes the work of original SGLM with reversed time step −h. ∗

∗

Lemma 1. The adjoint method of an SGLM (A, A,U, B, B,V ) is again an SGLM (A∗ , A ,U ∗ , B∗ , B ,V ∗ ) where ∗ A∗ = UV −1 B − A, A = A − UV −1 B, ∗

B = −V −1 B,

B∗ = V −1 B, ∗

U = UV

−1

∗

V =V

,

−1

(5)

.

Proof. The adjoint method is obtained form (3) by replacing h by −h and by exchanging y[n] and y[n−1] . If we apply these changes, we obtain

y

Y [n] = −hA f (Y [n] ) + h2Ag(Y [n] ) + Uy[n],

[n−1]

2

[n]

[n]

[n]

= −hB f (Y ) + h Bg(Y ) + V y .

(6) (7)

Extracting y[n] from the relation (7) and inserting it into the relation (6) gives Y [n] = h(UV −1 B − A) f (Y [n] ) + h2(A − UV −1 B)g(Y [n] ) + UV −1 y[n−1] , y[n] = hV −1 B f (Y [n] ) + h2(−V −1 B)g(Y [n] ) + V −1 y[n−1] ,

which is an SGLM with coefficients is given by (5). 4



⋆

Lemma 2. The SGLM (A, A,U, B, B,V ) is equivalent to the method (A⋆ , A , ⋆ U ⋆ , B⋆ , B ,V ⋆ ) where ⋆

U ⋆ = P−1UL,

⋆

V ⋆ = L−1V L.

A⋆ = P−1 AP,

A = P−1 AP,

B⋆ = L−1 BP,

B = L−1 BP,

(8)

where P is a permutation matrix and L is an invertible matrix. Proof. We consider the change of variables y[n] = Lyˆ[n] ,Y [n] = PYˆ [n] in the method (3). We have f (PY ) = P f (Y ), g(PY ) = Pg(Y ), then the method becomes PYˆ [n] = hAP f (Yˆ [n] ) + h2APg(Yˆ [n] ) + ULyˆ[n−1], Lyˆ[n] = hBP f (Yˆ [n] ) + h2BPg(Yˆ [n] ) + V Lyˆ[n−1] . The method is then defined by the tableau  −1 P AP P−1 AP L−1 BP L−1 BP

P−1UL L−1V L



.

(9)

Since P is permutation matrix, it does not affect the result of numerical scheme. Also, L only changes the coordinate basis. Thus, the tableau (9) is as an equivalent tableau for the method (A, A,U, B, B,V ).  Theorem 3. Let L ∈ Rr×r be an involution matrix such that L2 = I and P ∈ Rs×s a permutation matrix defined as Pi j = δi+ j,s+1 for i, j = 1, . . . , s. An SGLM (3) is symmetric if PAP = UV −1 B − A, PAP = A − UV −1 B, UL = PUV −1 , (10) BP = −V LB, V L = LV −1 . BP = V LB, Proof. By applying Lemmas 1 and 2, the conditions of theorem are sufficient to ensure the SGLM is auto-adjoint, and hence symmetric.  The multivalue nature of SGLMs, as the GLMs, exposes them to a type of instability referred to as parasitism. Parasitism is a phenomenon that describes the unacceptable growth of perturbations made to the components of the numerical solution associated with the non-principal eigenvalues of V . In [9] it was demonstrated that if parasitic growth parameters are zero, acceptable behaviour is obtained by the numerical solution without being distorted by the parasitic solutions over long time intervals. In [14] the authors have derived sharp estimates for the parasitic components and the error in the Hamiltonian numerically computed by a multi-value method using the idea of backward error analysis. They have proved that for carefully constructed methods (symmetric and zero growth parameters) the error in the parasitic components typically grows like h p+4 exp (h2 Lx), where L depends on the problem and on the coefficients of the method. Here, we are going to construct methods to be parasitism-free in which the influence of the parasitic components is largely removed. So, we would expect the methods have acceptable behaviour on the long term integrations what it will be shown in the numerical experiments.

5

We recall from [7] that in the canonical form of the methods based on V -diagonalization, V = V0 ⊕V−1 ⊕V1 ⊕ · · ·⊕V0, where V0 = diag(1, 1, . . . , 1), V−1 = diag(−1, −1, . . . , −1), Vi = diag(ζi , ζi , . . . , ζi , ζi , ζi , . . . , ζi ), i = 1, 2, . . . , k. The number of diagonal elements in these blocks are respectively m0 , m−1 and 2mi . Because of consistency of the method m0 ≥ 1, but it is possible that m−1 = 0, indicating that this block is missing. It is assumed that mi ≥ 1, although it is possible that k = 0 indicating that the final k blocks in V do not exist. In a similar way to [7] for GLMs, to construct symmetric SGLMs in the canonical coordinates to be parasitism-free, we can state the following proposition. Proposition 4. A symmetric SGLM in the canonical coordinates is parasitism-free if the (i, i) component for 2 ≤ i ≤ r of the matrix product BU is zero. 3. Construction of symmetric SGLMs In this section we construct two symmetric methods of order p = 4 and free of parasitism. We first recall the order conditions for a method of order p and stage order q. We assume that components of the input vector y[n−1] for the next step satisfy the conditions [n−1]

yi

p

=

∑ αik hk y(k) (xn−1 ) + O(h p+1),

i = 1, 2, . . . , r,

(11)

k=0

for some real parameters αik , i = 1, 2, . . . , r, k = 0, 1, . . . , p. Then, we require the internal [n] stage vector Y [n] = [Yi ]si=1 to be an approximation of order q to the vector y(xn−1 + ch) = [y(xn−1 + ci h)]si=1 , i.e. [n]

Yi

= y(xn−1 + ci h) + O(hq+1),

i = 1, 2, . . . , s,

and the output vector y[n] satisfy [n]

yi =

p

∑ αik hk y(k) (xn ) + O(h p+1),

i = 1, 2, . . . , r,

(12)

k=0

for the same numbers αik . Considering these assumptions, an SGLM has order p and stage order q with the stability order p if [4]  k Ack−1 c Ack−2    − − − U αk = 0, k = 0, 1, . . . , q,   k! (k − 1)! (k − 2)!    k α (13) Bck−2 Bck−1 k−l  k = 0, 1, . . . , p, − − − V αk = 0, ∑   l=0 l! (k − 1)! (k − 2)!      p(exp(z), z) = O(z p+1 ),

where

αk := [α1k α2k . . . αrk ]T ∈ Rr , 6

k = 0, 1, . . . , p,

and p(w, z) is the stability polynomial of SGLM (3) defined by p(w, z) = det wI −
M(z) in which M(z) = V + (zB + z2 B)(I − zA − z2 A)−1U is the stability matrix of the method. Because of lower cost computing, in the construction of the methods, we assume that the matrices A and A have the lower triangular form. We will consider the methods with r = s = 2 and r = 2, s = 3 and assume that the first component of the external vector approximates the exact solution, i.e. we consider # " y(xn ) [n] . y ≈ [n] y2 (xn ) Under these assumptions and applying preconsistent conditions [3], the structure of the methods would be   λ1 0 · · · 0 µ1 0 · · · 0 1 u12    a21 λ2 · · · 0 a21 µ2 · · · 0 1 u22      .  .. . . .. .. .. ..  .. . . .. A A U  .. . . . . . . . .  . = .  as1 as2 · · · λs as1 as2 · · · µs 1 us2  B B V      b11 b12 · · · b1s b11 b12 · · · b1s 1 v12  b21 b22 · · · b2s

b21 b22 · · · b2s

0 v22

We present two examples of symmetric SGLMs of order 4 with two possible options L = diag(1, 1) and L = diag(1, −1) which one of them is G-symplectic (cf. [25]) too. 3.1. A symmetric method with p = 4, q = 2 We consider two-stage symmetric SGLMs of order p = 4 and stage-order q = 2 with L = diag(1, −1). By applying symmetric conditions (10), we obtain SGLM   −a22 + b1 0 b2 u1 + a22 − b1 0 1 u1      b2 u1 − b1 a22 1 u1  b1 a22 A A U   = .   B B V b b −b b 1 0 1 1 1 1   0 0 −b2 b2 0 −1

In order to achieve order 4, due to symmetry, we only impose conditions of order 3 [7]. Since the method is in the canonical coordinates, by Proposition 4, the parasitism-free conditions reduce to (BU)22 = 0. After applying the order conditions, (BU)22 = 0, and 1 and a22 = 12 , the coefficient matrices of the setting the free parameters as b2 = u1 = 10 method take the following form   1 0 0 0 0 1 10     1  7 7  1 1 1 10 A A U  2 2  75 − 75 = (14) , 1 1  1 1 B B V 1 0  12 − 12  2 2  1 1 0 0 − 10 0 −1 10 7

 T with c = 0, 1 and the external vector 

y(xn )



 + O(h5 ). y[n] =  1 h3 y′′′ (xn ) 20

This method is symmetric with L = diag(1, −1) and # " 0 1 , P= 1 0 and has only one implicit stage.

3.2. A G-symplectic symmetric method with p = 4, q = 1 We derive a method of order 4 with s = 3 and q = 1. In addition to being symmetric, the method will be designed to possess G-symplecticity. Applying the order, G-symplecticity [25], symmetric and parasitism control conditions [25], and setting 1 1 and α24 = 30 , the coefficients matrices of the method the free parameters as b12 = 1000 take the following form   A A U = B B V √   1 3 1 1 − 3200 0 0 0 0 0 1 6000 2 6γ   √ 3 1 1 1 1  −6δ 0 0 0 0 1 6000 − 3200 2  3γ    √  1 1 1 1 1 3  −3δ 0 0 0 1 6000 − 3200 2  ,  3γ 6γ   √ √ 3   3 1 1 1 2 1 1 1 1 √ √ − 1 2 − − −   −√ 3 3 2000 1000 2000 800 1600 2−2 3 2−2 2−2   8 16 8 √ 0 0 0 − √ 0 −1 − √ 3 3 3 5( 2−2)

5( 2−2)

5( 2−2)

(15)

√ √ √ 3 where γ = 2 + 24 + 3 2 and δ = (1 + 3 2)2 . This method is G-symplectic with respect to the matrices

G=

"

1 1 6000

−

√ 1 3 3200 2

√

1 1 3 6000 − 3200 2 √ 3 1 2 10240000 ( 2 − 2)

and symmetric with L = diag(1, 1) and 

0 0

 P=  0 1 1 0 8

#

1

,



1 3γ

 D=  0 0



 0  . 0

0 − 31 δ 0

0



 0  , 1 3γ

Here the external vector takes the form " # y(xn ) [n] y = + O(h5). 1 4 (4) 30 h y (xn ) 3.3. Starting procedures To approximate the initial value vector, y[0] , a starting procedure is required. We define the internal stage values of the starting method, and the first and second derivatives at these internal points as       Ye1 f (Ye1 ) g(Ye1 ) Ye :=  Ye2  , f (Ye ) :=  f (Ye2 )  , g(Ye ) :=  g(Ye2 )  , Ye3 f (Ye3 ) g(Ye3 ) respectively. We assume that the starting procedure as follows

e Ye ) + Uy e f (Ye ) + h2 Ag( e 0, Ye = h A e Ye ) + Vy e . y[0] = h Be f (Ye ) + h2 Bg( 0

The starting procedure can be represented by a partitioned matrix " # e U e A e A . e V e Be B

e and V e need to be chosen as U e= Since we have only y0 available, the matrices U e as the [1, 1, 1]T and Ve = [1, 0]T . Also, we take the abscissa vector c˜ and the matrix A same in the fourth-order Lobatto IIIA method. So, the structure of the starting method is   0 0 0 0 0 0 1  5 1  1    24 3 − 24 0 0 0 1     1 2 1 0 0 0 1   6 3 6 .    0 0 0 0 0 0 1    b˜ 1 b˜ 2 b˜ 3 0 0 0 0 The b˜ i values can be found by solving the following equations 3

∑ b˜ i = 0,

i=1 3

∑ b˜ ic˜i = θ1 ,

i=1 3

1 ∑ b˜ i c˜2i = θ2 , 2 i=1 9

1 where for the method (14), θ1 = 20 , θ2 = 0 and for the method (15), θ1 = 0, θ2 = This implies     3 − 20 b˜ 1      b˜ 2  =  1  ,   5   1 b˜ 3 − 20

for the method (14), and

for the method (15).



b˜ 1





    b˜ 2  =     b˜ 3

2 15 4 − 15 2 15

1 30 .



 , 

4. Numerical experiments This section presents the numerical results of the constructed methods in Section 3 for some selected Hamiltonian problems. Our aim is to prepare and compare the numerical results of the following methods: S-SGLM4: SGLM of order 4 (14) which is symmetric but not G-symplectic; G-SGLM4: SGLM of order 4 (15) which is symmetric and G-symplectic; GLM4: GLM of order 4 with only one implicit stage which is symmetric and free of parasitism but not G-symplectic. This method is given in [7] with the coefficients matrices   0 0 0 0 1 1    0 12 0 0 1 1       1 1   3 6 0 0 1 −1  A U  . = (16)  B V  31 16 21 0 1 −1     1 1 1 1   6 3 3 6 1 0  1 6

− 16

1 6

− 16

0 −1

G-GLM4: GLM of order 4 which is symmetric and G-symplectic. This method is given in [13] with the coefficients matrices  1 1  0 0 1 24 6γ  1 1 1   γ  − δ 0 1 3 6 24      1  A U 1 1  1  −3δ γ 1 24  , =  3γ (17) 6 B V  √ √  3 3  1ϕ −1 − 2 2 − 4 1ϕ 1 1   6 4 3 3 6 12  1

−2

10

1

0 −1

where γ = 2 +

√ 3 4 2

√ √ + 3 2, δ = (1 + 3 2)2 and ϕ =

15 4

√ √ + 2 3 2 + 3 4.

Lobatto IIIB: 3-stage Lobatto IIIB method of order 4 which is symmetric but not symplectic (see[8]). Computational experiments are doing by applying the mentioned methods on the following three Hamiltonian problems: – The Kepler problem [20]: This problem is a special case of the two-body problem, which describes the motion of a single planet moving around a heavy sun. The equations of motion are defined by the Hamiltonian 1 1 . H = (p21 + p22 ) − q 2 q2 + q2 1

When started at

q1 (0) = 1 − e, q2 (0) = 0,

p1 (0) = 0,

2

p2 (0) =

r

1+e , 1−e

it has an elliptic periodic orbit of period 2π and eccentricity e ∈ [0, 1). The system has two invariants namely, the Hamiltonian H and angular momentum L given as L = q1 p2 − q2 p1 . The exact solution of the Kepler problem is given in [31] which is p q1 = cos(E) − e, q2 = 1 − e2 sin(E),

where E is the eccentric anomaly given by the Kepler formula x = E − e sin(E).

This problem has several symmetries and satisfies condition (1) with

ρ (p1 , p2 , q1 , q2 ) = (−p1 , −p2 , q1 , q2 ), (−p1 , p2 , q1 , −q2 ), or (p1 , −p2 , −q1 , q2 ). – The double pendulum [20]: The equations of motion a double pendulum are defined by the non-separable Hamiltonian H=

p21 + 2p22 − 2p1 p2 cos(q1 − q2) − cos(q2 ) − 2 cos(q1 ). 2(1 + sin2 (q1 − q2 ))

The initial conditions are taken to be q1 (0) = 3.14, q2 (0) = −3.1, p1 (0) = 0, p2 (0) = 0. For ρ (p1 , p2 , q1 , q2 ) = (−p1 , −p2 , q1 , q2 ) or (p1 , p2 , −q1 , −q2 ), the system is ρ -reversible.

11

–A non-reversible problem [18]: The equations of motion are defined by the separable Hamiltonian H(p, q) =

p3 p q6 q4 q3 1 − + + − + . 3 2 30 4 3 6

The Hamiltonian is not an even function in p and therefore it does not lead to a reversible system. The initial conditions are taken to be p(0) = 1, q(0) = 0. The solution is a non-symmetric orbit. To compare the accuracy of the methods, and at the same time to verify the order behaviour, we solve the Kepler problem with different stepsizes. Plots of global error versus stepsize for the G-SGLM4, S-SGLM4, G-GLM4, and Lobatto IIIB are given in Figure 1. This figure confirms that the order of convergence of the methods match the theoretical predictions and their accuracy are competitive. As a reasonable measure of the efficiency of the constructed methods of order 4, we compare the number of function evaluations, nfe, versus the accuracy of the methods. In Figures 2 and 3, the plot of nfe versus the maximum error in the Hamiltonian is given for the Kepler and the double pendulum problems using the symmetric methods and G-symplectic symmetric methods, respectively, which clearly show the superiority of the methods S-SGLM4 and G-SGLM4 over the other methods.

G-SGLM4 S-SGLM4 G-GLM4 Lobatto IIIB slope of order 4

log10 (kerrork)

-04

-06

-08

-10

-2.4

-2.1

-1.8

-1.5

-1.2

log10 h Figure 1: Numerical results for the Kepler problem (e = 0.5) over half period.

12

10 -3

10 -3 Lobatto IIIB GLM4 S-SGLM4

Lobatto IIIB GLM4 S-SGLM4

10 -4

error in Hamiltonian

error in Hamiltonian

10 -4

10 -5

10 -6

10 -5

10 -6

10 -7

10 -8 10

-7

10 -9

10 -8

10 5

10 5.5

10 -10 10 5

10 6

10 5.5

10 6

nfe

10 6.5

10 7

nfe

Figure 2: The number of function evaluations versus maximum error in the Hamiltonian for the S-SGLM (14), GLM4 (16) and Lobatto IIIB applied to the Kepler problem with e = 0.5 (left) and the double pendulum problem (right).

10-5

10-3 G-GLM4 G-SGLM4

10

error in Hamiltonian

error in Hamiltonian

10

G-GLM4 G-SGLM4

-6

10-7

10-8

10-9

10-10

10-11

-4

10-5

10-6

10-7

10-8

105.5

106

106.5

10-9

107

nfe

105.5

106

106.5

107

nfe

Figure 3: The number of function evaluations versus maximum error in the Hamiltonian for the G-SGLM (15) and G-GLM4 (17) applied to the Kepler problem with e = 0.5 (left) and the double pendulum problem (right).

To investigate possible drift in the Hamiltonian, for each of the three test problems, Figures 4, 6 and 7 are presented. Also, the deviations in the angular momentum for the Kepler problem are presented in Figure 5.

13

3e-07 2e-07 1e-07 0

10-2

10-1

100

101

102

103

101

102

103

(a) S-SGLM4 4e-07 3e-07 2e-07 1e-07 0 10-2

10-1

100

(b) G-SGLM4 Figure 4: The variation in the Hamiltonian for the Kepler problem (e = 0.6) using S-SGLM4 (14), G-SGLM4 (15) with h = 0.01.

6e-08 4e-08 2e-08 0 10-2

10-1

100

101

102

103

101

102

103

(a) S-SGLM4 2e-09 1e-09 0 -1e-09 -2e-09 10-2

10-1

100

(b) G-SGLM4 Figure 5: The variation in the angular momentum for the Kepler problem (e = 0.6) using S-SGLM4 (14), G-SGLM4 (15) with h = 0.01.

14

1e-05 8e-06 6e-06 4e-06 2e-06 0 10-2

10-1

100

101

102

103

104

102

103

104

(a) S-SGLM4 4e-05 3e-05 2e-05 1e-05 0 10-2

10-1

100

101

(b) G-SGLM4 Figure 6: The variation in the Hamiltonian for the double pendulum problem using S-SGLM4 (14), GSGLM4 (15) with h = 0.01.

2e-08 1e-08 0 -1e-08 10-1

100

101

102

103

104

102

103

104

(a) S-SGLM4 6e-08 4e-08 2e-08 0 10-1

100

101

(b) G-SGLM Figure 7: The variation in the Hamiltonian for the Non-reversible problem using S-SGLM4 (14), G-SGLM4 (15) with h = 0.02.

In Figures 4, 5 and 6, we observe that two methods nearly preserve the invariants of motion along the numerical solution and there is no apparent drift for the Kepler and double pendulum problems. Figure 7 shows the error in the Hamiltonian for Nonreversible problem. For this problem, visible energy drift appears in the numerical results of the S-SGLM4 while, as shown in Figure 7b, it does not appear with the case of 15

G-SGLM4 which has G-symplecticity property. The lack of symmetry in this problem and the lack of G-symplecticity in the S-SGLM4 cause the poor energy conservation. However, the results for the G-SGLM4 show that good conservation is possible for SGLMs. 5. Conclusion This paper is concerned with the actual possibility to employ SGLMs for the numerical treatment of reversible Hamiltonian systems. To do this, the SGLMs have been equipped with symmetry feature. Examples of symmetric SGLM have been provided, tested and compared with known symmetric methods in point of view of computational cost. The results for the symmetric SGLM showed that such methods have good energy conservation when applied to the reversible Hamiltonian systems. For nonsymmetric Hamiltonian system, we observed that the symmetric SGLM is not good enough while the G-symplectic symmetric method with zero parasitic growth factors can give favourable long-time behaviour. The future work will involve to deriving higher order efficient G-symplectic symmetric methods for Hamiltonian systems and also the extensive analysis of the SGLMs to get estimates for the parasitic solution components and error in the Hamiltonian. References [1] A. Abdi, Construction of high-order quadratically stable second-derivative general linear methods for the numerical integration of stiff ODEs, J. Comput. Appl. Math. 303 (2016) 218–228. [2] A. Abdi, M. Bra´s, G. Hojjati, On the construction of second derivative diagonally implicit multistage integration methods for ODEs, Appl. Numer. Math. 76 (2014) 1–18. [3] A. Abdi, G. Hojjati, An extension of general linear methods, Numer. Algor. 57 (2011) 149–167. [4] A. Abdi, G. Hojjati, Maximal order for second derivative general linear methods with Runge–Kutta stability, Appl. Numer. Math. 61 (2011) 1046–1058. [5] A. Abdi, G. Hojjati, Implementation of Nordsieck second derivative methods for stiff ODEs, Appl. Numer. Math. 94 (2015) 241–253. [6] J.C. Butcher, R. D’Ambrosio, Partitioned general linear methods for separable Hamiltonian problems, Appl. Numer. Math. 117 (2017) 69–86. [7] J.C. Butcher, A. Hill, T. Norton, Symmetric general linear methods, BIT Numer. Math. 56 (2016) 1189–1212. [8] J.C. Butcher, Numerical methods for ordinary differential equations, Wiley, New York, 2016.

16

[9] J.C. Butcher, Y. Habib, A. T. Hill, T. J. Norton, The control of parasitism in Gsymplectic methods, SIAM J. Numer. Anal. 52 (2014) 2440–2465. [10] J.C. Butcher, G. Hojjati, Second derivative methods with RK stability, Numer. Algor. 40 (2005) 415–429. [11] B. Cano, J. Sanz-Serna, Error growth in the numerical integration of periodic orbits by multistep methods, with application to reversible systems, IMA J. Numer. Anal. 18 (1998) 57–75. [12] R.P. Chan, A. Y. Tsai, On explicit two-derivative Runge–Kutta methods, Numer. Algor. 53 (2010) 171–194. [13] R. D’Ambrosio, G. De Martino, B. Paternoster, Numerical integration of Hamiltonian problems by G-symplectic methods, Adv. Comput. Math. 40 (2014) 553– 575. [14] R. D’Ambrosio, E. Hairer, Long-term stability of multi-value methods for ordinary differential equations, J. Sci. Comput. 60 (2014) 627–640. [15] T. Eirola, J. Sanz-Serna, Conservation of integrals and symplectic structure in the integration of differential equations by multistep methods, Numer. Math. 61 (1992) 281–290. [16] W. Enright, Second derivative multistep methods for stiff ordinary differential equations, SIAM J. Numer. Anal. 11 (1974) 321–331. [17] A. K. Ezzeddine, G. Hojjati, A. Abdi, Perturbed second derivative multistep methods, J. Numer. Math. 23 (2015) 235–245. [18] E. Faou, E. Hairer, T.-L. Pham, Energy conservation with non-symplectic methods: examples and counter-examples, BIT Numer. Math. 44 (2004) 699–709. [19] E. Hairer, C. Lubich, Symmetric multistep methods over long times, Numer. Math. 97 (2004) 699–723. [20] E. Hairer, C. Lubich, G. Wanner, Geometric numerical integration: structurepreserving algorithms for ordinary differential equations, vol. 31, Springer Science & Business Media, 2006. [21] E. Hairer, D. Stoffer, Reversible long-term integration with variable stepsizes, SIAM J. Sci. Comput. 18 (1997) 257–269. [22] G. Hojjati, M.Y. Rahimi Ardabili, S. Hosseini, New second derivative multistep methods for stiff systems, Appl. Math. Model. 30 (2006) 466–476. [23] R.I. McLachlan, On the numerical integration of ordinary differential equations by symmetric composition methods, SIAM J. Sci. Comput. 16 (1995) 151–168. [24] A. Movahedinejad, G. Hojjati, A. Abdi, Second derivative general linear methods with inherent Runge–Kutta stability, Numer. Algor. 73 (2016) 371–389. 17

[25] M. Hosseini Nasab, G. Hojjati, A. Abdi, G-symplectic second derivative general linear methods for Hamiltonian problems, J. Comput. Appl. Math. 313 (2017) 486–498. [26] G.D. Quinlan, S. Tremaine, Symmetric multistep methods for the numerical integration of planetary orbits, Astron. J. 100 (1990) 1694–1700. [27] J.M. Sanz-Serna, L. Abia, Order conditions for canonical Runge–Kutta schemes, SIAM J. Numer. Anal. 28 (1991) 1081–1096. [28] D. Stoffer, Variable steps for reversible integration methods, Computing 55 (1995) 1–22. [29] M. Suzuki, Fractal decomposition of exponential operators with applications to many-body theories and monte carlo simulations, Phys. Lett. A 146 (1990) 319– 323. [30] A.Y. Tsai, R.P. Chan, S. Wang, Two-derivative Runge–Kutta methods for PDEs using a novel discretization approach, Numer. Algor. 65 (2014) 687–703. [31] H. Van de Vyver, An embedded exponentially fitted Runge–Kutta–Nystr¨om method for the numerical solution of orbital problems, New Astron. 11 (2006) 577–587. [32] G. Wanner, Runge–Kutta methods with expansion in even powers of h, Computing 11 (1973) 81–85. [33] H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett. A 150 (1990) 262–268.

18