Adjustment of force-gradient operator in symplectic methods

Many force-gradient explicit symplectic integration algorithms have been designed for the Hamiltonian $H=T (\mathbf{p})+V(\mathbf{q})$ with kinetic energy $T(\mathbf{p})=\mathbf{p}^2/2$ in the existing references. When the force-gradient operator is appropriately adjusted as a new operator, they are still suitable for a class of Hamiltonian problems $H=K(\mathbf{p},\mathbf{q})+V(\mathbf{q})$ with \emph{integrable} part $K(\mathbf{p},\mathbf{q}) = \sum_{i=1}^{n} \sum_{j=1}^{n}a_{ij}p_ip_j+\sum_{i=1}^{n} b_ip_i$, where $a_{ij}=a_{ij}(\textbf{q})$ and $b_i=b_i(\textbf{q})$ are functions of coordinates $\textbf{q}$. The newly adjusted operator is not a force-gradient operator but is similar to the momentum-version operator associated to the potential $V$. The newly extended (or adjusted) algorithms are no longer solvers of the original Hamiltonian, but are solvers of slightly modified Hamiltonians. They are explicit symplectic integrators with symmetry or time-reversibility. Numerical tests show that the standard symplectic integrators without the new operator are generally poorer than the corresponding extended methods with the new operator in computational accuracies and efficiencies. The optimized methods have better accuracies than the corresponding non-optimized counterparts. Among the tested symplectic methods, the two extended optimized seven-stage fourth-order methods of Omelyan, Mryglod and Folk exhibit the best numerical performance. As a result, one of the two optimized algorithms is used to study the orbital dynamical features of a modified H\'{e}non-Heiles system and a spring pendulum. These extended integrators allow for integrations in Hamiltonian problems, such as the spiral structure in self-consistent models of rotating galaxies and the spiral arms in galaxies.


Introduction
In some cases, many trajectories of nonlinear ordinary differential dynamical systems exhibit chaotical behavior; namely, the separations between the trajectories and their nearby trajectories display exponentially sensitive dependence on initial conditions. The trajectories are not analytically integrable, and therefore their studies mainly rely on numerical integration schemes. In general, detecting the chaotical behavior needs long enough numerical integrations with reliable results. Thus, the adopted computational schemes are demanded to perform good stability and high precision. An eighth-and ninth-order Runge-Kutta-Fehlberg integrator [RKF8 (9)] with adaptive step sizes significantly improves the accuracies of the integrals of motion and the trajectories. However, the higher-precision solutions come at the expense of computational time. Furthermore, lower-order numerical methods that preserve geometric properties of the flow of differential equations, i.e., geometric numerical integration schemes [1], are often employed to achieve very accurate long-time determination of the trajectories. The geometric integrators include symplectic integrators for Hamiltonian systems [2][3][4], symmetric integrators for reversible systems, manifold correction schemes for the consistency of integrals (or quasi integrals) of motion [5,6], energy-preserving methods [7], etc.
When dealing with Hamiltonian systems, symplectic integrators are the most appropriate geometric solvers, which own the symplectic nature of Hamiltonian dynamics. The errors in the integrals of motion involving energy integral have no secular growth and tend to zero for infinitesimal time steps. In this sense, symplectic methods approximately conserve the integrals of motion. There are implicit symplectic schemes for inseparable Hamiltonian systems, and explicit symplectic schemes for integrable separable Hamiltonian systems. The implicit schemes were developed by Feng and Qin [8] using generating functions. Implicit symplectic Runge-Kutta methods including the second-order implicit midpoint rule were presented by Sanz-Serna [9]. The family of Gauss-Legndre Runge-Kutta methods [10] are also symplectic. These implicit methods are applied to full Hamiltonian systems. If a part of the Hamiltonian systems is explicitly solved and another part of the Hamiltonian systems is implicitly solved, then explicit and implicit combined symplectic methods [11][12][13] are obtained and can reduce the expense of computational time, compared with the implicit symplectic methods for the full Hamiltonian systems.
Are symplectic integrations always implicit for nonseparable Hamiltonians? No, they are not. In fact, explicit symplectic integrations have been possibly available for some nonseparable Hamiltonians [14]. Relatively recently, explicit symplectic integrators were designed for the Hamiltonian of Schwarzschild spacetime split into four integrable parts [15]. There are a class of explicit extended phase-space symplectic or symplectic-like integrations [16] for an arbitrary separable or inseparable Hamiltonian system. Of course, explicit symplectic integrations are less computationally expensive in general than same order implicit symplectic integrations.
Naturally, explicit symplectic integrations are extensively applied to separable Hamiltonians. Ruth [3] proposed second-and third-order explicit symplectic methods for Hamiltonian systems of the form H = T (p) + V (q). Along this direction, higher order standard explicit symplectic schemes were developed by many authors [17][18][19][20]. It is worth emphasizing that Ruth [3] also gave another third-order symplectic algorithm in which the computation of force gradient is included. Such construction is the force-gradient symplectic integrator. Based on this idea, a series of higher order explicit force-gradient symplectic integration algorithms were established and applied by several authors [21][22][23]. Positive time coefficients can be admissible in many of the algorithmic constructions. In view of this, the force-gradient algorithms with positive intermediate time-steps are suitable for solving time-irreversible problems, such as imaginary time Schrödinger equations. A force-gradient explicit symplectic integration needs less exponential functions of Lie operators than a same order standard explicit symplectic integration. The former algorithm is generally superior to the latter one at same order in accuracy.
When dealing with Hamiltonian systems with the integrable perturbation decomposition form H = H 0 (p, q) + εH 1 (q), the Wisdom-Holman symplectic map of second order [24] drastically improves the numerical accuracy, compared with the standard explicit symplectic method of second order for the Hamiltonian splitting form H = T (p) + V (q). According to the perturbation decomposition of Hamiltonian systems, a number of higher order explicit symplectic integrations were developed and generalized in some references [25,26].
It can be seen from the above presentations that the construction of symplectic integrators is closely related to Hamiltonian systems or their splitting forms. In this contribution, we plan to develop explicit forcegradient symplectic integration algorithms for Hamiltonian systems of the form Because this Hamiltonian is different from the Hamiltonian H = T (p) + V (q), the force-gradient operator for the latter Hamiltonian should be modified appropriately in the former Hamiltonian.
The remainder of the present paper is organized as follows. The force-gradient operator is extended or adjusted in Sect. 2. Taking two models, we estimate the numerical performance of some extended version fourth-order force-gradient symplectic integrators in Sect. 3. For comparison, the standard fourth-order symplectic method of Forest & Ruth [17] and its optimized methods [20] are employed. Possible differences in the dynamics of order and chaos among the tested symplectic integrators are compared, and the dynamics of the two systems is truly described. Finally, the main results are concluded in Sect. 4.

Extension of force-gradient operator
In this Section, some known force-gradient symplectic integrators are introduced. Then, their applications are extended. In what follows, the symplecticity of the extended algorithms is clearly shown.

Existing force-gradient symplectic integrators
A kinetic energy T (p) is a quadratic function of ndimensional momentum vector p = (p 1 , · · · , p n ) and a potential energy V (q) depends on coordinate q = (q 1 , · · · , q n ) only. T and V determine a Hamiltonian Lie derivative operators of T and V are defined as where is a position-version operator. In other words, T has an analytical solution as an explicit function of time.
Because Bp i = −V q i and Bq i = 0, B is a momentumversion operator and is easily analytically solved. The two operators A and B can symmetrically compose a second-order Verlet symplectic integrator [27] where τ is a step size, and W is written in terms of Baker-Campbell-Hausdroff (BCH) formula as The two operators can also symmetrically compose fourth-order explicit symplectic algorithms, such as the Forest-Ruth method [17] M4 = e ατA e β τB e ( 1 2 −α)τA e (1−2β )τB ×e ( 1 2 −α)τA e β τB e ατA , where β = 1/(2 − 3 √ 2) and α = β /2. An optimized Forest-Ruth-like explicit symplectic algorithm of order 4 in [20] is where time coefficients are In fact, only two of the three coefficients sufficiently satisfy the conditions for order 4. In this case, one of the three coefficients has a free choice. However, the free coefficient such as ξ can be determined when the norm of the leading term of fifth-order truncation errors is minimized. This is the so-called optimized method. Another optimized algorithm of order 4 is where time coefficients are Explicit symplectic algorithms to arbitrary even orders were given in [18,19]. Let us consider one of the truncation error terms about τ 3 in Eq.
Thus, C is still a momentum-version operator where V q i q j = ∂ 2 V /∂ q i ∂ q j , and f = ( f 1 , · · · , f n ) with f j = −V q j is a component of the force governed by potential V . That is to say, C is a momentum-version operator with respect to the gradient of the square of force. In this sense, C is called as a force-gradient operator in Refs. [21][22][23].
The above demonstrations show that the third-order truncation error term C in the second-order method M2 and B belong to momentum-version operators. This means that a combination of B and C can yield higherorder algorithms. For example, When A ↔ B in Eq. (13), there is another five-stage fourth-order scheme An optimized five-stage fourth-order method with the operator C [23] is where λ = 1/6, ξ = −17/18000, χ = 71/4500, and k = −4ξ − 1 2 + 3χ 2 . The choice of ξ minimizes the norm of the leading term of fifth-order truncation errors. Two optimized seven-stage fourth-order methods with the operator C in [23] are where where θ = 0.1159953608486416E + 00, λ = 0.2825633404177051E + 00, Note that the combinations of B and C in Eqs. (15)- (17) are slightly unlike those in Ref.
[23], but the time coefficients are the same as those in Ref. [23]. Although the three methods F4O, F4V and F4P are optimized algorithms, their differences are that a free coefficient for F4O minimizes the norm of the leading term of fifth-order truncation errors, but two free coefficients such as ξ and χ for F4V and F4P do. Algorithms (12)- (17) with the force-gradient operator C are the so-called force-gradient explicit symplectic methods in the known publications (e.g., Refs. [21][22][23]). Their constructions are based on one of the third-order truncation error terms in the second-order method M2 included in the potential energy. Only the time coefficient λ (6ξ + k) with respect to the operator C in Eq. (15) is negative, but the time coefficients with respect to the operators A and B in Eqs. (12)- (17) and the time coefficients with respect to the operator C in Eqs. (12)- (14), (16) and (17)

Adjustment of the force-gradient operator
Now, let us suppose a Hamiltonian where a i j = a i j (q) and b i = b i (q). K(p, q) is a formal kinetic energy unlike the kinetic energy in Eq. (1). It is a quadratic function of the momentum vector p, and also depends on the coordinate vector q. Lie derivative operator of K is expressed as where In this sense, A is a momentum-and position-version mixed operator. If K is nonlinear, it is not necessarily analytically solved.
Here, K is assumed to have an analytical solution as an explicit function of time, i.e., operator A is analyti- we obtain the results as follows: where Obviously, D is different from the force-gradient operator C. It is an extended or adjusted operator of the force-gradient operator C.
Algorithms (12)- (17) with the force-gradient operator C become useless for the Hamiltonian J without doubt. However, when the force-gradient operator C is replaced with the extended operator D, these integrators are still available for the Hamiltonian J. That is, they correspond in sequence to the following forms Higher-order force-gradient algorithms in Refs. [22,23] can become the extended methods similar to the constructions (25)- (29). The time coefficients of the composite operators A, B and D in each of the newly extended methods (24)- (29) correspond to those in the existing forcegradient algorithms (12)- (17). No N4O but N2, N4, N4 * , N4V and N4P can have positive intermediate time-steps. These extended methods are symmetric or time-reversible. The condition for symmetry or time-reversibility is, e.g., N4O(τ) × N4O(−τ) = id for the method N4O [1], where id denotes an identical map. The condition is easily checked from a theoretical point of view because the exponents of these operators are linear combinations of τ and τ 3 terms. If the exponents include even power terms of τ such as τ 2 and τ 4 terms, then the condition for symmetry or time-reversibility is not satisfied.

Preservations of symplecticity and volume of the phase space
The standard algorithms M2, M4, M4V and M4P are symplectic because A as a Lie operator with respect to the phase flow of a sub-Hamiltonian H 1 = T (p) is symplectic, and B as a Lie operator with respect to the phase flow of another sub-Hamiltonian H 2 = V (q) is also symplectic. A product of the two symplectic operators and its compositions are still symplectic for the Hamiltonian (2). That is to say, when operators A and B correspond to the phase flows of the two sub-Hamiltonians for the Hamiltonian (2) (namely, the Hamiltonian (2) is a symplectically separable Hamiltonian system), their composition products are naturally symplectic from a physical point of view [28]. Noting this idea, we can similarly show the symplecticity of the existing force-gradient algorithms (12)-(17) and the extended algorithms (24)- (29). Some discussions are given as follows.
At first, we investigate the construction mechanisms of the existing force-gradient algorithms (12)- (17). As is mentioned above, they are constructed by adding one of the third-order truncation error terms in the second-order method M2 to the potential energy. This is one path for understanding the construction mechanisms of the original force-gradient algorithms. Another path is that the force-gradient algorithms are not directly applied to solve the Hamiltonian (18) but are applied to solve some modified Hamiltonians. To show this, we take a Hamiltonian solved by the algorithm F2 as an example. Eq. (12) is a solver of a modified Hamiltonian Operator A acting on the sub-Hamiltonian T advancing a time of τ is the exact phase flow e τA of the sub-Hamiltonian T , and is symplectic. Operator B F2 = (B + τ 2 24 C) acting on the sub-Hamiltonian V F2 advancing a time of τ/2 is the exact phase flow e τ 2 B F2 of the sub-Hamiltonian V F2 , and is also symplectic. Naturally, a composition porduct of these symplectic operators such as e τA e τ 2 B F2 is symplectic [28]. In fact, the explicit symmetric composition scheme F2 for the symplectically separable modified Hamiltonian H F2 is the standard second-order symplectic method M2 with B → B F2 . Of course, different force-gradient symplectic algorithms such as F4 correspond to different modified Hamiltonians. Namely, the force-gradient symplectic algorithms for the original Hamiltonian (2) are the standard symplectic methods for the corresponding modified Hamiltonians.
The construction mechanisms of the extended algorithms (24)- (29) are similar to those of the existing force-gradient algorithms (12)- (17). To show this, we consider the algorithm N2 for the Hamiltonian problem J in Eq. (18). There is a modified potential . Although the function V D is difficulty written in general, it should exist from a mathematical point of view. In fact, we do not need to know what the function V D is, but we need to know only what the function ∂V D /∂ q i is. The modified Hamiltonian solved by the algorithm N2 is a symplectically separable Hamiltonian system J N2 = K + V N2 (q). Only one difference in the construction mechanisms between the forcegradient algorithm F2 and the extended algorithm N2 is that the modified potential V F2 for F2 is easily expressed, whereas the modified potential V N2 for N2 is not. This does not destroy the symplecticity of the algorithm N2 with the extended operator D for the Hamiltonian (18). Similarly, the modified Hamiltonians solved by the algorithms N4 and N4 * can be expressed as symplectically separable Hamiltonian sys- In this way, the symplecticity of the algorithms N4 and N4 * is shown sufficiently. The modified Hamiltonians for the extended algorithms (27)-(29) also exist and therefore the extended algorithms remain symplectic.
Precisely speaking, the symplecticity of each of the aforementioned algorithms should satisfy the condition where S and I are 2n × 2n matrixes. Here, we take four-dimensional phase-space variables Z = (x, y, p x , p y ) as an example to show the expressions of S and I. The at an (m − 1)th step advancing the time step τ for one of the aforementioned algorithms are expressed as where , .
Matrix S for satisfying the condition (30) is a symplectic matrix. Such symplectic matrixes are present for the above algorithms M2, · · · , F2, · · · , N2, · · · . This symplecticity means a symplectic structure described by a closed nondegenerate differential 2-form ω = dx ∧ d p x + dy ∧ d p y on a four-dimensional differential manifold M 4 (symbol ∧ denotes a wedge product). The condition (30) corresponds to the equality of the symplectic structure ω m−1 at the (m − 1)th step and the symplectic structure ω m at the mth step: ω m = ω m−1 . A differential 4-form comes from the product of original differential 2-forms: The preservation of differential 2-form naturally leads to that of differential 4-form where the determinant of the matrix S is det|S| = 1. From a geometric point of view, the differential 4-form corresponds to a volume of the phase space: Therefore, the differential 4-form and volume of the phase space are preserved when det|S| = 1 (see [28] for more details on the symplectic structure and volume of the phase space). The result on det|S| = 1 will be tested in later numerical experiments.
In short, the novel contribution of this paper is to extend the existing force-gradient symplectic algorithms for the Hamiltonian (2) in Refs. [21][22][23] to solve the Hamiltonian (18). Above all, the force-gradient operator C must be replaced by the new operator D.

Numerical simulations
A modified Hénon-Heiles system and a spring pendulum are taken as two models to check the numerical performance of the extended algorithms N4, N4O, N4V and N4P in accuracies of energy and position. Methods M4, M4V and M4P are compared with the extended algorithms. Possible regular and chaotic dynamical differences between the algorithms N4 and M4 are shown. Dynamical behavior of order and chaos in the two problems are described in terms of the best extended algorithm.

Modified Hénon-Heiles system
Let us consider a modified Hénon-Heiles system. Set V as the potential of Hénon-Heiles system [29] V = 1 2 The standard kinetic energy of Hénon-Heiles system is T = p 2 x + p 2 y /2. Here, it is slightly modified as K belongs to one of the forms given in Eq. (19). Obviously, the operator A for K, the operator B for V and the operator D are easily, analytically solvable. Thus, the algorithms N4, N4 * , N4O, N4V, N4P, M4, M4V and M4P are easily available for the system H = K +V (for convenience, J in Eq. (18) is replaced by H). The time step is τ = 0.1. Energy in the system is E = 1/120. Initial conditions are x = 0, y = −2.02 and p y = 0. The initial positive value of p x is determined by E = H. Fig. 1(a) plots energy errors |∆H| = E t − E for several algorithms, where E t denotes the numerical energy at time t. The errors have secular growth for the conventional fourth-order Runge-Kutta (RK4), whereas do not have for the algorithms N4, N4O, N4V, N4P, M4, M4V and M4P. The property without secular drift in the energy errors is due to the symplecticity of these methods. The symplectic methods according to the energy errors listed in Table 1 are mainly divided into three groups as follows. M4 has the poorest anergy accuracy, and N4P and N4V have the best anergy accuracies. N4P and N4V have almost the same accuracy, and their accuracies are about three orders of magnitude better than the accuracy for M4. The four methods N4, M4P, M4V and N4O have minor differences in the energy accuracies. The errors for large to small are N4, M4P, M4V and N4O. In particular, the accuracy for N4 is one order of magnitude better than that for M4. A position error at time t for each of the methods RK4, N4, N4O, N4V, N4P, M4, M4V and M4P is estimated by |∆r| = (x 2 − x 1 ) 2 + (y 2 − y 1 ) 2 , where the solutions (x 1 , y 1 ) are given by the method, and the solutions (x 2 , y 2 ) are obtained from the higherprecision integrator RKF8 (9). When the integration time reaches 10 4 corresponding to 10 5 steps, the position errors are shown in Fig. 1(b) and Table 2. RK4 has the largest error with an order of 10 0.47 . The position error for M4 is also larger than 1. The position errors for N4P and N4V are two orders of magnitude smaller than that for M4, and one order of magnitude smaller than those for the four methods N4, M4P, M4V and N4O. The four methods N4, M4P, M4V and N4O are almost the same in the position errors. These results in Tables 1 and 2 indicate that the standard symplectic integrators without the operator D are poorer than the corresponding extended methods with the operator D in energy accuracies (e.g., M4 inferior to N4, M4V and M4P inferior to N4O, N4V and N4P). The optimized methods are better than the corresponding nonoptimized methods (e.g., M4V and M4P superior to M4, and N4O, N4V and N4P superior to N4). The optimized methods with one free coefficient are inferior to those with two free coefficients (e.g., N4O inferior to N4V and N4P). Clearly, M4 performs the poorest accuracy, and N4V and N4P exhibit the best accura-cies.
The smaller the time step gets, the higher the accuracy of an integrator is. When the time step is appropriately smaller, e.g., τ = 0.01, the energy errors are described in Fig. 1(c) and Table 1. The errors are 10 −6.75 for M4, 10 −8 for N4, M4V, M4P, and N4O, and 10 −9.7 for N4V and N4P. That is, the energy errors typically decrease for each symplectic algorithm. Although RKF8(9) as a non-symplectic method has a secular drift in the energy errors, it has the best energy accuracy and is a good reference integrator for evaluating the performance of the other methods. The position errors in Fig. 1(d) and Table 2 are about 10 −2 for M4 and RK4, 10 −5 ∼ 10 −4 for M4V, M4P, N4 and N4O, and 10 −6 for N4V and N4P. That is to say, as far as the accuracies of energy and position for the smaller time step τ = 0.01 are concerned, M4 is still the poorest one, and N4V and N4P are the best ones among the symplectic methods. Table 3 lists CPU times for each algorithm in Figs. 1 (a) and (c). RK4 is the fastest, and RKF8(9) is the slowest. The computational efficiencies of the symplectic integrators from high to low are N4, N4O, M4, N4P, N4V, M4V and M4P. Of course, there are no relatively dramatic differences in the computational cost between N4 (or N4O) and N4P (or N4V). Fig. 2 is used to check whether the determinants of the matrix S, det|S|, are 1 for the three methods RK4, M2 and N2. The error det|S| − 1 for RK4 grows with time spanning 300. This means that RK4 does not satisfy Eqs. (34) and (35). However, the errors det|S| − 1 for M2 and N2 almost remain at the machine precision in the double-precision environment. This sufficiently supports the preservations of the differential 4-form ω and the phase-space volume in the standard symplectic method M2 and the extended algorithm N2. In principle, det|S| = 1 for RKF8(9) and det|S| = 1 for the fourth-order methods M4, M4V, M4P, N4, N4 * , N4O, N4V and N4P can be confirmed numerically. However, the matrixes S for these algorithms have such long expressions (with over 1000 pages outputted by Matlab) that their determinants are difficultly computed. Fig. 3 shows that the two symplectic integrators N4 and M4 provide different dynamical phase-space structures to the same orbit in Fig. 1 because they have different numerical accuracies for the time step τ = 0.1. A single Kolmogorov-Arnold-Moser (KAM) torus on the Poincaré section x = 0 with p x > 0 is described by the Forest-Ruth method M4 in Fig. 3(a), but many islands are given by the newly extended method N4 in Fig. 3(b). Although the two tori are regular, they are different. In fact, a many-islands torus becomes easier for the occurrence of chaos than a single torus. Which of the algorithms M4 and N4 can provide correct results? N4 can because the results of N4 are consistent with those of higher-precision method RKF8(9) in Fig. 3(c), and are also the same as those of the five methods M4V, M4P, N4O, N4V and N4P. The different KAM tori described by the two methods N4 and M4 are because N4 is one order of magnitude better than M4 in the accuracies of energy and position, as shown in Figs. 1 (a) and (b). In other words, the time step τ = 0.1 is too large to be chosen for M4. However, the phase-space structures can be truly described by M4 for the smaller time step τ = 0.01 because M4 has better accuracies for the smaller time step τ = 0.01 than for the larger time step τ = 0.1 (see Fig. 1 and Tables 1 and 2 for more information).
Only when the initial value y = −1.108 with the initial value p x is altered, does M4 with the larger time step τ = 0.1 indicate that the orbit considered in Fig. 4(a) seems to be a many-islands torus, but exhibits the chaoticity because many discrete points are randomly filled with small regions. Unlike M4, N4 with the larger time step τ = 0.1 and RKF8 (9) show the regularity of the same orbit in Figs. 4 (b) and (c). As claimed above, N4 is superior to M4 in accuracy, therefore, the KAM torus is physically given by N4, but the non-physically spurious chaoticity is caused by M4. These results are also confirmed by fast Lyapunov indicators (FLIs) in Fig. 4(d). The FLIs are from a modified form of Lyapunov exponents [30]. They are originally defined in terms of the lengths of tangent vectors by Froeschlé & Lega [31]. Using the phasespace distances between two adjacent orbits at times 0 and t, d(0) and d(t), Wu et al. [32] suggested the computation of FLI according to the following form A bounded orbit is ordered if its FLI grows algebraically with time log 10 t, but chaotic when its FLI increases exponentially. In other words, the method of complete different growth rates of FLIs with time is faster to distinguish between the two cases of order and chaos than the technique of Lyapunov exponents. When the integration time reaches 3000 in Fig.  4(d), the FLI is 25 for M4, and smaller than 2.5 for N4. Clearly, M4 and N4 for the larger time step τ = 0.1 in-deed give the orbit chaotic and regular dynamical behaviors, respectively. Of course, the chaoticity for M4 is spurious because of M4 performing the poor accuracy. If the smaller time step τ = 0.01 is adopted, M4 like N4 can give the true dynamical behavior to the orbit. Given the initial value y = −1.654 and the larger time step τ = 0.1, the methods of Poincaré section and FLI in Fig. 5 show that the integrated orbit is a regular single torus for M4, whereas a figure-eight orbit for N4 and RKF8 (9). The figure-eight orbit seems to be regular, but is in fact chaotic due to the existence of a hyperbolic point, which has a stable direction and a unstable direction. The chaoticity for N4 is physical, but the regular dynamical information given by M4 is not physical. This is because M4 does not give high enough accuracy to the numerical solutions, but N4 does. M4 can also provide the reliable results for the smaller time step τ = 0.01. In addition, the energy error (not plotted) of the chaotic orbit for N4 in Fig. 5 is similar to that of the regular orbit for N4 in Fig. 1. Namely, the energy accuracy for N4 is independent of the regularity or chaoticity of orbits.
Several main results can be concluded from the above demonstrations. The standard symplectic integrators without the operator D are inferior to the corresponding extended methods with the operator D in computational accuracies and efficiencies. The optimized methods have better accuracies than the corresponding non-optimized methods. N4V and N4P exhibit the best accuracies. Although N4 can provide reliable results on the orbital dynamical behavior for the larger time step τ = 0.1, one of the two optimized extended methods N4P and N4V should be the best integrator from the computational accuracies and efficiencies. Now, we apply the optimized extended method N4P with the large time step τ = 0.1 to explore the dynamics of the modified Hénon-Heiles system. As the magnitude of negative initial value of y decreases, there is a dynamical transition from physical many-islands tori (Figs. 6 (a)-(d)), to single torus (Fig. 6(e)) and to chaotic orbits (Fig. 6(f)). This seems to show that the strength of chaos is enhanced with a decrease of the magnitude of negative initial value of y. However, chaos is not always enhanced, as can be seen from the dependence of FLI on the initial value of y in Fig. 7 that displays the dynamical transition from order to chaos. Chaos mainly occurs for the initial values of y in the vicinity of -2.25∼-2.1, -1.6, and -1.2∼-1. Here, the initial value x = 0 is fixed, and each value of FLI is obtained after integration time t = 3000. Numerical tests show that the threshold of FLIs between the ordered and chaotic cases is 4. The FLIs larger than the threshold determine the chaoticity of bounded orbits, but the FLIs less than the threshold indicate the regularity of bounded orbits.

Spring pendulum
A spring pendulum in polar coordinates is described by the Hamiltonian [1] The Hamiltonian is divided into kinetic energy K and potential energy V as follows: K is one of the expressions in Eq. (19) and is analytically solved. In this case, the extended algorithms such as N4 can be suitable for integrating the spring pendulum problem. Taking the step size τ = 0.1 and energy E = 1/12, we choose initial conditions r = 1.15, p r = 0 and ϕ = 0.05π. The initial value p ϕ > 0 is solved from the energy equation E = H. Fig. 8 and Table 4 show that RK4 has the largest errors in the energy and position, and RKF8(9) exhibits the smallest energy error. The errors in the energy have secular drifts for the non-symplectic methods RK4 and RKF8(9), but do not have for the symplectic methods M4, M4V, M4P, N4, N4O, N4V and N4P. The errors in the energy and position for N4V and N4P are about three orders of magnitude smaller than those for M4. Table  5 lists CPU times of the methods. The standard symplectic methods are slower than the corresponding extended schemes in computational efficiencies. The optimized methods N4V and N4P need small additional cost compared with the corresponding non-optimized method N4.
The tested orbit in Fig. 8 is a regular single closed torus on the Poincaré section ϕ = 0 and p ϕ > 0 in Fig 9(a). Unlike in the modified Hénon-Heiles system, M4 in the present problem is the same as anyone of the methods M4V, M4P, N4, N4O, N4V, N4P and RKF8 (9) in the description of phase-space structures. This is because the accuracies in the energy and position for M4 in Fig. 8 are higher than those in Fig.  1. In fact, the energy for M4 is accurate to an order of 10 −3 , and the position for M4 is accurate to an order of 1 in Fig. 1 when the integration time t = 10 4 . However, the energy for M4 is accurate to an order of 10 −5 , and the position for M4 is accurate to an order of 0.1 in Fig. 8 when the integration time t = 10 4 . Because of this, M4 can truly describe the phase-space structures, as the methods M4V, M4P, N4, N4O, N4V, N4P and RKF8 (9) can. For given integrator and time step, the numerical accuracy closely depends on model Hamiltonians, and particularly depends on the periods of orbits in the Hamiltonians. The larger the periods are, the better the numerical accuracy is.
Considering the best performance in numerical accuracies and computational efficiency, we use N4P to trace different phase-space structures when various initial values are given to ϕ. For example, single-torus orbits exist for ϕ = 0.19π in Fig. 9(b) and ϕ = 0.358π in Fig. 9(e). There are many islands for ϕ = 0.35π in Fig. 9(d), ϕ = 0.361π in Fig. 9(f), ϕ = 0.366π in Fig. 9(g), and ϕ = 0.39π in Fig. 9(i). Chaos occurs for ϕ = 0.2π in Fig. 9(c) and 0.376π in Fig. 9(h). No universal rule for the dynamical transition from order to chaos seems to be given to a variation of the initial value ϕ. However, Fig. 10(a) for the description of the initial values of ϕ and their corresponding FLIs shows that chaos mainly occurs for the initial values of ϕ in the vicinity of 0.25, 1.75, and 2.25. This is because the spring pendulum suffers from strong perturbations in the vicinity of ϕ = π/4, 7π/4. Scanning the initial values of r and their corresponding FLIs in Fig. 10(b) displays that chaos mainly occurs for the initial values of r in the vicinity of 0.5, and 0.75∼2.25.
Besides the method of FLIs, the 0-1 test for chaos [33] can be applied to explore the transition from order to chaos with the initial value ϕ or r varying. The 0-1 test chaos indicator is described here. Set Z = (r, ϕ, p r , p ϕ ) as a solution of the system (39) at time t. ψ(Z) is a function of Z; for example, ψ(Z) = r is a simple choice. The authors of [33] defined two functions where c > 0 is an arbitrarily constant. The mean-square displacement of q(t) is Finally, the asymptotic growth rate of the mean-square displacement is expressed as Based on ergodic theory, Λ = 0 signifies regular dynamics, but Λ = 1 signifies chaotic dynamics. Taking c = 1.8, T = 1000000, t = 1000 and time step τ = 0.1, we calculate the values of Λ with respect to the orbits in Fig. 10 and obtain the dependence of Λ on the initial values ϕ in Fig. 11. The values of Λ in Fig. 11 That is to say, the dynamical properties described by the 0-1 indicator Λ in Fig. 11 are consistent with those given by the FLIs in Fig. 10. However, because T is large enough, computations of L(t) are relatively expensive. In fact, CPU time is 284.41 seconds for each of the Λ values in Fig.  11, and 0.36 seconds for each of the FLIs in Fig. 10. Thus, the FLIs are quicker to distinguish between the ordered and chaotic two cases than the 0-1 test indicator.

Conclusions and discussions
Many force-gradient explicit symplectic integration algorithms with the force-gradient operator C for the Hamiltonian (2) with the kinetic energy (1) have been in Refs. [21][22][23]. However, these algorithms become useless for the Hamiltonian (18) with the integrable kinetic energy (19) if the force-gradient operator C is not altered. We find that the existing integrators are still available for the Hamiltonian (18) when the forcegradient operator C gives place to a new operator D. This new operator is not a force-gradient operator but is similar to the momentum-version operator associated to the potential V . The extended algorithms are no longer solvers of the original Hamiltonian but are solvers of slightly modified Hamiltonians. They are explicit symplectic integrators with symmetry or timereversibility.
Numerical tests show that the standard symplectic integrators without the operator D are generally inferior to the corresponding extended methods with the operator D in computational accuracies and efficiencies. For example, the fourth-order Forest-Ruth symplectic scheme cannot provide reliable results to the description of regular and chaotic dynamical features of the modified Hénon-Heiles system for the use of some appropriately large time steps, but the corresponding extended ones can. The optimized methods have better accuracies than the corresponding non-optimized methods. Among the tested symplectic methods, the two extended optimized sevenstage fourth-order methods of Omelyan, Mryglod and Folk (N4V and N4P) exhibit the best numerical performance, and their accuracies are about three orders of magnitude better than the accuracies of the Forest-Ruth symplectic scheme. Finally, one of the two optimized algorithms is used to study the orbital dynamical features of the modified Hénon-Heiles system and the spring pendulum.
The proposed extended algorithms are suitably applicable to any Hamiltonian systems like the Hamiltonian system (18). The kinetic energies like Eq. (19) can be found in some references. For example, the Hamiltonian in Eq. (5.21) of Ref. [34] is The kinetic energies given in Eq. (40) are universal for two-dimensional Hamiltonian systems H = m(ẋ 2 + y 2 )/2 + V(x, y) in polar coordinates (r, ϕ), where potentials V (x, y) are non-axisymmetric. In addition, three-dimensional Hamiltonian problems H = m(ẋ 2 + y 2 +ż 2 )/2 +V (x, y, z) in spherical coordinates (r, θ , ϕ) are expressed as which resembles the Hamiltonian (18) with the integrable kinetic energy (19). In other words, the Hamiltonian (18) is not restricted to several special examples, but is often met in many situations. In this sense, this extension has wide applications. As an example, the kinetic energies in the Hamiltonians of the spiral structure in self-consistent models of rotating galaxies [35] and the spiral arms in galaxies [36] are Eq.
(40). Thus, the newly extended force-gradient explicit symplectic methods should be suitable for the study of chaotic spiral arms in the models of rotating galaxies. This problem will be considered in a future work.        Tables 1 and 2.     Fig. 1. M4 gives a single torus to the orbit, but N4 like RKF8(9) exhibits many islands. The result is reliable for N4, but is not for M4 because the larger time step τ = 0.1 causes M4 to yield poor accuracy to the numerical solutions. The orbit for M4 in (a) seems to be a chaotic many-islands torus, but the orbits for N4 and RKF8(9) in (b) and (c) are regular many-islands tori only.
(d) Fast Lyapunov indicators (FLIs) for the two methods M4 and N4 integrating the orbit. The FLIs given by RKF8 (9) are almost the same as those given by N4. When the integration ends, the FLI is smaller than 2.5 for N4, but reaches 25 for M4. The FLI for N4 indicates the regularity, while the FLI for M4 shows the chaoticity. The results are correct for N4, whereas chaos for M4 is spurious in panels (a) and (d) due to the large time step τ = 0.1 leading to M4 with poor accuracy to the numerical solutions.     Table 4.  Fig. 7 is a regular single-torus orbit in Fig. 8(a). Regular single-torus orbits also appear in (b) and (e). The orbits in (d), (f), (g) and (i) are many-islands tori, but are chaotic in (c) and (h). These results are consistently given by the methods M4, M4P, M4V, N4, N4O, N4P, N4V and RKF8 (9).   Fig. 11.-Same as Fig. 10 but the dependence of the 0-1 test chaos indicator Λ on the initial value ϕ or r. The regular and chaotic properties described by the 0-1 test are the same as those given by the FLIs in Fig. 10.