Symplectic Effective Order Numerical Methods for Separable Hamiltonian Systems

Abstract

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.

1. 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,

$$\begin{array}{c}\hfill {\displaystyle \frac{dp}{dt}}=-{\displaystyle \frac{\partial H}{\partial q}},{\displaystyle \frac{dq}{dt}}={\displaystyle \frac{\partial H}{\partial p}},\end{array}$$

(1)

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,

$$\begin{array}{c}\hfill {\displaystyle \frac{dH}{dt}}={\displaystyle \frac{\partial H}{\partial p}}·{\displaystyle \frac{\partial p}{\partial t}}+{\displaystyle \frac{\partial H}{\partial q}}·{\displaystyle \frac{\partial q}{\partial t}}=0.\end{array}$$

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 $w{v}^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.

2. Algebraic Structure of PRK Methods

We are concerned with the numerical solution of separable Hamiltonian systems,

$${\left(\begin{array}{c}u\\ v\end{array}\right)}^{\prime }=\left(\begin{array}{c}k\left(v\right)\\ r\left(u\right)\end{array}\right),u\left(t_0\right)=u_0,v\left(t_0\right)=v_0,k,r:{\mathrm{I}\mathrm{R}}^N\to {\mathrm{I}\mathrm{R}}^N,$$

(2)

using two s-stages RK methods $M=[abc]$ and $\tilde{M}=[\tilde{a}\tilde{b}\tilde{c}]$,

$$\begin{array}{cccccc}\hfill U_i& =u_n+h\sum _{j=1}^{s-1}{a}_{ij}k\left(V_j\right),\hfill & \hfill V_i& =v_n+h\sum _{j=1}^s{\tilde{a}}_{ij}r\left(U_j\right),\hfill & \hfill i& =1,2,\dots ,s,\hfill \\ \hfill u_{n+1}& =u_n+h\sum _{j=1}^s{b}_{j}k\left(V_j\right),\hfill & \hfill v_{n+1}& =v_n+h\sum _{j=1}^s{\tilde{b}}_{j}r\left(U_j\right),\hfill \end{array}$$

where $U_i$ and $V_i$ are the stages for the u and v variables, $b_i$ and ${\tilde{b}}_i$ are quadrature weights, $c_i$ and ${\tilde{c}}_i$ are quadrature nodes, $A={\left(a_{ij}\right)}_{s\times s}$ and $\tilde{A}={\left({\tilde{a}}_{ij}\right)}_{s\times s}$ are matrices of s-stage PRK methods M and $\tilde{M}$, respectively. Here M is an explicit RK method and $\tilde{M}$ is an implicit RK method. Similarly, we can define two PRK methods $S=[ABC]$ and $\tilde{S}=[\tilde{A}\tilde{B}\tilde{C}]$, termed as starting methods. The four PRK methods can be represented by Butcher tableaux as,

$$\begin{array}{cc}\hfill \mathrm{c}& \mathrm{a}\hfill \\ \hfill & b^T\hfill \end{array},\begin{array}{cc}\hfill \tilde{c}& \tilde{a}\hfill \\ & {\tilde{b}}^T\end{array},\begin{array}{cc}\hfill C& A\hfill \\ \hfill & B^T\hfill \end{array},\mathrm{and}\begin{array}{cc}\hfill \tilde{C}& \tilde{A}\hfill \\ & {\tilde{B}}^T\end{array}.$$

We need to find the composed methods $MS$ and $\tilde{M}\tilde{S}$, which are given by the Butcher tableaux,

$$\begin{array}{ccc}c& A& 0\\ C+\sum _{i=1}^s{b}_i& b& A\\ & b^T& B^T\end{array},\begin{array}{ccc}\tilde{c}& \tilde{A}& 0\\ \tilde{C}+\sum _{i=1}^s{\tilde{b}}_i& \tilde{b}& \tilde{A}\\ & {\tilde{b}}^T& {\tilde{B}}^T\end{array}.$$

(3)

The composition asserts to carry out the calculations with starting methods S and $\tilde{S}$ firstly and applying the PRK methods M and $\tilde{M}$ to the output. To explain the composition (3), we consider four PRK methods with two stages each given as,

$$\begin{array}{ccc}\hfill c_1& a_{11}\hfill & \hfill a_{12}\\ \hfill c_2& a_{21}\hfill & \hfill a_{22}\\ \hfill & b_1\hfill & \hfill b_2\end{array},\begin{array}{ccc}\hfill {\tilde{c}}_1& {\tilde{a}}_{11}& {\tilde{a}}_{12}\\ {\tilde{c}}_2& {\tilde{a}}_{21}& {\tilde{a}}_{22}\\ & {\tilde{b}}_1& {\tilde{b}}_2\end{array},\begin{array}{ccc}\hfill C_1& A_{11}\hfill & \hfill A_{12}\\ \hfill C_2& A_{21}\hfill & \hfill A_{22}\\ \hfill & B_1\hfill & \hfill B_2\end{array},\begin{array}{ccc}\hfill {\tilde{C}}_1& {\tilde{A}}_{11}& {\tilde{A}}_{12}\\ {\tilde{C}}_2& {\tilde{A}}_{21}& {\tilde{A}}_{22}\\ & {\tilde{B}}_1& {\tilde{B}}_2\end{array}.$$

Solving the differential equation $u^{{}^{\prime }}=k\left(v\right)$, 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,

$$\begin{array}{cccc}\hfill U_1& =u_0+a_{11}hk\left(V_1\right)+a_{12}hk\left(V_2\right),\hfill & \hfill {\overline{U}}_1& =u_1+A_{11}hk\left({\overline{V}}_1\right)+A_{12}hk\left({\overline{V}}_2\right),\hfill \\ \hfill U_2& =u_0+a_{21}hk\left(V_1\right)+a_{22}hk\left(V_2\right),\hfill & \hfill {\overline{U}}_2& =u_1+A_{21}hk\left({\overline{V}}_1\right)+A_{22}hk\left({\overline{V}}_2\right),\hfill \\ \hfill u_1& =u_0+b_{1}hk\left(V_1\right)+b_{2}hk\left(V_2\right),\hfill & \hfill u_2& =u_1+B_{1}hk\left({\overline{V}}_1\right)+B_{2}hk\left({\overline{V}}_2\right).\hfill \end{array}$$

The composition $MS$ means that we are going from $u_0$ to $u_2$ using the composed method $MS$ given as,

$$\begin{array}{cc}\hfill U_1& =u_0+a_{11}hk\left(V_1\right)+a_{12}hk\left(V_2\right),\hfill \\ \hfill U_2& =u_0+a_{21}hk\left(V_1\right)+a_{22}hk\left(V_2\right),\hfill \\ \hfill {\overline{U}}_1& =u_0+b_{1}hk\left(V_1\right)+b_{2}hk\left(V_2\right)+A_{11}hk\left({\overline{V}}_1\right)+A_{12}hk\left({\overline{V}}_2\right),\hfill \\ \hfill {\overline{U}}_2& =u_0+b_{1}hk\left(V_1\right)+b_{2}hk\left(V_2\right)+A_{21}hk\left({\overline{V}}_1\right)+A_{22}hk\left({\overline{V}}_2\right),\hfill \\ \hfill u_2& =u_0+b_{1}hk\left(V_1\right)+b_{2}hk\left(V_2\right)+B_{1}hk\left({\overline{V}}_1\right)+B_{2}hk\left({\overline{V}}_2\right).\hfill \end{array}$$

(4)

The Butcher table for the above composed PRK methods is as follows

$$\begin{array}{ccccc}{c}_1& a_{11}& a_{12}& 0& 0\\ c_2& a_{21}& a_{22}& 0& 0\\ C_1+\sum b_i& b_1& b_2& A_{11}& A_{12}\\ C_2+\sum b_i& b_1& b_2& A_{21}& A_{22}\\ & b_1& b_2& B_1& B_2\end{array}$$

Similarly, the composition $\tilde{M}\tilde{S}$ can be obtained and is represented by the following Butcher table

$$\begin{array}{ccccc}{\tilde{c}}_1& {\tilde{a}}_{11}& {\tilde{a}}_{12}& 0& 0\\ {\tilde{c}}_2& {\tilde{a}}_{21}& {\tilde{a}}_{22}& 0& 0\\ {\tilde{C}}_1+\sum {\tilde{b}}_i& {\tilde{b}}_1& {\tilde{b}}_2& {\tilde{A}}_{11}& {\tilde{A}}_{12}\\ {\tilde{C}}_2+\sum {\tilde{b}}_i& {\tilde{b}}_1& {\tilde{b}}_2& {\tilde{A}}_{21}& {\tilde{A}}_{22}\\ & {\tilde{b}}_1& {\tilde{b}}_2& {\tilde{B}}_1& {\tilde{B}}_2\end{array}$$

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\left(h^{p+1}\right)$. 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:

$$\begin{array}{cc}\hfill \sum _i{b}_i& =1,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \sum _i{\tilde{b}}_i& =1,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \sum _i{b}_i{\tilde{c}}_i& ={\displaystyle \frac{1}{2}},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \sum _i{\tilde{b}}_i{c}_i& ={\displaystyle \frac{1}{2}},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \sum _i{b}_i{\tilde{c}}_i^2& ={\displaystyle \frac{1}{3}},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \sum _i{\tilde{b}}_i{c}_i^2& ={\displaystyle \frac{1}{3}},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \sum _{i,j}{b}_{i,j}{\tilde{a}}_{ij}{c}_j& ={\displaystyle \frac{1}{6}},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \sum _{i,j}{\tilde{b}}_i{a}_{ij}{\tilde{c}}_j& ={\displaystyle \frac{1}{6}}.\hfill \end{array}$$

(12)

3. 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 $\tilde{t}$. The order of a tree is total number of vertices in a tree. Whereas, the density $\gamma \left(t\right)$ is the product of number of vertices of a tree and their sub-trees, when we remove the root.

Example 1.

Consider a bi-color tree with order 5 and $\gamma \left(t\right)=5\times 4\times 2=40$.

The repeated differentiation of Equation (2) gives

$$\begin{array}{cccc}\hfill u^{\left(1\right)}& =k\left(v\right),\hfill & \hfill v^{\left(1\right)}& =r\left(u\right),\hfill \\ \hfill u^{\left(2\right)}& ={\displaystyle \frac{\partial k}{\partial v}}r,\hfill & \hfill v^{\left(2\right)}& ={\displaystyle \frac{\partial r}{\partial u}}k,\hfill \\ \hfill u^{\left(3\right)}& ={\displaystyle \frac{{\partial }^{2}k}{\partial v\partial v}}(r,r)+{\displaystyle \frac{\partial k}{\partial v}}{\displaystyle \frac{\partial r}{\partial u}}k,\hfill & \hfill v^{\left(3\right)}& ={\displaystyle \frac{{\partial }^{2}r}{\partial u\partial u}}(k,k)+{\displaystyle \frac{\partial r}{\partial u}}{\displaystyle \frac{\partial k}{\partial v}}r,\hfill \\ \hfill ⋮& & \hfill ⋮& \end{array}$$

(13)

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

Table 1. Composition rule for trees with black vertex root.

${\mathit{t}}_{\mathit{i}}$	Tree	$(\mathit{\beta }\mathit{\alpha })({\mathit{t}}_{\mathit{i}})$	$(\mathit{E}\mathit{\beta })({\mathit{t}}_{\mathit{i}})$
$t_1$		${\beta }_1+{\alpha }_1$	${\beta }_1+1$
$t_2$		${\beta }_2+{\tilde{\beta }}_1{\alpha }_1+{\alpha }_2$	${\beta }_2+{\beta }_1+{\displaystyle \frac{1}{2}}$
$t_3$		${\beta }_3+{\tilde{\beta }}_1^2{\alpha }_1+2{\tilde{\beta }}_1{\alpha }_2+{\alpha }_3$	${\beta }_3+2{\beta }_2+{\beta }_1+{\displaystyle \frac{1}{3}}$
$t_4$		${\beta }_4+{\beta }_1{\alpha }_2+{\tilde{\beta }}_2{\alpha }_1+{\alpha }_4$	${\beta }_4+{\beta }_2+{\displaystyle \frac{1}{2}}{\beta }_1+{\displaystyle \frac{1}{6}}$

and

Table 2. Composition rule for trees with white vertex root.

${\mathit{t}}_{\mathit{i}}$	Tree	$(\tilde{\mathit{\beta }}\tilde{\mathit{\alpha }})({\mathit{t}}_{\mathit{i}})$	$(\mathit{E}\tilde{\mathit{\beta }})({\mathit{t}}_{\mathit{i}})$
${\tilde{t}}_1$		${\tilde{\beta }}_1+{\tilde{\alpha }}_1$	${\tilde{\beta }}_1+1$
${\tilde{t}}_2$		${\tilde{\beta }}_2+{\beta }_1{\tilde{\alpha }}_1+{\tilde{\alpha }}_2$	${\tilde{\beta }}_2+{\tilde{\beta }}_1+{\displaystyle \frac{1}{2}}$
${\tilde{t}}_3$		${\tilde{\beta }}_3+{\beta }_1^2{\tilde{\alpha }}_1+2{\beta }_1{\tilde{\alpha }}_2+{\tilde{\alpha }}_3$	${\tilde{\beta }}_3+2{\tilde{\beta }}_2+{\tilde{\beta }}_1+{\displaystyle \frac{1}{3}}$
${\tilde{t}}_4$		${\tilde{\beta }}_4+{\tilde{\beta }}_1{\tilde{\alpha }}_2+{\beta }_2{\tilde{\alpha }}_1+{\tilde{\alpha }}_4$	${\tilde{\beta }}_4+{\tilde{\beta }}_2+{\displaystyle \frac{1}{2}}{\tilde{\beta }}_1+{\displaystyle \frac{1}{6}}$

. The elementary differentials can be represented by bi-color rooted trees as shown in

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

t	Elementary Differentials	$\mathbf{\Phi }\left(\mathit{t}\right)$	$\tilde{\mathit{t}}$	Elementary Differentials	$\mathbf{\Phi }\left(\tilde{\mathit{t}}\right)$
	k	$b_i$		r	${\tilde{b}}_i$
	${\displaystyle \frac{\partial k}{\partial v}}r$	$b_i{\tilde{c}}_i$		${\displaystyle \frac{\partial r}{\partial u}}k$	${\tilde{b}}_i{c}_i$
	${\displaystyle \frac{{\partial }^{2}k}{\partial v\partial v}}(r,r)$	$b_i{\tilde{c}}_i^2$		${\displaystyle \frac{{\partial }^{2}r}{\partial u\partial u}}(k,k)$	${\tilde{b}}_i{c}_i^2$
	${\displaystyle \frac{\partial k}{\partial v}}{\displaystyle \frac{\partial r}{\partial u}}k$	$b_i{\tilde{a}}_{ij}{c}_j$		${\displaystyle \frac{\partial r}{\partial u}}{\displaystyle \frac{\partial k}{\partial v}}r$	${\tilde{b}}_i{a}_{ij}{\tilde{c}}_j$

.

The terms on left hand side from Equations (5)–(12) are called elementary weights $\Phi \left(t\right)$ and $\Phi \left(\tilde{t}\right)$. 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

$$\Phi \left(t\right)={\displaystyle \frac{1}{\gamma \left(t\right)}},\mathrm{and}\Phi \left(\tilde{t}\right)={\displaystyle \frac{1}{\gamma \left(\tilde{t}\right)}},$$

for all bi-color rooted trees t and $\tilde{t}$ of order less than or equal to p [13].

Let the elementary weight functions of the methods M and $\tilde{M}$ are $\alpha$ and $\tilde{\alpha }$, respectively such that $\alpha$ maps trees t and $\tilde{\alpha }$ maps trees $\tilde{t}$ to expression in terms of method coefficients as follows

Similarly, we can define elementary weight functions $\beta$ and $\tilde{\beta }$ for the starting methods S and $\tilde{S}$, respectively. The composition of PRK methods in terms of their elementary weights is $\beta \alpha$ and $\tilde{\beta }\tilde{\alpha }$ given as where the right hand side of Equation (14) contains terms from

Table 4. Calculation for the term $\beta \alpha \left(t_3\right)$.

Tree	No Cut	First Cut	Second Cut	Third Cut	All Cuts
t
u
$t\backslash u$
term	${\beta }_3$	${\tilde{\beta }}_1{\alpha }_2$	${\tilde{\beta }}_1{\alpha }_3$	${\tilde{\beta }}_1^2{\alpha }_1$	${\alpha }_3$

which has tree u as a sub-tree of tree t and $t\backslash u$ is the remaining tree when u is removed from t.

The order condition related to the tree for the composed method (4) is

4. 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 $SM{S}^{-1}$ has order q [6]. For PRK methods, we are interested to construct methods M, $\tilde{M}$, S, $\tilde{S}$ together with the inverse methods $S^{-1}$ and ${\tilde{S}}^{-1}$ such that the compositions $SM{S}^{-1}$ and $\tilde{S}\tilde{M}{\tilde{S}}^{-1}$ have the effective order q which implies $\beta \alpha \left(t\right)=E\beta \left(t\right)$ and $\tilde{\beta }\tilde{\alpha }\left(\tilde{t}\right)=E\tilde{\beta }\left(\tilde{t}\right)$ 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 Table 1 and Table 2.

A method of effective order 3 is obtained by comparing two columns of Table 1 and Table 2 each for trees up to order 3. Thus we have

$$\begin{array}{cc}\hfill {\alpha }_1=& 1,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\tilde{\alpha }}_1=& 1,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\alpha }_2=& {\displaystyle \frac{1}{2}},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\tilde{\alpha }}_2=& \frac{1}{2},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\alpha }_3=& {\displaystyle \frac{1}{3}}+2{\beta }_2,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\tilde{\alpha }}_3=& {\displaystyle \frac{1}{3}}+2{\tilde{\beta }}_2,\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill {\alpha }_4=& {\displaystyle \frac{1}{6}}+{\beta }_2-{\tilde{\beta }}_2,\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill {\tilde{\alpha }}_4=& {\displaystyle \frac{1}{6}}+{\tilde{\beta }}_2-{\beta }_2.\hfill \end{array}$$

(22)

By eliminating $\beta$ and $\tilde{\beta }$ values, Equations (19)–(22) become

$$\begin{array}{cc}\hfill 6{\alpha }_4+3{\tilde{\alpha }}_3-3{\alpha }_3=& 1,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill 6{\tilde{\alpha }}_4+3{\alpha }_3-3{\tilde{\alpha }}_3=& 1.\hfill \end{array}$$

(24)

5. 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 $\tilde{M}$ for separable Hamiltonian system (2) is symplectic if the following condition is satisfied [14]

$$\begin{array}{c}\hfill \mathrm{diag}(b)\tilde{a}+a^T\mathrm{diag}(\tilde{b})-b\tilde{b}=0.\end{array}$$

(25)

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 $\widetilde{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$, $\widetilde{t_3}$ and $\widetilde{t_4}$ are non-superfluous, we can either take ${\alpha }_3$ or $\widetilde{{\alpha }_4}$ and also ${\alpha }_4$ or $\widetilde{{\alpha }_3}$ resulting in reducing Equations (23) and (24) to ${\alpha }_3=\frac{1}{3}$ and $\widetilde{{\alpha }_3}=\frac{1}{3}$. Now consider the following Butcher table for methods M and $\tilde{M}$ which satisfy the symplectic condition (25)

$$\begin{array}{cc}0& 0\\ b_1& b_1\\ b_1+b_2& b_1& b_2\\ & b_1& b_2& b_3\end{array}\begin{array}{cc}{\tilde{b}}_1& {\tilde{b}}_1\\ {\tilde{b}}_1+{\tilde{b}}_2& {\tilde{b}}_1& {\tilde{b}}_2\\ {\tilde{b}}_1+{\tilde{b}}_2+{\tilde{b}}_3& {\tilde{b}}_1& {\tilde{b}}_2& {\tilde{b}}_3\\ & {\tilde{b}}_1& {\tilde{b}}_2& {\tilde{b}}_3\end{array}$$

(26)

The Equations (15)–(17) and (23)–(24) after simplification can be written in terms of elementary weights as

$$\begin{array}{cc}\hfill {\displaystyle \sum _{i=1}^3{b}_i}& =1,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\displaystyle \sum _{i=1}^3{\tilde{b}}_i}& =1,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill b_1{\tilde{b}}_1+b_2({\tilde{b}}_1+{\tilde{b}}_2)+b_3& ={\displaystyle \frac{1}{2}},\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\tilde{b}}_1^2{b}_1+b_2{({\tilde{b}}_1+{\tilde{b}}_2)}^2+b_3{({\tilde{b}}_1+{\tilde{b}}_2+{\tilde{b}}_3)}^2& ={\displaystyle \frac{1}{3}},\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\tilde{b}}_2{b}_1^2+{\tilde{b}}_3{(b_1+b_2)}^2& ={\displaystyle \frac{1}{3}}.\hfill \end{array}$$

(31)

To get the values of 6 unknowns from 5 equations, we have one degree of freedom. Let us take ${\tilde{b}}_1={\displaystyle \frac{2}{3}}$ and solve Equations (27)–(31) to get M and $\tilde{M}$ methods as follows

$$\begin{array}{cccc}& 0& 0& 0\\ & {\displaystyle \frac{13+\sqrt{205}}{12}}& 0& 0\\ & {\displaystyle \frac{13+\sqrt{205}}{12}}& {\displaystyle \frac{5}{6}}& 0\\ & {\displaystyle \frac{13+\sqrt{205}}{12}}& {\displaystyle \frac{5}{6}}& {\displaystyle \frac{-11-\sqrt{205}}{12}}\end{array}\begin{array}{cccc}& {\displaystyle \frac{2}{3}}& 0& 0\\ & {\displaystyle \frac{2}{3}}& {\displaystyle \frac{5+\sqrt{205}}{30}}& 0\\ & {\displaystyle \frac{2}{3}}& {\displaystyle \frac{5+\sqrt{205}}{30}}& {\displaystyle \frac{5-\sqrt{205}}{30}}\\ & {\displaystyle \frac{2}{3}}& {\displaystyle \frac{5+\sqrt{205}}{30}}& {\displaystyle \frac{5-\sqrt{205}}{30}}\end{array}$$

Derivation of Starting Method

For the starting method, we have the following equations

$$\begin{array}{cc}\hfill {\beta }_1=& 0,\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\tilde{\beta }}_1=& 0,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\beta }_2=& {\displaystyle \frac{1}{2}}{\alpha }_3-{\displaystyle \frac{1}{6}},\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill {\tilde{\beta }}_2=& {\displaystyle \frac{1}{2}}{\tilde{\alpha }}_3-{\displaystyle \frac{1}{6}}.\hfill \end{array}$$

(35)

The starting methods should be symplectic [16]. The solution of (32)–(35) leads us to the following symplectic staring PRK methods S and $\tilde{S}$ as

$$\begin{array}{cccc}& 0& 0& 0\\ & {\displaystyle \frac{1}{3}}& 0& 0\\ & {\displaystyle \frac{1}{3}}& {\displaystyle \frac{2}{5}}& 0\\ & {\displaystyle \frac{1}{3}}& {\displaystyle \frac{2}{5}}& -{\displaystyle \frac{11}{15}}\end{array},\begin{array}{cccc}& {\displaystyle \frac{1}{3}}& 0& 0\\ & {\displaystyle \frac{1}{3}}& -{\displaystyle \frac{11}{18}}& 0\\ & {\displaystyle \frac{1}{3}}& -{\displaystyle \frac{11}{18}}& {\displaystyle \frac{5}{18}}\\ & {\displaystyle \frac{1}{3}}& -{\displaystyle \frac{11}{18}}& {\displaystyle \frac{5}{18}}\end{array}.$$

6. 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_1={\tilde{b}}_3$, $b_2={\tilde{b}}_2$ and $b_3={\tilde{b}}_1$ in Equations (27)–(31), so we have

$$\begin{array}{cc}\hfill {\displaystyle \sum _{i=1}^3{\tilde{b}}_i}& =1,\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill {\displaystyle \sum _{i=1}^3{\tilde{b}}_i}& =1,\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill {\tilde{b}}_3{\tilde{b}}_1+{\tilde{b}}_2({\tilde{b}}_1+{\tilde{b}}_2)+{\tilde{b}}_1& ={\displaystyle \frac{1}{2}},\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill {\tilde{b}}_2{\tilde{b}}_3^2+{\tilde{b}}_3{({\tilde{b}}_3+{\tilde{b}}_2)}^2& ={\displaystyle \frac{1}{3}}.\hfill \end{array}$$

(39)

Sanz-Serna suggested in [14] to take ${\tilde{b}}_3$ = 0.91966152, which leads us to the following main methods M and $\tilde{M}$ with effective order 3 with just 3 stages

$$\begin{array}{cccc}& 0& 0& 0\\ & 0.91966152& 0& 0\\ & 0.91966152& -0.18799162& 0\\ & 0.91966152& -0.18799162& 0.26833010\end{array},\begin{array}{cccc}& 0.26833010& 0& 0\\ & 0.26833010& -0.18799162& 0\\ & 0.26833010& -0.18799162& 0.91966152\\ & 0.26833010& -0.18799162& 0.91966152\end{array}.$$

The starting methods S and $\tilde{S}$ for mutually adjoint symplectic effective order PRK method are constructed by using $B_1={\tilde{B}}_3$, $B_2={\tilde{B}}_2$, $B_3={\tilde{B}}_1$ in Equations (32) to (34) to get

$$\begin{array}{cc}{\tilde{B}}_1+{\tilde{B}}_2+{\tilde{B}}_3\hfill & =0,\hfill \end{array}$$

(40)

$$\begin{array}{cc}{\tilde{B}}_1{\tilde{B}}_3+{\tilde{B}}_2({\tilde{B}}_1+{\tilde{B}}_2)+{\tilde{B}}_1({\tilde{B}}_1+{\tilde{B}}_2+{\tilde{B}}_3)\hfill & \simeq 0.\hfill \end{array}$$

(41)

By solving Equations (40) and (41) with ${\tilde{B}}_3={\displaystyle \frac{1}{2}}$, we get the following S and $\tilde{S}$ methods:

$$\begin{array}{cccc}& 0& 0& 0\\ & {\displaystyle \frac{1}{2}}& 0& 0\\ & {\displaystyle \frac{1}{2}}& {\displaystyle -\frac{1}{4}}& 0\\ & {\displaystyle \frac{1}{2}}& -{\displaystyle \frac{1}{4}}& -{\displaystyle \frac{1}{4}}\end{array}\mathrm{and}\begin{array}{cccc}& -{\displaystyle \frac{1}{4}}& 0& 0\\ & -{\displaystyle \frac{1}{4}}& -{\displaystyle \frac{1}{4}}& 0\\ & -{\displaystyle \frac{1}{4}}& -{\displaystyle \frac{1}{4}}& {\displaystyle \frac{1}{2}}\\ & -{\displaystyle \frac{1}{4}}& -{\displaystyle \frac{1}{4}}& {\displaystyle \frac{1}{2}}\end{array}$$

7. 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.

7.1. Kepler’s Two Body Problem

Consider Kepler’s two body problem given as

$$\begin{array}{cc}\hfill u_1^{{}^{\prime }}=v_1,& u_1\left(0\right)=1-e,\hfill \\ \hfill u_2^{{}^{\prime }}=v_2,& u_2\left(0\right)=0,\hfill \\ \hfill v_1^{{}^{\prime }}={\displaystyle \frac{-u_1}{r^3}},& v_1\left(0\right)=0,\hfill \\ \hfill v_2^{{}^{\prime }}={\displaystyle \frac{-u_2}{r^3}},& v_2\left(0\right)=\sqrt[]{{\displaystyle \frac{1+e}{1-e}}}.\hfill \end{array}$$

where $r=\sqrt[]{u_1^2+u_2^2}$. The energy is given by,

$$\begin{array}{c}\hfill H={\displaystyle \frac{1}{2}}(v_1^2+v_2^2)-{\displaystyle \frac{1}{\sqrt{u_1^2+u_2^2}}}.\end{array}$$

The exact solution after half revolution is

$$\begin{array}{ccc}\hfill u_1\left(\pi \right)=1+e,& u_2\left(\pi \right)=0,v_1\left(\pi \right)=0,\hfill & \hfill v_2\left(\pi \right)=\sqrt{{\displaystyle \frac{1-e}{1+e}}}.\end{array}$$

To verify the effective order 3 behavior, we apply the starting methods S and $\tilde{S}$ to perturb initial values to ${\left({\tilde{u}}_1\right)}_0$, ${\left({\tilde{u}}_2\right)}_0$ and ${\left({\tilde{v}}_1\right)}_0$, ${\left({\tilde{v}}_2\right)}_0$, respectively. Then we apply the main methods M and $\tilde{M}$ for n number of iterations to ${\left({\tilde{u}}_1\right)}_0$, ${\left({\tilde{u}}_2\right)}_0$ and ${\left({\tilde{v}}_1\right)}_0$, ${\left({\tilde{v}}_2\right)}_0$, respectively and obtain the numerical solutions ${\left({\tilde{u}}_1\right)}_n$, ${\left({\tilde{u}}_2\right)}_n$ and ${\left({\tilde{v}}_1\right)}_n$, ${\left({\tilde{v}}_2\right)}_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\left(t_n\right)$, $u_2\left(t_n\right)$ and $v_1\left(t_n\right)$, $v_2\left(t_n\right)$ and perturb them using starting methods S and $\tilde{S}$ to obtain ${\tilde{u}}_1\left(t_n\right)$, ${\tilde{u}}_2\left(t_n\right)$ and ${\tilde{v}}_1\left(t_n\right)$, ${\tilde{v}}_2\left(t_n\right)$. Finally, to obtain global error, we take the difference between numerical and exact solutions, i.e., $\left||{\tilde{u}}_n-\tilde{u}\left(t_n\right)|\right|$. Effective order 3 behavior for symplectic and mutually adjoint symplectic PRK is confirmed from

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

h	n	Global Error	Ratio
${\displaystyle \frac{\pi }{225}}$	225	$7.7741637102284\times {10}^{-04}$
			$8.927465$
${\displaystyle \frac{\pi }{450}}$	450	$8.7081425338124\times {10}^{-05}$
			$8.475956$
${\displaystyle \frac{\pi }{900}}$	900	$1.02739357736254\times {10}^{-05}$
			$8.063877$
${\displaystyle \frac{\pi }{1800}}$	1800	$1.27406905909534\times {10}^{-06}$

and

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

h	n	Global Error	Ratio
${\displaystyle \frac{\pi }{40}}$	40	$4.635890382086\times {10}^{-03}$
			$7.963129$
${\displaystyle \frac{\pi }{80}}$	80	$5.82169473987045\times {10}^{-04}$
			$8.062262$
${\displaystyle \frac{\pi }{160}}$	160	$7.22091971134495\times {10}^{-05}$
			$7.836579$
${\displaystyle \frac{\pi }{320}}$	320	$9.21437776243649\times {10}^{-06}$

.

We calculate the ratio of global errors calculated at step length h, ${\displaystyle \frac{h}{2}}$, ${\displaystyle \frac{h}{4}}$ and ${\displaystyle \frac{h}{8}}$. For method of order p, ratio should approximately be $2^p$ [12].

The next experiment is to verify the energy conservation behaviour of the symplectic methods. In this experiment, the step-size $h=2\pi /1000$ is used for ${10}^5$ iterations. Figure 1 and Figure 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.

7.2. Harmonic Oscillator

The motion of a unit mass attached to a spring with momentum u and position co-ordinates v defines the Hamiltonian system

$$v^{{}^{\prime }}=u{u}^{{}^{\prime }}=v$$

The energy is given by

$$H={\displaystyle \frac{u^2}{2}}+\frac{v^2}{2}$$

The exact solution is

$$\left[\begin{array}{c}u\left(t\right)\\ v\left(t\right)\end{array}\right]=\left[\begin{array}{cc}cos\left(t\right)& -sin\left(t\right)\\ sin\left(t\right)& cos\left(t\right)\end{array}\right]\left[\begin{array}{c}u\left(0\right)\\ v\left(0\right)\end{array}\right]$$

We have applied the symplectic PRK and mutually adjoint symplectic PRK methods with step-size $h=2\pi /1000$ with ${10}^5$ iteration in this experiment. We have obtained good energy conservation as shown in Figure 3 and Figure 4 by symplectic effective order PRK and mutually adjoint symplectic effective order PRK methods, respectively.

8. 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.