Novel high-order energy-preserving diagonally implicit Runge–Kutta schemes for nonlinear Hamiltonian ODEs

Applied Mathematics Letters 102 (2020) 106091

Contents lists available at ScienceDirect

Applied Mathematics Letters www.elsevier.com/locate/aml

Novel high-order energy-preserving diagonally implicit Runge–Kutta schemes for nonlinear Hamiltonian ODEs Hong Zhang, Xu Qian ∗, Songhe Song Department of Mathematics, College of Liberal Arts and Science, National University of Defense Technology, Changsha, 410073, PR China

article

info

Article history: Received 25 September 2019 Received in revised form 9 October 2019 Accepted 9 October 2019 Available online 15 October 2019 Keywords: Nonlinear Hamiltonian ODEs Energy-preserving property Invariant energy quadratization High-order diagonally implicit Runge–Kutta schemes

abstract We utilize the invariant energy quadratization approach to transform the nonlinear Hamiltonian ODE into an equivalent form which has a quadratic energy. Then the reformulation is discretized using a class of diagonally implicit Runge–Kutta schemes. The proposed schemes are proved to conserve a modified quadratic energy conservation law and converge with high-order accuracy. Numerical experiments for several ODEs are provided to illustrate the accuracy, convergence and conservative properties of the proposed method. Comparisons with the symplectic Runge– Kutta scheme and energy-preserving average vector field method are presented to demonstrate the energy-preserving property and accuracy of the proposed method. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction In this letter, we are concerned with numerical integrations of nonlinear Hamiltonian ordinary differential equations. For the past several decades, there have been a great deal of theoretical and numerical studies on the structure preserving methods for this kind of problems [1–3]. The increased interests in this subject could mainly be attributed to the superior qualitative behavior over long time integration of such methods. Among the proposed theories, the energy preserving methods have attracted much attention in recent years. In order to preserve the evolution of the Hamiltonian energy of nonlinear differential equations, a variety of numerical methods have been developed in literature. These include the discrete gradient method [4], the projection method [5], the average vector field (AVF) method [6], the partitioned AVF method [7], the Hamiltonian boundary value method (HBVM) [8], the Kahn method [9], etc. Among the proposed theories, two energy quadratization approaches have been developed widely in recent years. One is the invariant energy quadratization approach which was first proposed by Yang et al. [10,11] to preserve the evolution ˆ for the gradient flows. The essential idea behind the IEQ method is to transform the of a modified energy H energy into a quadratic form of a new variable via a change of variables. This new, equivalent system still ∗ Corresponding author. E-mail addresses: [email protected] (H. Zhang), [email protected] (X. Qian), [email protected] (S. Song).

https://doi.org/10.1016/j.aml.2019.106091 0893-9659/© 2019 Elsevier Ltd. All rights reserved.

2

H. Zhang, X. Qian and S. Song / Applied Mathematics Letters 102 (2020) 106091

retains a similar energy law in terms of new variables. In discretization, the advantage of such a reformulation is that all nonlinear terms could be treated semi-explicitly, which in turn leads to a linear system. This IEQ approach has also been applied in many other problems since then. For example, the molecular beam epitaxial growth model [12], the moving contact line model [13], the Sine–Gordon equation [14], the Lorentz force system [15] and so on. Recently, with the aid of the IEQ technique, Jiang et al. [16] developed a novel class of arbitrarily high-order energy-preserving Runge–Kutta/Gauss collocation methods for the Camassa–Holm equation. The other quadratization technique is called the scalar auxiliary variable (SAV) approach which was proposed by Shen et al. [17] as a modification of the IEQ and overcame several shortcomings of it. For more details and applications of SAV, we refer the interested readers to the Refs. [18–20]. Since the structure preserving methods greatly outperform other methods and the energy quadratization approach could be easily generalized to developing high order energy preserving schemes, the objective of the present work is to develop energy preserving schemes for the nonlinear Hamiltonian ODEs. This method is based on the IEQ approach which works for both conservative and dissipative systems, but we only focus on the conservative system in this work. Different from the approaches in Refs. [16,20], we will adopt a specific class of diagonally implicit Runge–Kutta schemes to conserve the reformulated quadratic energy. We take the classical mechanics as an example. In its Hamiltonian formulation, the generalized coordinates q1 , q2 , . . . , qn and generalized momenta p1 , p2 , . . . , pn are used to represent the stage of a mechanical system. The system of motion is defined in terms of a Hamiltonian function H(p1 , p2 , . . . , pn , q1 , q2 , . . . , qn ) by the equations ⎧ ∂H dp ⎪ ⎪ i =− , ⎨ dt ∂qi i = 1, 2, . . . , n, (1) ∂H dqi ⎪ ⎪ ⎩ = , dt ∂pi where n ≥ 1 is the number of degrees of freedom. Denote pi = yi , qi = yi+n , y = (y1 , y2 , . . . , y2n )T and ∂H T ∂H ∂H , , . . . , ∂y ) . Then the system (1) is usually written as ∇y H = ( ∂y 1 ∂y2 2n y ′ = J −1 ∇y H(y) = f (y),

(2)

with ( J=

0 −In

) In , J−1 = −J = JT , 0

(3)

where In ∈ Rn×n is an identity matrix, and 0 ∈ Rn×n is a zero matrix. One intrinsic property of the Hamiltonian system (2) is the energy conservation law dy dH = ∇y H(y) = ∇y HJ−1 ∇y H = 0. dt dt

(4)

A numerical method is said to be energy-preserving if H(y) is constant along the numerical solution. In this work we will incorporate the well-known IEQ technique together with the quadratic invariant conserving Runge–Kutta scheme to preserve the Hamiltonian energy of high degree. The other parts of the letter are organized as follows. Section 2 introduces the IEQ approach, by which we transform the nonlinear Hamiltonian ODE into a reformulation. Section 3 is devoted to the discretization and analysis of the transformed ODE, there we will present the energy-preserving diagonally implicit Runge– Kutta schemes up to fifth order. In Section 4, several numerical experiments are carried out to demonstrate the effectiveness of the proposed schemes. Finally, Section 5 ends with conclusion and further comments.

H. Zhang, X. Qian and S. Song / Applied Mathematics Letters 102 (2020) 106091

3

2. Reformulation of the nonlinear Hamiltonian ODE using the IEQ approach In this section, we recall that the main challenge to develop efficient, energy preserving schemes for solving the nonlinear Hamiltonian ODEs lies in the discretization of the nonlinear terms. As we know, all RK methods preserve arbitrary linear invariants [21], and the symplectic RK methods preserve arbitrary quadratic invariants [22]. However, no RK methods can preserve arbitrary polynomial invariants of degree higher than two [23]. In order to utilize the quadratic invariant-preserving property of the symplectic RK methods, we will transform the Hamiltonian function of high degree to degree of two by adopting the newly developed IEQ technique. We illustrate the IEQ approach by applying it to a venerable example, the Kepler two-body problem. This is a Hamiltonian system with two degrees of freedom p˙1 = −

q1 (q12

+

3 q22 ) 2

, p˙2 = −

q2 (q12

3

+ q22 ) 2

, q˙1 = p1 , q˙2 = p2 ,

(5)

where p = (p1 , p2 )T and q = (q1 , q2 )T are the velocity and position vectors, respectively. And the Hamiltonian energy reads 1 1 H(p, q) = (p21 + p22 ) − (6) 1 . 2 2 (q + q 2 ) 2 1

2

In the IEQ approach, in order to quadratize the Hamiltonian function, we introduce a new variable r(q) = (q2 +q12 )1/4 and denote 1

2

R(q) = ∇q r(q) = [R1 , R2 ]T ,

where R1 =

q1 ∂r =− 2 , ∂q1 2(q1 + q22 )5/4

Then the system (5) could be transformed to ⎧ ⎪ ⎪ p˙1 ⎪ ⎪ ⎪ ⎪ p˙2 ⎪ ⎨ q˙1 ⎪ ⎪ ⎪ ⎪ q˙2 ⎪ ⎪ ⎪ ⎩ r˙

R2 =

(7)

q2 ∂r =− 2 . ∂q2 2(q1 + q22 )5/4

(8)

= 2R1 r, = 2R2 r, = p1 ,

(9)

= p2 , = R1 p1 + R2 p2 ,

with the initial condition of r given by r(t = 0) =

1 . (q1 (t = 0)2 + q2 (t = 0)2 )1/4

For convenience, we rewrite the above system (9) as { yt = k(y, r), rt = l(y, r).

(10)

(11)

Remark 1. For the nonlinear Hamiltonian ODEs, the SAV approach will provide the same reformulations as the IEQ approach. Proposition 2.1. The transformed system (9) admits the energy conservation law in terms of y and r, 1 ˆ H(y, r) = (p21 + p22 ) − r2 = [y; r]T S[y; r] ≡ Const, 2

(12)

4

H. Zhang, X. Qian and S. Song / Applied Mathematics Letters 102 (2020) 106091

where ⎡1 2

0

⎢0 ⎢ S=⎢ ⎢0 ⎣0 0

0 0 0

1 2

⎤ 0 0 0 0 0 0⎥ ⎥ 0 0 0⎥ ⎥. 0 0 0⎦ 0 0 −1

(13)

Proof . ˆ dH = p1 p1t + p2 p2t − 2rrt dt = 2p1 R1 r + 2p2 R2 r − 2r(R1 p1 + R2 p2 ) = 0. This completes the proof.

□

Remark 2. If Ω is the domain of definition of k(y, r) and l(y, r), differentiation affirms that the conservation law (12) is equal to [k(y, r); l(y, r)]T S[y; r] = 0.

(14)

3. Energy-preserving diagonally implicit Runge–Kutta scheme In this section, we describe the energy-preserving diagonally implicit Runge–Kutta scheme for the reformulated system using IEQ approach. ∑s n n Let bi , aij (i, j = 1, . . . , s) be real numbers and ci = j=1 aij . For given y , r , when an s-stage RK scheme is applied to (9), the following intermediate values are first calculated by ⎧ s ∑ ⎪ n ⎪ y = y + ∆t aij kj , ⎪ i ⎪ ⎨ j=1

s ∑ ⎪ ⎪ n ⎪ ⎪ r = r + ∆t aij lj , i ⎩

(15)

j=1

where kj = k(yj , rj ), and lj = l(yj , rj ). Then (y n+1 , rn+1 ) is updated by ⎧ s ∑ ⎪ n+1 n ⎪ ⎪ y = y + ∆t bi k i , ⎪ ⎨ i=1

s ∑ ⎪ ⎪ n+1 n ⎪ ⎪ r = r + ∆t bi l i . ⎩

(16)

i=1

Theorem 3.1. Let M = diag(b)A + AT diag(b) − bbT , the IEQ-DIRK scheme satisfying the condition M = 0 conserves the energy of the system (9), i.e., ˆ n+1 , rn+1 ) = H(u ˆ n , rn ). H(u

(17)

H. Zhang, X. Qian and S. Song / Applied Mathematics Letters 102 (2020) 106091

5

Table 1 The values of bi for different stages of DIRK satisfying M = 0 (a = 1.351207). Stage (order)

b1 b4

b2 b5

b3 b6

1 (2) [2]

1

2 (2) [2]

1 2

1 2

3 (3) [2]

a

a

1 − 2a

4 (4) [2]

2.70309412 1.8606818856

−0.53652708

2.37893931

5 (4) [28]

−2.150611289942181 1.452223059167718

1.452223059167718 −2.150611289942181

2.3967764615489258

6 (5) [28]

0.5080048194000274 0.5685658926458251

1.360107162294827 −1.4598520495864393

2.0192933591817224 −1.9961191839359627

Proof . By virtue of (14), we have ˆ n+1 , rn+1 ) = [y n+1 ; rn+1 ]T S[y n+1 ; rn+1 ] H(y s ∑ bi [k(yi , ri ); l(yi , ri )]T S[yi ; ri ] = [y n ; rn ]T S[y n ; rn ] + 2∆t    i=1

+ ∆t2

s ∑

=0

(bi bj − bi aij − bj aji )[k(yi , ri ); l(yi , ri )]T S[k(yj , rj ); l(yj , rj )]

(18)

i,j=1

= [y n ; rn ]T S[y n ; rn ] − ∆t2

s ∑

mij [k(yi , ri ); l(yi , ri )]T S[k(yj , rj ); l(yj , rj )].

i,j=1

Since M = 0, then the IEQ-DIRK preserves the quadratic Hamiltonian, i.e. ˆ n+1 , rn+1 ) = H(y ˆ n , rn ). □ H(y

(19)

Lemma 3.2 ([24]). The diagonally implicit Runge–Kutta (DIRK) schemes satisfying the condition M = 0 can be expressed by the Butcher table,

c

A = b

where ci =

∑s

j=1

T

b1 2

c1 c2 .. .

b1 .. .

cs

b2 2

b1

.. . b2

. ···

bs 2

b1

b2

···

bs

..

(20)

aij , and bi ̸= 0 (i = 1, 2, . . . , s), aij = 0 (i < j).

The quadratic invariant-preserving Runge–Kutta schemes have been further discussed by Sanz-Serna and Verwer [25], Franco and G´ omez [26] and so on. Recently some high order diagonally implicit symplectic RK schemes have been proposed by Cong [27] and Kalogiratou [28]. In Table 1 we present the coefficients of several DIRK schemes satisfying M = 0 as will be tested in our work. 4. Numerical experiments In this section, we present some numerical results computed using the IEQ-DIRK methods described in the previous section. First, we demonstrate the numerical convergence orders of the schemes. Then we show

6

H. Zhang, X. Qian and S. Song / Applied Mathematics Letters 102 (2020) 106091

Table 2 Time accuracy tests at T = 10 for different stages of the IEQ-DIRK schemes. Stage (order)

∆t

Error

Order

Stage (order)

∆t

Error

Order

1 (2)

0.1 0.05 0.025 0.0125

3.999e−2 9.854e−3 2.454e−3 6.130e−3

– 2.029 2.007 2.002

2 (2)

0.1 0.05 0.025 0.0125

9.854e−3 2.454e−3 6.130e−4 1.532e−4

– 2.007 2.002 2.000

3 (3)

0.1 0.05 0.025 0.0125

1.067e−3 1.594e−4 2.146e−5 2.778e−6

– 2.679 2.875 2.942

4 (4)

0.1 0.05 0.025 0.0125

3.431e−3 2.340e−4 1.486e−5 9.584e−7

– 3.872 3.877 3.953

5 (4)

0.1 0.05 0.025 0.0125

5.731e−4 3.111e−5 1.990e−6 1.251e−7

– 4.084 4.236 4.176

6 (5)

0.1 0.05 0.025 0.0125

2.101e−5 7.548e−7 2.376e−8 7.784e−10

– 4.784 4.993 5.012

Fig. 1. Conservation of the quadratic Hamiltonian energy (12) from t = 0 to t = 100 for different IEQ-DIRK schemes with ∆t = 0.01.

the accuracy and conservation properties of the proposed schemes. The test cases are very typical and have been enormously studied by researchers in [29,30]. Example 1. First, we test the convergence order and conservation property of the proposed schemes. Consider the Kepler two-body problem (5), when the initial conditions are taken as p1 (0) = 0, p2 (0) = 1, q1 (0) = 1, q2 (0) = 0,

(21)

the exact solution can be given by p1 (t) = − sin(t), p2 (t) = cos(t), q1 (t) = cos(t), q2 (t) = sin(t).

(22)

Choosing the time step size ∆t = 0.1, 0.05, 0.025, 0.0125, the corresponding errors and orders for the IEQDIRK scheme (15), (16) with stages k = 1, 2, 3, 4, 5, 6 are presented in Table 2, from which the desired rates of convergence are observed. In Fig. 1, we display the evolution of the solutions obtained by the IEQ-DIRK schemes with different stages. The results clearly show that the orbits are preserved accurately and the IEQ-DIRK schemes conserve the modified energy (12) to roundoff errors. Example 2. In the second example, we consider the H´enon–Heiles system with Hamiltonian function H(p1 , p2 , q1 , q2 ) =

1 2 1 (p1 + p22 + q12 + q22 ) + q12 q2 − q23 . 2 3

(23)

H. Zhang, X. Qian and S. Song / Applied Mathematics Letters 102 (2020) 106091

7

The equations of motion read ⎧ p˙1 ⎪ ⎪ ⎪ ⎪ ⎨ p˙2 ⎪ ⎪ q˙1 ⎪ ⎪ ⎩ q˙2 This system has a critical energy value Hc = bounded to unbounded orbits.

1 6

= −q1 − 2q1 q2 , = −q2 − q12 + q22 , = p1 ,

(24)

= p2 . at which the qualitative nature of the solution changes from

√ In the H´enon–Heiles system, by introducing a new variable r = q12 q2 − 13 q23 + c0 (c0 > 23 ), we can get an equivalent reformulation ⎧ p˙1 = −q1 − 2rR1 , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ p˙2 = −q2 − 2rR2 , ⎪ ⎨ q˙1 = p1 , (25) ⎪ ⎪ ⎪ ⎪ q˙2 = p2 , ⎪ ⎪ ⎪ ⎩ r˙ = R1 p1 + R2 p2 , where ∂r ∂r q1 q2 q12 − q22 R1 = , R2 = . (26) =√ = √ ∂q1 ∂q2 q12 q2 − 13 q23 + c0 2 q12 q2 − 13 q23 + c0 ˆ = 1 (p2 + p2 + q 2 + q 2 ) + r2 − c0 , the constant c0 > 2 Such that the Hamiltonian is transformed to H 2 1 2 2 1 3 1 3 2 ensures q1 q2 − 3 q2 + c0 > 0 and could be obtained by finding the minimum value of q12 q2 − 31 q23 in the region [−1, 1] × [−1, 1]. Remark 3. For general nonlinear Hamiltonian ODEs, when the nonlinear term of the energy is unbounded from below, a modified stabilized IEQ approach which uses the convex splitting idea could be utilized to obtain a quadratic reformulation. For more detail, we refer to [31,32]. In the numerical experiments, we choose the initial condition p1 = p2 = 0, q1 = 0.1, q2 = −0.5 at a point located on the boundary of the critical triangular region (see Fig. 2) and c0 = 1. For comparisons we solve this system using the second order IEQ-DIRK(1, 2) scheme, the symplectic RK(1, 2) (SRK) scheme and the second order AVF scheme. All the three schemes are implicit. In Fig. 2 we compare the results of the three schemes over the time interval t ∈ [0, 1000]. The demonstrated results show that all three orbits are constrained in the critical triangular region. Fig. 3 shows that the conservation of the modified energy obtained by the IEQ-DIRK scheme is comparative to the results obtained by the AVF scheme and is much better than the symplectic RK(1, 2) scheme. 5. Conclusion Throughout this letter, we studied a class of energy-preserving diagonally implicit Runge–Kutta schemes for the nonlinear Hamiltonian ODEs. First, we transform the original ODE to an equivalent formulation where the Hamiltonian energy is transformed to a quadratic invariant. Based on this reformulation, a special class of diagonally implicit Runge–Kutta schemes satisfying M = 0 are utilized to obtain the energypreserving integrators. By comparing with the classical symplectic RK scheme and the AVF scheme, the numerical results support the theoretical analysis. Future work would extend the invariant energy quadratization technique to the nonlinear Hamiltonian PDEs and investigate efficient explicit energy-preserving Runge–Kutta schemes based on the IEQ reformulation. We also mention that the scalar auxiliary variable (SAV) approach could also be utilized in those cases.

H. Zhang, X. Qian and S. Song / Applied Mathematics Letters 102 (2020) 106091

8

Fig. 2. Hénon–Heiles system: configuration space orbits obtained using the IEQ-DIRK(1, 2), SRK(1, 2) and the AVF integrators. Left: IEQ-DIRK; middle: SRK; right: AVF. The initial state marked on the lower edge of the triangle at q1 = 0.1, q2 = −0.5, p1 = p2 = 0 corresponds to the critical energy H = 1/6. Parameters: ∆t = 0.1, t ∈ [0, 1000].

Fig. 3. Hénon–Heiles system: energy conservation tests from t = 0 to t = 1000 for IEQ-DIRK(1, 2), SRK(1, 2) and AVF integrators.

Acknowledgments This work was supported by the National Natural Science Foundation of China (Nos. 11571366, 11901577, 11971481), the National Natural Science Foundation, PR China of Hunan (No. S2017JJQNJJ0764), the fund from Hunan Provincial Key Laboratory, PR China of Mathematical Modeling and Analysis in Engineering (No. 2018MMAEZD004), the funds from NNW, PR China (No. NNW2018-ZT4A08) and the Research Fund of NUDT, PR China (No. ZK17-03-27). References [1] E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, Vol. 31, Springer, 2006. [2] K. Feng, M. Qin, Symplectic difference schemes for hamiltonian systems, in: Symplectic Geometric Algorithms for Hamiltonian Systems, Springer, 2010, pp. 187–211. [3] D. Furihata, T. Matsuo, Discrete Variational Derivative Method: A Structure-Preserving Numerical Method for Partial Differential Equations, Chapman and Hall/CRC, 2010. [4] B.O. Morten Dahlby, A general framework for deriving integral preserving numerical methods for PDEs, SIAM J. Sci. Comput. 33 (5) (2011) 2318–2340. [5] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II, Springer Berlin Heidelberg, 1996. [6] E. Celledoni, V. Grimm, R.I. McLachlan, D. McLaren, D. ONeale, B. Owren, G. Quispel, Preserving energy resp. dissipation in numerical pdes using the average vector field method, J. Comput. Phys. 231 (20) (2012) 6770–6789.

H. Zhang, X. Qian and S. Song / Applied Mathematics Letters 102 (2020) 106091

9

[7] W. Cai, H. Li, Y. Wang, Partitioned averaged vector field methods, J. Comput. Phys. 370 (2018) 25–42. [8] L. Brugnano, F. Iavernaro, D. Trigiante, Hamiltonian boundary value methods (energy preserving discrete line integral methods), J. Numer. Anal. Ind. Appl. Math. 5 (1–2) (2010) 17–37. [9] E. Celledoni, R.I. McLachlan, B. Owren, G.R.W. Quispel, Geometric properties of kahan’s method, J. Phys. A 46 (2) (2012) 025201. [10] X. Yang, L. Ju, Linear and unconditionally energy stable schemes for the binary fluid–surfactant phase field model, Comput. Methods Appl. Mech. Engrg. 318 (2017) 1005–1029. [11] X. Yang, L. Ju, Efficient linear schemes with unconditional energy stability for the phase field elastic bending energy model, Comput. Methods Appl. Mech. Engrg. 315 (2017) 691–712. [12] X. Yang, J. Zhao, Q. Wang, Numerical approximations for the molecular beam epitaxial growth model based on the invariant energy quadratization method, J. Comput. Phys. 333 (2017) 104–127. [13] X. Yang, H. Yu, Efficient second order unconditionally stable schemes for a phase field moving contact line model using an invariant energy quadratization approach, SIAM J. Sci. Comput. 40 (3) (2018) B889–B914. [14] C. Jiang, W. Cai, Y. Wang, A linearly implicit and local energy-preserving scheme for the Sine-Gordon equation based on the invariant energy quadratization approach, J. Sci. Comput. 80 (3) (2019) 1629–1655. [15] H. Li, Q. Hong, An efficient energy-preserving algorithm for the lorentz force system, Appl. Math. Comput. 358 (2019) 161–168. [16] C. Jiang, Y. Wang, Y. Gong, Arbitrarily high-order energy-preserving schemes for the Camassa–Holm equation, arXiv preprint arXiv:1906.07342. [17] J. Shen, J. Xu, J. Yang, The scalar auxiliary variable (SAV) approach for gradient flows, J. Comput. Phys. 353 (2018) 407–416. [18] J. Shen, J. Xu, Convergence and error analysis for the scalar auxiliary variable (SAV) schemes to gradient flows, SIAM J. Numer. Anal. 56 (5) (2018) 2895–2912. [19] W. Cai, C. Jiang, Y. Wang, Y. Song, Structure-preserving algorithms for the two-dimensional Sine-Gordon equation with Neumann boundary conditions, J. Comput. Phys. 395 (2019) 166–185. [20] Y. Gong, J. Zhao, Q. Wang, Arbitrarily high-order unconditionally energy stable schemes for gradient flow models using the scalar auxiliary variable approach, arXiv preprint arXiv:1907.04254. [21] D.F. Griffiths, D.J. Higham, Numerical Methods for Ordinary Differential Equations: Initial Value Problems, Springer Science & Business Media, 2010. [22] G. Cooper, Stability of runge–kutta methods for trajectory problems, IMA J. Numer. Anal. 7 (1) (1987) 1–13. [23] M. Calvo, A. Iserles, A. Zanna, Numerical solution of isospectral flows, Math. Comput. Amer. Math. Soc. 66 (220) (1997) 1461–1486. [24] M. Qin, M. Zhang, Symplectic runge–kutta schemes for hamiltonian systems, JCM Suppl. Issue (1992) 205–215. [25] J. Sanz-Serna, L. Abia, Order conditions for canonical runge–kutta schemes, SIAM J. Numer. Anal. 28 (4) (1991) 1081–1096. [26] J. Franco, I. G´ omez, Fourth-order symmetric DIRK methods for periodic stiff problems, Numer. Algorithms 32 (2–4) (2003) 317–336. [27] C. Jiang, Y. Cong, A sixth order diagonally implicit symmetric and symplectic runge–kutta method for solving hamiltonian systems, J. Appl. Anal. Comput. 5 (1) (2015) 159–167. [28] Z. Kalogiratou, T. Monovasilis, T. Simos, Diagonally implicit symplectic runge–kutta methods with special properties, in: AIP Conference Proceedings, Vol. 1479, AIP, 2012, pp. 1387–1390. [29] L. Brugnano, F. Iavernaro, D. Trigiante, A. two step, Fourth-order method with energy preserving properties, Comput. Phys. Comm. 183 (9) (2012) 1860–1868. [30] H. Li, Y. Wang, M. Qin, A sixth order averaged vector field method., J. Comput. Math. 34 (5) (2016) 479–498. [31] Z. Liu, Efficient invariant energy quadratization and scalar auxiliary variable approaches without bounded below restriction for phase field models, arXiv preprint arXiv:1906.03621. [32] Z. Liu, X. Li, Efficient modified techniques of invariant energy quadratization approach for gradient flows, Appl. Math. Lett. 98 (2019) 206–214.