Symplectic Effective Order Numerical Methods for Separable Hamiltonian Systems

A family of explicit symplectic partitioned Runge-Kutta methods are derived with effective order 3 for the numerical integration of separable Hamiltonian systems. The proposed explicit methods are more efficient than existing symplectic implicit Runge-Kutta methods. A selection of numerical experiments on separable Hamiltonian system confirming the efficiency of the approach is also provided with good energy conservation.


Introduction
Over the last few decades, a lot of progress has been made in developing numerical methods for ordinary differential equations (ODEs), which can produce efficient, reliable and qualitatively correct numerical solutions by preserving some qualitative features of the exact solutions [1,2].This field of research is called geometric numerical integration.In this paper, we are oriented to obtain numerical approximation of Hamiltonian differential equations of the form, where H is known as the Hamiltonian or the total energy of the system and q and p are generalized coordinates and generalized momenta respectively.The autonomous Hamiltonian systems have two important properties; one is that the total energy remains constant, The second important property is that the phase flow is symplectic which imply that the motion along the phase curve retains the area of a bounded sub-domain in the phase space.We need such numerical methods which can mimic both properties of the Hamiltonian systems.For this we use symplectic numerical methods to solve system of Equation (1).The symplectic methods are numerically more efficient than non-symplectic methods for integration over long interval of time [3].Among the class of one step symplectic methods, the implicit symplectic Runge-Kutta (RK) methods were developed and presented in [4].
For general explicit RK methods, it is known that up to order 4, number of stages required in a RK method are equal to the order of the method whereas, for order ≥5, number of stages become greater than the order of the method.Thus for example, an RK method of order 5 needs at least 6 stages and an RK method of order 6 requires at least 17 stages and so on.The more the number of stages are, the more costly the method is.
Butcher [5] tried to overcome this complexity of order barrier by presenting the idea of effective order.He implemented his idea on RK method of order 5 and was able to construct explicit RK methods of effective order 5 with just 5 internal stages [5].Later on, the idea was extended to construct diagonally extended singly implicit RK methods for the numerical integration of stiff differential equations [6][7][8][9].Butcher also used the effective order technique on symplectic RK methods for the numerical integration of Hamiltonian systems [10,11].
A RK method v has an effective order q if we have another RK method w, known as the starting method, such that wv n w −1 has an order q.The method w is used only once in the beginning followed by n iterations of the main method v and the method w −1 , known as the finishing method is used only once at the end.
For separable Hamiltonian systems, it is advantageous to solve some components of the differential-equation with one RK method and solve other components of differential equation with another RK method and collectively they are termed as partitioned Runge-Kutta (PRK) methods.We have extended the idea of Butcher to construct symplectic effective order PRK methods which are explicit in nature and hence are less costly than symplectic implicit RK methods.For the effective order PRK methods, we construct two main PRK methods, two starting and two finishing PRK methods such that the two starting methods are applied once at the beginning followed by n iterations of the main PRK methods and the two finishing methods are used at the end.

Algebraic Structure of PRK Methods
We are concerned with the numerical solution of separable Hamiltonian systems,  .
We need to find the composed methods MS and M S , which are given by the Butcher tableaux, The composition asserts to carry out the calculations with starting methods S and S firstly and applying the PRK methods M and M to the output.To explain the composition (3), we consider four PRK methods with two stages each given as, , Solving the differential equation u = k(v), we take the first step going from u 0 to u 1 using starting method S. We then take the second step going from u 1 to u 2 using main method M. The equations are given as, The composition MS means that we are going from u 0 to u 2 using the composed method MS given as, The Butcher table for the above composed PRK methods is as follows Similarly, the composition M S can be obtained and is represented by the following Butcher table A numerical scheme has order p, if after one iteration, the numerical solution matches with the Taylor's series of the exact solution to the extent that the remainder term has O(h p+1 ).The comparison of the numerical solution with the Taylor's series of the exact solution provides order conditions which must be satisfied to obtain a particular order numerical method.Thus for order 3 method, we have the following order conditions:

Rooted Trees for PRK Methods and Order Conditions
There is a graph theoretical approach to study order conditions of RK methods due to Butcher [12].We first study the basic concepts from graph theory and then relate the graphs to the order conditions of RK methods.
A graph with non-cyclic representation consisting of vertices and edges with one vertex acting as root is called as rooted tree.The rooted trees whose vertices are either black or white are known as bi-color rooted trees.The trees having black color root is termed as t, while the white rooted trees are represented by t.The order of a tree is total number of vertices in a tree.Whereas, the density γ(t) is the product of number of vertices of a tree and their sub-trees, when we remove the root.
The repeated differentiation of Equation (2) gives =r(u), The terms on the right hand side of Equation ( 13) are called elementary differentials and can be represented graphically with the help of bi-color rooted trees.Thus k is represented by a black vertex and r is represented by a white vertex, whereas the differentiation is represented by an edge between two vertices.According to the nature of the differential Equation (2) in which k depends only on v and r depends only on u, we only consider trees in which a black vertex has a white child and vice versa [13].Such trees are given in Tables 1 and 2. The elementary differentials can be represented by bi-color rooted trees as shown in Table 3.
Table 1.Composition rule for trees with black vertex root.
Table 2. Composition rule for trees with white vertex root.
The terms on left hand side from Equations ( 5)-( 12) are called elementary weights Φ(t) and Φ( t).The elementary weights are nonlinear expressions of the coefficients of PRK methods and can be related to bi-color rooted trees as shown in Table 3 [12,13].A PRK method is of order p iff for all bi-color rooted trees t and t of order less than or equal to p [13].
Let the elementary weight functions of the methods M and M are α and α, respectively such that α maps trees t and α maps trees t to expression in terms of method coefficients as follows α( Similarly, we can define elementary weight functions β and β for the starting methods S and S, respectively.The composition of PRK methods in terms of their elementary weights is βα and β α given as where the right hand side of Equation ( 14) contains terms from Table 4 which has tree u as a sub-tree of tree t and t \ u is the remaining tree when u is removed from t.
The order condition related to the tree for the composed method (4) is = βα( ).

Effective Order PRK Methods
Let M and S be two RK methods together with an inverse method S −1 .The effective order q means that the composition SMS −1 has order q [6].For PRK methods, we are interested to construct methods M, M, S, S together with the inverse methods S −1 and S −1 such that the compositions SMS −1 and S M S −1 have the effective order q which implies βα(t) = Eβ(t) and β α( t) = E β( t) for all trees of order q [12], where E is the exact flow for which corresponding order conditions are satisfied.The terms are given in Tables 1 and 2.
A method of effective order 3 is obtained by comparing two columns of Tables 1 and 2 each for trees up to order 3. Thus we have By eliminating β and β values, Equations ( 19)-( 22) become

Symplectic PRK Methods with Effective Order 3
The flow of a Hamiltonian system ( 2) is symplectic and it is a well known fact that the discrete flow by symplectic Runge-Kutta methods is symplectic [14].The PRK method M and M for separable Hamiltonian system (2) is symplectic if the following condition is satisfied [14] diag Moreover, the composition of two symplectic RK methods is symplectic [10,15].For symplectic PRK methods, the trees related to order conditions can be divided into superfluous and non-superfluous bi-color trees.Unlike RK methods, the superfluous bi-color trees of PRK methods also contribute one order condition together with one condition from non-superfluous bi-color tree [14].Since, the underlying bi-color tree of the rooted trees t 2 and t 2 is superfluous.So, we can ignore the condition (18) because it is automatically satisfied.Moreover, the underlying bi-color trees of t 3 , t 4 , t 3 and t 4 are non-superfluous, we can either take α 3 or α 4 and also α 4 or α 3 resulting in reducing Equations ( 23) and (24) to α 3 = 1 3 and α 3 = 1 3 .Now consider the following Butcher table for methods M and M which satisfy the symplectic condition (25) The Equations ( 15)-( 17) and ( 23)-( 24) after simplification can be written in terms of elementary weights as To get the values of 6 unknowns from 5 equations, we have one degree of freedom.Let us take

Derivation of Starting Method
For the starting method, we have the following equations The starting methods should be symplectic [16].The solution of (32)-( 35) leads us to the following symplectic staring PRK methods S and S as .

Mutually Adjoint Symplectic Effective Order PRK Methods
A separable Hamiltonian system remains unchanged by changing the role of kinetic and potential energies, position, and momentum and also inverting the direction of time.For the two PRK method tableaux (26) being mutually adjoint, we have b Sanz-Serna suggested in [14] to take b 3 = 0.91966152, which leads us to the following main methods M and M with effective order 3 with just 3 stages The starting methods S and S for mutually adjoint symplectic effective order PRK method are constructed by using 32) to (34) to get By solving Equations ( 40) and (41) with B 3 = 1 2 , we get the following S and S methods: and

Numerical Experiments
The symplectic RK methods should be implicit and hence their computational cost is higher due to large number of function evaluations.On the other hand, we can use symplectic explicit PRK methods for Hamiltonian systems with lesser computational cost because of the explicit nature of the stages.Our derived effective order symplectic PRK methods are explicit and we have used MATLAB to implement them on Hamiltonian systems for the energy conservation and order confirmation of these methods.

Kepler's Two Body Problem
Consider Kepler's two body problem given as where r = u 2 1 + u 2 2 .The energy is given by, The exact solution after half revolution is To verify the effective order 3 behavior, we apply the starting methods S and S to perturb initial values to ( u 1 ) 0 , u 2 ) 0 and ( v 1 ) 0 , ( v 2 ) 0 , respectively.Then we apply the main methods M and M for n number of iterations to ( u 1 ) 0 , ( u 2 ) 0 and ( v 1 ) 0 , ( v 2 ) 0 , respectively and obtain the numerical solutions ( u 1 ) n , ( u 2 ) n and ( v 1 ) n , ( v 2 ) n computed at t n = t 0 + nh where h is the step-size.We then evaluate exact solutions at t n to get u 1 (t n ), u 2 (t n ) and v 1 (t n ), v 2 (t n ) and perturb them using starting methods S and S to obtain u 1 (t n ), u 2 (t n ) and v 1 (t n ), v 2 (t n ).Finally, to obtain global error, we take the difference between numerical and exact solutions, i.e., || u n − u(t n )||.Effective order 3 behavior for symplectic and mutually adjoint symplectic PRK is confirmed from Tables 5 and 6.The next experiment is to verify the energy conservation behaviour of the symplectic methods.In this experiment, the step-size h = 2π/1000 is used for 10 5 iterations.Figures 1 and 2 depict good conservation of energy with symplectic and mutually adjoint symplectic effective order PRK methods, respectively.Whereas, the energy error was bounded above by 10 −13 and 10 −14 , respectively.

Harmonic Oscillator
The motion of a unit mass attached to a spring with momentum u and position co-ordinates v defines the Hamiltonian system The energy is given by The exact solution is u(t) v(t) = cos(t) − sin(t) sin(t) cos(t) We have applied the symplectic PRK and mutually adjoint symplectic PRK methods with step-size h = 2π/1000 with 10 5 iteration in this experiment.We have obtained good energy conservation as shown in Figures 3 and 4 by symplectic effective order PRK and mutually adjoint symplectic effective order PRK methods, respectively.

Conclusions
In this paper, the composition of two PRK methods is elaborated in terms of Butcher tableaux as well as in terms of their elementary weight functions.The effective order conditions are provided for PRK methods and 3 stage effective order 3 symplectic PRK methods are constructed and successfully applied to separable Hamiltonian systems with good energy conservation.It is worth mentioning that we are able to construct mutually adjoint symplectic effective order 3 PRK methods with just 3 stages, whereas an equivalent method of order 3 with 4 stages is given in [16].In dynamical solar systems, we deal with many problems which are described by Hamiltonian systems, like, Kepler's two body problem and more realistic Jovian five body problem.These symplectic methods are much useful for such types of physical phenomena to observe their positions and energy conservation.
using two s-stages RK methods M = [a b c] and M = [ ã b c], where U i and V i are the stages for the u and v variables, b i and b i are quadrature weights, c i and c i are quadrature nodes, A = (a ij ) s×s and A = ( a ij ) s×s are matrices of s-stage PRK methods M and M, respectively.Here M is an explicit RK method and M is an implicit RK method.Similarly, we can define two PRK methods S = [A B C] and S = [ Ã B C], termed as starting methods.The four PRK methods can be represented by Butcher tableaux as,

c 1 a 11 a 12 c 2 a 21 a 22 b 1 b 2 , c 1 a 11 a 12 c 2 a 21 a 22 b 1 b 2
(27)-(31) to get M and M methods as follows

Figure 4 .
Figure 4. Energy error for the Harmonic oscillator problem (e = 0) with mutually adjoint symplectic effective PRK method using step size h = 2π/1000 for 10 5 steps.

Table 3 .
Elementary differentials and elementary weights of bi-color rooted trees.

Table 5 .
Global errors and their comparison: Symplectic Effective order PRK method.

Table 6 .
Global errors and their comparison: Mutually adjoint symplectic Effective order PRK method.