Explicit Symplectic Runge–Kutta–Nyström Methods Based on Roots of Shifted Legendre Polynomial

: To date, all explicit symplectic Runge–Kutta–Nyström methods of order ﬁve or above are derived by numerical solutions of order condition equations and symplectic condition. In this paper, we derive 124 sets of seven-stage ﬁfth-order explicit symplectic Runge–Kutta–Nyström methods with closed-form coefﬁcients in the Butcher tableau using the roots of a degree-3 shifted Legendre polynomial. One method is analyzed and its P-stable interval is derived. Numerical tests on the two newly discovered methods are performed, showing their long-time stability and large step size stability over some existing methods.


Introduction
In recent decades, geometric numerical methods that preserve at least one geometric property of a given system have received significant attention in various fields, such as astrophysics, computational fluid dynamics and molecular dynamics; see monographs [1][2][3][4][5] and references therein.For other traditional numerical methods for ordinary differential equations, such as Runge-Kutta methods, we refer to the book [6].
Consider the classical second-order differential equation [2,4] q = − U(q), q(t 0 ) = q 0 , q(t 0 ) = p 0 , where U : R d → R is assumed to be sufficiently smooth.Introducing p = q, wherein q and p are vectors, represents a function related to potential energy.An s-stage explicit symplectic Runge-Kutta-Nyström (ESRKN) method applied to the above equation is [7] where h is the step size, bi = b i (1 − c i ), i = 1, . . ., s and āij = 0, i ≤ j, b j (c i − c j ), i > j.
The ESRKN method (1) is expressed in a Butcher tableau by where the symbols match their locations in (1).For more information about Butcher tableaus, we refer to a Butcher book [6].
Therefore, an ESRKN method is determined only by nodes c i and quadrature coefficients b i .There are different strategies used to obtain c i and b i ; we mention (without being exhaustive) the following references [5,[8][9][10][11][12][13][14][15][16].For ESRKN methods with an order higher than four, none of the aforementioned works present c i and b i in exact closed form.In other words, those methods obtained do not satisfy order conditions and symplectic conditions exactly, but within some digits of accuracy.Hence, an interesting question is whether or not there exists a closed form of c i and b i for ESRKN with an order higher than four.In this paper, we present a positive answer to this question.
In an earlier work, authors in [8] constructed a fifth-order ESRKN method by solving 13 equations of order conditions, while constructing a sixth-order method requires solving 23 equations of order conditions.Obviously, as the order increases, so does the computational cost of solving the nonlinear order equations to produce c i and b i .Although the research by Calvo and Hairer has significantly reduced the number of independent order conditions in the ESRKN method [17], the coupled order conditions equations still comprise a considerable proportion of the total order conditions equations.On the other hand, despite modern computers' tremendous performance, it is nearly impossible to solve order conditions directly to obtain c i and b i with a higher order.One possibility is to assign some suitable c i first and then solve order conditions to yield b i .But what is the preassigned acceptable c i ?It is worth noting that the c i nodes of both three-stage fourth-order ESRKN methods in Tables 1 and 2 are the roots of the degree-2 shifted Legendre polynomial [18].
Numerous researchers have already successfully utilized polynomials or functions to design numerical methodologies.For instance, Revelli and Ridolfi [19] developed an interpolation collocation method based on the sine function.Additionally, a fourth-order algorithm was designed based on sinc-collocation in [20].The successful instances have sparked a great deal of interest in our investigation for constructing higher-order ESRKN methods based on shifted Legendre polynomials.
We investigate the potential of using the degree-3 shifted Legendre polynomial roots as the c i nodes and come up with a five-stage fourth-order ESRKN method, which further confirms that it is feasible to construct ESRKN using the shifted Legendre polynomial roots as the c i nodes.We then conduct further research to derive fifth-order ESRKN methods.The main contribution of this paper is to show that this is an effective technique for deriving ESRKN methods and we demonstrate this by deriving 124 seven-stage fifth-order ESRKN methods with c i and b i in closed form.To the best of our knowledge, this is the first work carried out to derive closed-form ESRKN methods with order five.We will explore closed-form ESRKN methods with order 6 or higher in a future paper.
This paper is organized as follows.In Section 2, we introduce order conditions for the ESRKN method (1) and demonstrate two closed-form ESRKN methods with order four and five, respectively.Their P-stable intervals are also provided.Section 3 presents numerical experiments.In Section 4, some conclusions are drawn.Appendix A lists 124 sets of preassigned acceptable c i , where the corresponding closed-form real solution b i could be determined by order conditions.

ESRKN Methods with Order 4 and 5 2.1. Order Conditions for ESRKN Method
Given the assumption bi = b i (1 − c i ), the order conditions for RKN method (1) are as follows [8,21]: If only condition (1) is satisfied, the method is of order one.If both conditions (1) and ( 2) are satisfied, the method is of order two.If conditions (1) to (4) are satisfied, the method is of order three.If conditions (1) to (7) are satisfied, the method is of order four, and if conditions (1) to ( 13) are satisfied, the method is of order five.It has been demonstrated that, for ESRKN methods, conditions ( 7), ( 12) and ( 13) are redundant [8].

A Five-Stage Fourth-Order ESRKN Method
Let the degree-l normalized shifted Legendre polynomial P l (x) be [21], The roots of degree-3 shifted Legendre polynomial By placing g 1 , g 2 , g 3 in five positions as a set of c i , we find a closed-form five-stage fourth-order ESRKN method whose Butcher tableau is presented in Table 3. Theorem 1.The explicit five-stage method defined in Table 3 is a fourth-order symplectic RKN method and its P-stable interval [22] is (0, 7.75342...).

A Seven-Stage Fifth-Order ESRKN Method
We choose the stage number as seven and place g 1 , g 2 and g 3 of (3) in seven positions of {c i }.There are 3 7 = 2187 permutations for c i .By solving the fifth-order conditions of 13 equations, we find real solution b i for 124 sets of c i .For each set of nodes c i , we present its adjoint counterpart in Appendix A. Therefore, we find in total 124 seven-stage fifth-order ESRKN methods in closed form, where c i is chosen among g 1 , g 2 and g 3 of (3) and b i is solved by fifth-order conditions in (2).In particular, we display in Table 4 such an ESRKN method.Here, in Theorem 2. The seven-stage explicit method defined in Table 4 is a fifth-order symplectic RKN method and its P-stable interval is (0, 9.22575...).
Proof.We verify the 13 equations in (2) by plugging the coefficients c i , b i , bi and āij from Table 4.We compute the spectrum of M in (4) with c, b, b and ā from Table 4 to obtain

Numerical Experiments
We compute two numerical examples.The first example has a known solution that the computation verifies, determining the symplecticity-preserving and the order convergence of the new methods.The second example shows a strong stability of the new methods over some existing methods.
We solve the Kepler problem, cf.[4], The Hamiltonian is The exact solution is , where e = 3/10, T = 2πa 3/2 , a = 40/7 and E(t) is the solution of the Kepler equation We solve the Kepler problem in t ∈ [0, 5T] using the five-stage fourth-order method defined in Table 3, and using the seven-stage fifth-order method defined in Table 4, with grid size h = 5T/2 7 , . . ., 5T/2 12 .In Figure 1, we plot the error of the computed Hamiltonian (6).It indicates that both methods are symplectic.The error is roughly caused by the computer round-off.The error of computed Hamiltonian, H(q, q) − H(q 0 , q0 ), on [0, 5T], using the five-stage fourth-order method (top) and using the seven-stage fifth-order method (bottom).
In Table 5, we list the computed error at t = 5T and the computed order of convergence using the fifth-stage fourth-order method in Table 3, and using the seven-stage fifth-order method defined in Table 4.
Table 5. Error profile using the 5-stage 4th order method in Table 3 (left), and using the 7-stage 5th order method in Table 4 (right).We next solve the Hénon-Heiles Hamiltonian system [23,24] with initial values (p 1 , p 2 , q 1 , q 2 ) = (0, 0, 0.1, −0.5).The corresponding Hamiltonian function is denoted by In the following experiments, we apply four ESRKN methods to (7): the five-stage fourth-order ESRKN method listed in Table 3; the seven-stage fifth-order ESRKN method listed in Table 4; the five-stage fifth-order ESRKN method from [8] listed in the first part of Table 6, which satisfies the order condition with an error tolerance of 10 −15 ; the sevenstage fifth-order RKN method from [10] listed in the second part of Table 6, which satisfies the order condition with an error tolerance of less than 3×10 −16 .In short, we will denote these four methods as NEW4, NEW5, OS5 and CS5, cf. the abbreviation table at the end.In all the following figures, the cyan, red, magenta and blue lines denote the results derived by NEW4, NEW5, OS5 [8] and CS5 [10], respectively.Table 6.Top: the five-stage fifth-order method (OS5) from [8].Bottom: the seven-stage fifth-order method (CS5) from [10].In Table 7, we present the convergence rates of the NEW4, NEW5, OS5 [8] and CS5 [10] methods, where global errors for q and p corresponding to these four methods at five step sizes h in the time interval [0, 1] are provided.It verifies the convergence rates of all four methods.The research findings from [25] revealed that particles of (7) manifest chaotic behavior when the system energy H > 1 12 .Furthermore, it was proposed that when H is less than 1  6 , particle trajectories are confined within the equilateral triangle defined by H = 1 6 .For the sake of this study, an initial value of H = 1 6 is selected, resulting in chaotic numerical behavior while constraining trajectories within the triangular region.This property is expected to be retained by symplectic numerical methods.Numerical computations have demonstrated, however, that those methods such as OS5 [8] and CS5 [10], which do not fully meet the order condition equations and symplectic condition, are incapable of retaining long-term stability.After a certain time, the variations in particle trajectories escalate rapidly, eventually exceeding the boundaries of the triangular region.
Figure 2 illustrates the trajectory plots derived from numerical simulations of the Hénon-Heiles system within the interval [0, 1000] using a step size of h = 0.5.All four methods display trajectories that are confined within the equilateral triangle area in this situation.Nonetheless, Figure 3 shows that when the OS5 [8] method is used to generate numerical trajectories within the interval [0, 17,575], the trajectories escape directly from the upper region.Figure 4 looks into the causes of the OS5 [8] method's escape phenomena, attributing them to the numerical inability to maintain system energy.The capacity of symplectic algorithms to maintain correct numerical behavior over lengthy periods of time is a significant advantage.Therefore, we continue to increase the time interval.In Figure 5, it is evident that the NEW4, NEW5 and CS5 [10] methods maintain accurate numerical behavior even at t = 100,000.With such extensive time intervals, the need for further comparisons with bigger time intervals appears to be decreasing for the NEW4, NEW5 and CS5 [10] methods.As a result, the step size is steadily increased in the interval [0, 10,000], beginning with h = 0.50 and gradually increasing until it reaches h = 0.86.The CS5 [10] method begins to gradually deviate from the triangular region at this point.However, the NEW4 and NEW5 methods continue to display reliable numerical performance in the face of such variances.These phenomena are illustrated in Figure 6.Finally, Figure 7 depicts the variation in Hamiltonian energy over time using the NEW4, NEW5 and CS5 [10] methods, emphasizing that the inability to maintain system energy is the cause of the escape phenomena.

Discussion
This article introduces an approach to constructing ESRKN methods that precisely satisfy the order conditions and the symplectic conditions, resulting in 124 sets of sevenstage fifth-order ESRKN methods.Based on the results of numerical experiments, it is clear that using ESRKN methods that exactly meeting the order conditions and the symplectic conditions has advantages over approaches that have inherent errors in satisfying these criteria.In light of these findings, a future research will be dedicated to developing sixthorder ESRKN methods.

Figure 2 .
Figure 2. Solving the numerical orbit diagram of the Hénon-Heiles system over the interval [0, 1000] with a step size of h = 0.5.

Figure 3 .Figure 4 .
Figure 3. Solving the numerical orbit diagram of the Hénon-Heiles system over the interval [0, 17,575] with a step size of h = 0.5.

Figure 5 .
Figure 5. Solving the numerical orbit diagram of the Hénon-Heiles system over the interval [0, 100,000] with a step size of h = 0.5.

Figure 6 .Figure 7 .
Figure 6.Solving the numerical orbit diagram of the Hénon-Heiles system over the interval [0, 10,000] with a step size of h = 0.86.

Table A1 .
A total of 124 preassigned c i for which real solutions b i are solved by order conditions (2).