Symmetric and symplectic exponential integrators for nonlinear Hamiltonian systems

# Symmetric and symplectic exponential integrators for nonlinear Hamiltonian systems

Yajun Wu 111School of Mathematical Sciences, Qufu Normal University, Qufu 273165, P.R.China. E-mail: 1921170786@qq.com    Bin Wang School of Mathematical Sciences, Qufu Normal University, Qufu 273165, P.R.China; Mathematisches Institut, University of Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany. E-mail: wang@na.uni-tuebingen.de
###### Abstract

This letter studies symmetric and symplectic exponential integrators when applied to numerically computing nonlinear Hamiltonian systems. We first establish the symmetry and symplecticity conditions of exponential integrators and then show that these conditions are extensions of the symmetry and symplecticity conditions of Runge-Kutta methods. Based on these conditions, some symmetric and symplectic exponential integrators up to order four are derived. Two numerical experiments are carried out and the results demonstrate the remarkable numerical behavior of the new exponential integrators in comparison with some symmetric and symplectic Runge-Kutta methods in the literature.

Keywords: exponential integrators; symmetric methods; symplectic methods; Hamiltonian systems

MSC (2000): 65L05, 65P10

## 1 Introduction

In this letter, we explore efficient symmetric and symplectic methods for solving the initial value problems expressed in the following from

 y′(t)=My(t)+f(y(t)),t∈[t0,T],y(t0)=y0, (1)

where is assumed to be a linear operator on a Banach space with a norm , is the infinitesimal generator of a strongly continuous semigroup on and the function is analytic (see, e.g. ). It follows from the assumption of that there exist two constants and such that

 ∥∥etM∥∥X←X≤Ceωt,     t≥0. (2)

We note that the linear operator can be a matrix if is chosen as or . Under this situation, is accordingly the matrix exponential function. It is known that the exact solution of (1) can be represented by the variation-of-constants formula

 y(t)=etMy0+∫t0e(t−τ)Mf(y(τ))dτ. (3)

Problems of the form (1) often arise in a wide range of practical applications such as quantum physics, engineering, flexible body dynamics, mechanics, circuit simulations and other applied sciences (see, e.g. [5, 7, 18, 21]). Some highly oscillatory problems, Schrödinger equations, parabolic partial differential equations with their spatial discretisations all fit the form. In order to solve (1) effectively, many researches have been done and the readers are referred to [9, 10, 11, 12, 13] for example. Among them, a standard form of exponential integrators is formulated and these integrators have been studied by many researchers. We refer the reader to [1, 2, 3, 4, 5, 8, 15, 17, 19, 20] for some examples on this topic and a systematic survey of exponential integrators is referred to .

On the other hand, it can be observed that the problem (1) can become a nonlinear Hamiltonian system if

 f(y)=J−1∇U(y),  M=J−1Q,

where is a smooth potential function, is a symmetric matrix, and with the identity . The energy of this Hamiltonian system is . For this system, symplectic exponential integrators (see ) are strongly recommended since they can preserve the symplecticity of the original problems and provide good long time energy preservation and stability. Besides, it is shown in  that symmetric methods also have excellent long time behaviour when applied to reversible differential equations and symmetric exponential integrators have been considered for solving Schrödinger equations in . However, it seems that symmetric exponential integrators have not been used for ODEs and moreover symmetric and symplectic exponential integrators have never been studied so far, which motives this letter.

The main contribution of this letter is to analyse and derive symmetric and symplectic exponential integrators. The integrators have symmerty and symplecticity simultaneously. The letter is organized as follows. In Section 2, we present the scheme of exponential integrators and derive its properties including the symmetry and symplecticity conditions. Then we are devoted to the construction of some practical symmetric and symplectic exponential integrators in Section 3. In Section 4, we carry out two numerical experiments and the numerical results demonstrate the remarkable efficiency of the new integrators in comparison with some existing methods in the scientific literature. The last section is concerned with conclusions.

## 2 Exponential integrators and their properties

In this section, we first present the scheme of exponential integrators and then analyze their symmetry and symplecticity conditions.

###### Definition 2.1

(See .) An -stage exponential integrator (EI) for the problem (1) is defined by

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩Yi=ecihMy0+hs∑j=1¯aij(hM)f(Yj), i=1,…,s,y1=ehMy0+hs∑i=1¯bi(hM)f(Yi), (4)

where are constants and and are matrix-valued functions of .

###### Remark 2.2

It is worth mentioning that if an exponential integrator reduces to a classical Runge-Kutta (RK) method with the coefficients for .

The next theorem gives the symmetry conditions of the EI method.

###### Theorem 2.3

An s-stage exponential integrator is symmetric if and only if its coefficients satisfy the following conditions:

 ci=1−cs+1−i,i=1,2,…,s,¯aij(hM)=ecihM¯bs+1−j(−hM)−¯as+1−i,s+1−j(−hM),i,j=1,2,…,s,¯bi(hM)=ehM¯bs+1−i(−hM),i=1,2,…,s. (5)

Proof  For the exponential integrator , exchanging and replacing by yields

 ˆYi=e−cihMy1−hs∑j=1¯aij(−hM)f(ˆYj), i=1,…,s,y0=e−hMy1−hs∑i=1¯bi(−hM)f(ˆYi). (6)

Then we have

 y1=ehMy0+hs∑i=1ehM¯bi(−hM)f(ˆYi),ˆYi=e(1−ci)hMy0+hs∑j=1e(1−ci)hM¯bj(−hM)f(ˆYj)−hs∑j=1¯aij(−hM)f(ˆYj)=e(1−cihMy0+hs∑j=1(e(1−ci)hM¯bj(−hM)−¯aij(−hM))f(ˆYj). (7)

For the second formula of (7) and the first formula of , it is required that the following conditions are true

 Y1=ˆYs,  Y2=ˆYs−1,  ⋯,  Ys=ˆY1.

Based on these conditions, we obtain that the following two formulae are equal

 Yi=ecihMy0+hs∑j=1¯aij(hM)f(Yj),ˆYs+1−i=e(1−cs+1−i)hMy0+hs∑j=1(e(1−cs+1−i)hM¯bj(−hM)−¯as+1−i,j(−hM))f(ˆYj). (8)

This implies that

 ci=1−cs+1−i,i=1,2,…,s,¯aij(hM)=e(1−cs+1−i)hM¯bs+1−j(−hM)−¯as+1−i,s+1−j(−hM),i,j=1,2,…,s. (9)

Again, according to the second formula of and the first formula of , we obtain the third result of (5). Therefore, the exponential integrator is symmetric if and only if the conditions (5) hold.

###### Remark 2.4

It is noted that when , these symmetry conditions become

 ci=1−cs+1−i,i=1,2,…,s,¯aij(0)=¯bs+1−j(0)−¯as+1−i,s+1−j(0),i,j=1,2,…,s,¯bi(0)=¯bs+1−i(0),i=1,2,…,s, (10)

which are the exact symmetry conditions of -stage RK methods.

About the symplecticity conditions of exponential integrators, we have the following result.

###### Theorem 2.5

(See .) If the coefficients of an -stage exponential integrator satisfy

 ¯bi(hM)TJSS−1i=S−TiSTJ¯bi(hM)=γJ,γ∈R,  i=1,2,…,s,¯bi(hM)TJ¯bj(hM)=¯bi(hM)TJSS−1i¯aij(hM)+¯aji(hM)TS−TjSTJ¯bj(hM),i,j=1,2,…,s, (11)

where and for then the integrator is symplectic.

###### Remark 2.6

We also remark that when , these conditions reduce to

 ¯bi(0)¯bj(0)=¯bi(0)¯aij(0)+¯bj(0)¯aji(0), i,j=1,2,…,s, (12)

which are the exact symplecticity conditions of -stage RK methods.

## 3 Symmetric and symplectic EI

In this section, we derive a class of symmetric and symplectic exponential integrators. To this end, we consider the following special exponential integrators.

###### Definition 3.1

(See ) Define a special kind of -stage exponential integrators by

 ¯aij(hM)=aije(ci−cj)hM,  ¯bi(hM)=bie(1−ci)hM,  i,j=1,…,s, (13)

where

 c=(c1,…,cs)⊺,  b=(b1,…,bs)⊺,  A=(aij)s×s (14)

are the coefficients of an -stage RK method. We denote this class of exponential integrators by SEI.

For these special exponential integrators, the following properties can be derived.

###### Theorem 3.2
• If the -stage RK method (14) is symmetric, then the -stage SEI (13) is also symmetric.

• The SEI (13) is symplectic if the RK method (14) is symplectic.

• The SEI (13) is symmetric and symplectic if the RK method (14) is symmetric and symplectic.

• If the -stage RK method (14) is of order , then the SEI (13) is also of order .

Proof  Inserting (13) into the symmetry conditions of (5) yields

 ci=1−cs+1−i,aije(ci−cj)hM=ecihMbs+1−je−(1−cs+1−j)hM−as+1−i,s+1−je−(cs+1−i−cs+1−j)hM=(bs+1−j−as+1−i,s+1−j)e(ci−cj)hM,bie(1−ci)hM=ehMbs+1−ie−(1−cs+1−i)hM=bs+1−ie(1−ci)hM,

which can be simplified as the symmetry conditions (10) of RK methods. Thus the first statement is true. The second statement can be obtained immediately by considering Theorem 3.2 of . Based on the above two results, the third one holds. The last result comes from Theorem 3.1 of .

In what follows, based on Theorem 3.2 we construct some practical symmetric and symplectic SEI integrators.

### 3.1 One-stage symmetric and symplectic SEI

First consider an one-stage RK method with the coefficients:

 c1 a11 b1

According to (10) and (12), this method is symmetric and symplectic if

 c1=1/2, a11=b1−a11, b21=2b1a11. (15)

From these formulae, it follows that

 c1=1/2,  a11=12b1. (16)

This gives a class of symmetric and symplectic exponential integrators by considering (13). As an example, we choose and denote the method as SSSEI1s2. It can be checked that this RK method is implicit midpoint rule. Thus the symmetric and symplectic SEI is of order two.

### 3.2 Two-stage symmetric and symplectic SEI

Consider a two-stage RK method whose coefficients are given by a Butcher tableau:

c1 c2 a11 a12 a21 a22 b1 b2

The symmetry conditions of this method are

 c1=1−c2, b1=b2,  a11=b2−a22, a12=b1−a21, a21=b2−a12, a22=b1−a11. (17)

The RK method is symplectic if

 b21=2b1a11,  b1b2=b1a12+b2a21,  b22=2b2a22. (18)

According to (17) and (18), we obtain

 c1=1−c2, b1=b2,  a11=a22=b1/2, a12+a21=b1. (19)

In the light of the third-order conditions of RK methods (see )

 a11+a12=c1,a21+a22=c2, b1+b2=1, b1c1+b2c2=1/2,b1c21+b2c22=1/3,b1(a11c1+a12c2)+b2(a21c1+a22c2)=1/6, (20)

we choose the parameters by

 c1=3−√36, b1=1/2, a12=3−2√312. (21)

This choice as well as (19) and (13) gives an symmetric and symplectic SEI. Moreover, it can be seen that the corresponding RK method is Gauss method of order four (see ). Therefore, the SEI is also of order four, which is denoted by SSSEI2s4.

### 3.3 Three-stage symmetric and symplectic SEI

We turn to considering three-stage symmetric and symplectic integrators. The following Butcher tableau describes a three-stage RK method:

c1 c2 c3 a11 a12 a13 a21 a22 a23 a31 a32 a33 b1 b2 b3

This method is symmetric if the following conditions are true

 c1=1−c3,c2=1/2,b1=b3,b1=a31+a13=a21+a23=a11+a33,b2=a32+a12=2a22,b3=a33+a11=a23+a21=a13+a31. (22)

The RK methods are symplectic if

 b1a11+b1a11=b21,    b2a21+b1a12=b1b2,   b3a31+b1a13=b1b3,b2a22+b2a22=b22,    b3a33+b3a33=b23,      b3a32+b2a23=b2b3. (23)

By the formulae (22) and (23), the coefficients can be given as

c1 12 1−c1 b12 0 0 b1 b22 0 b1 b2 b12 b1 b2 b1

This result as well as (13) yields a class of symmetric and symplectic SEI. As an example and following , we consider

 c1=8−23√2−3√412,  b1=4+23√2+3√46,  b2=−1−23√2−3√43

and denote the method as SSSEI3s4. According to the fourth-order conditions of RK methods (see )

 a11+a12+a13=c1,  a21+a22+a23=c2,  a31+a32+a33=c3,  b1+b2+b3=1,b1c1+b2c2+b3c3=1/2,  b1c21+b2c22+b3c23=1/3,  b1c31+b2c32+b3c33=1/4,b1(a11c1+a12c2+a13c3)+b2(a21c1+a22c2+a23c3)+b3(a31c1+a32c2+a33c3)=1/6,b1c1(a11c1+a12c2+a13c3)+b2c2(a21c1+a22c2+a23c3)+b3c3(a31c1+a32c2+a33c3)=1/8,b1(a11c21+a12c22+a13c23)+b2(a21c21+a22c22+a23c23)+b3(a31c21+a32c22+a33c23)=1/12,b1a11(a11c1+a12c2+a13c3)+b1a12(a21c1+a22c2+a23c3)+b1a13(a31c1+a32c2+a33c3)+b2a21(a11c1+a12c2+a13c3)+b2a22(a21c1+a22c2+a23c3)+b2a23(a31c1+a32c2+a33c3)+b3a31(a11c1+a12c2+a13c3)+b3a32(a21c1+a22c2+a23c3)+b3a33(a31c1+a32c2+a33c3)=1/24, (24)

it can be checked that the coefficients of the RK method satisfy all the conditions. Thus this RK method is of order four and the symmetric and symplectic SEI has same order. This symmetric and symplectic SEI is denoted by SSSEI3s4.

## 4 Numerical experiments

This section presents two numerical experiments to show the remarkable efficiency of the new integrators as compared with some existing RK methods. The integrators for comparisons are chosen as:

• SSSEI1s2: the one-stage symmetric and symplectic EI of order two presented in this letter;

• SSSEI2s4: the two-stage symmetric and symplectic EI of order four presented in this letter;

• SSSEI3s4: the three-stage symmetric and symplectic EI of order four presented in this letter;

• SSRK1s2: the one-stage symmetric and symplectic RK method of order two obtained by letting for SSSEI1s2 (implicit midpoint rule);

• SSRK2s4: the two-stage symmetric and symplectic RK method of order four obtained by letting for SSSEI2s4 (Gauss method of order four);

• SSRK3s4: the three-stage symmetric and symplectic RK method of order four obtained by letting for SSSEI3s4 (the method was given in ).

Problem 1. As the first numerical example, we consider the Duffing equation defined by

 (qp)′=(0               1−ω2−k2   0)(qp)+(02k2q3),   (q(0)p(0))=(0ω).

It is a Hamiltonian system with the Hamiltonian The exact solution of this system is with the Jacobi elliptic function . For this problem, we choose , , and for The efficiency curves are shown in Figure 1 (i). We integrate this problem with a fixed stepsize in the interval for . The results of energy conservation are presented in Figure 1 (ii). Figure 1: (i): The logarithm of the global error (GE) over the integration interval against tend/h. (ii): The logarithm of the maximum global error of Hamiltonian GEH=max|Hn−H0| against log10(tend).

Problem 2. The second numerical example is the following averaged system in wind-induced oscillation

 (x1x2)′=(−ζ−λλ−ζ)(x1x2)+(x1x212(x21−x22)),

where is a detuning parameter and is a damping factor with The first integral (when ) or Lyapunov function (when ) of this system is

 H=12r(x21+x22)−12sin(θ)(x1x22−13x31)+12cos(θ)(−x21x2+13x32).

We choose the initial values Firstly we consider and solve the problem on the interval with for . The efficiency curves are shown in Figure 2 (i). Then this problem is integrated with on the interval See Figure 2 (ii) for the energy conservation for different methods. Secondly we choose and the efficiency curves are shown in Figure 2 (iii) on with Figure 2: (i): The logarithm of the global error (GE) over the integration interval against tend/h. (ii): The logarithm of the maximum global error of Hamiltonian GEH=max|Hn−H0| against log10(tend). (iii): The logarithm of the global error (GE) over the integration interval against tend/h.

From the numerical results, it follows clearly that the symmetric and symplectic exponential integrators behave much better than symmetric and symplectic RK methods.

## 5 Conclusions and discussions

In this letter, in order to solve the differential equations (1) by using symmetric and symplectic methods, we present the symmetry and symplecticity conditions for exponential integrators. Then based on these conditions, we consider a special kind of exponential integrators and construct some practical symmetric and symplectic exponential integrators. The remarkable efficiency of the new integrators is shown by the numerical results from two numerical experiments in comparison with some existing RK methods in the literature.

## References

•  H. Berland, B. Owren, B. Skaflestad, B-series and order conditions for exponential integrators, SIAM J. Numer. Anal. 43 (2005) 1715-1727.
•  M. Caliari, A. Ostermann, Implementation of exponential Rosenbrock-type integrators, Appl. Numer. Math. 59 (2009) 568-581.
•  M.P. Calvo, C. Palencia, A class of explicit multistep exponential integrators for semilinear problems, Numer. Math. 102 (2006) 367-381.
•  E. Celledoni, D. Cohen, B. Owren, Symmetric exponential integrators with an application to the cubic Schrödinger equation, Found. Comput. Math. 8 (2008) 303-317.
•  V. Grimm, M. Hochbruck, Error analysis of exponential integrators for oscillatory second-order differential equations, J. Phys. A: Math. Gen. 39 (2006) 5495-5507.
•  E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd edn. Springer-Verlag, Berlin, Heidelberg, 2006.
•  M. Hochbruck, A. Ostermann, Exponential integrators, Acta Numer. 19 (2010) 209-286.
•  M. Hochbruck, A. Ostermann, J. Schweitzer, Exponential rosenbrock-type methods, SIAM J. Numer. Anal. 47 (2009) 786-803.
•  A. Iserles, On the global error of discretization methods for highly-oscillatory ordinary differential equations, BIT Num. Math. 42 (2002) 561-599.
•  A. Iserles, Think globally, act locally: solving highly-oscillatory ordinary differential equations, Appl. Num. Anal. 43 (2002) 145-160.
•  A.K. Kassam, L.N.Trefethen, Fourth-order time-stepping for stiff PDEs, SIAM J. Sci. Comput. 26 (2005) 1214-1233.
•  M. Khanamiryan, Quadrature methods for highly oscillatory linear and nonlinear systems of ordinary differential equations: part I, BIT Num. Math. 48 (2008) 743-762.
•  S. Krogstad, Generalized integrating factor methods for stiff PDEs, J. Comput. Phys. 203 (2005) 72-88.
•  L. Mei, X. Wu, Symplectic exponential Runge-Kutta methods for solving nonlinear Hamiltonian systems, J. Comput. Phys. 338 (2017) 567–584.
•  A. Ostermann, M. Thalhammer, W.M. Wright, A class of explicit exponential general linear methods, BIT Numer. Math. 46 (2006) 409-431.
•  J.M. Sanz-Serna, L. Abia, Order conditions for canonical Runge–Kutta schemes, SIAM J. Numer. Analy. 28 (1991) 1081-1096.
•  B. Wang, A. Iserles, X. Wu, Arbitrary–order trigonometric Fourier collocation methods for multi-frequency oscillatory systems, Found. Comput. Math. 16 (2016) 151-181.
•  B. Wang, T. Li, Y. Wu, Arbitrary-order functionally fitted energy-diminishing methods for gradient systems, Appl. Math. Lett. 83 (2018) 130-139
•  B. Wang, X. Wu, The formulation and analysis of energy-preserving schemes for solving high-dimensional nonlinear Klein-Gordon equations. IMA. J. Numer. Anal. DOI: 10.1093/imanum/dry047 (2018)
•  B. Wang, H. Yang, F. Meng, Sixth order symplectic and symmetric explicit ERKN schemes for solving multi-frequency oscillatory nonlinear Hamiltonian equations, Calcolo, 54 (2017) 117-140
•  X. Wu, B. Wang, Recent Developments in Structure-Preserving Algorithms for Oscillatory Differential Equations. Springer Nature Singapore Pte Ltd, 2018
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters   