© The Authors 2025

1 Introduction

Symplecticity is characteristic of Hamiltonian systems and numerical results obtained by symplectic integrators maintain long-term stability compared to nonsymplectic integrators. Thus, symplectic integrators have proven to be powerful tools in studies of gravitational dynamics (Vogelaere 1900; Ruth 1983; Feng 1985; Wang et al. 2021).

An important approach to constructing symplectic integrators is based on the use of generating functions (Channell & Scovel 1990; Miesbach & Pesch 1992; Chen et al. 2003; Makazaga & Murua 2009). Gao et al. (2012) constructed a high-order and fixed-step symplectic integrator SQQ for Hamiltonian systems based on the principle of least action and generating function. In the acronym SQQ, S denotes the symplectic nature of the integrator, while the two Qs indicate that generalized coordinates are used as independent variables at both sides in the generating function (Gao et al. 2012). Overall, SQQ is a general-purpose integrator, and has been applied to trajectory planning (Peng et al. 2023), kinodynamic planning, optimal control of cable-driven tensegrity manipulators (Song et al. 2020), and multibody dynamics problems (Peng et al. 2019, 2020; Li et al. 2021). The numerical results demonstrate that SQQ exhibits high accuracy and effective energy conservation behavior.

However, in the numerical computation of multi-body gravitational problems, traditional fixed-step symplectic algorithms face two primary challenges. The first challenge arises from the computational efficiency issues caused by the multi-scale nature of time in such systems. For systems with strong variations in interaction, such as high-eccentricity or hyperbolic Kepler problems, the fixed-step symplectic integrator becomes inefficient. In contrast, adaptive time-step integrators can handle these scenarios effectively by employing smaller time steps in rapidly changing regions and larger time steps in slowly varying regions (Dehnen & Read 2011; Springel 2005). Unfortunately, traditional step-size control strategies destroy the long-term conservation properties of symplectic integrators (Preto & Tremaine 1999; Emel’yanenko 2007; Blanes & Iserles 2012).

One approach to constructing adaptive time-step symplectic integrators involves the application of time transformation, which decouples the time step from the integration step (Hairer 1997; Mikkola & Tanikawa 1999; Preto & Tremaine 1999; Wang et al. 2020; Wang & Nitadori 2020). However, time transformation often leads to a non-separable Hamiltonian system, which, in turn, transforms an explicit algorithm into an implicit one (Hairer 1997; Huang & Leimkuhler 1997). In general, implicit algorithms require solutions to nonlinear equations. Due to the complexity introduced by the time-transformed Hamiltonian, computing the Jacobian matrix for solving this nonlinear system can become computationally expensive. This leads to the second challenge: how to efficiently solve the nonlinear equations at each time step. To address this challenge, we employed a quasi-Newton method based on Broyden’s method (Broyden 1965) to accelerate the inversion of the Jacobian matrix, which has been widely used in many research fields (e.g. Bulín & Hajžman 2021; Chen et al. 2022; Sherman & Morrison 1950). Broyden’s method iteratively updates an approximation to the Jacobian inverse during the process of solving the nonlinear equations. By avoiding the explicit computation of the Jacobian matrix, this approach significantly reduces computational costs and improves efficiency.

To address the two important challenges mentioned above, our proposed method combines the time transformation and the quasi-Newton approach, offering a robust framework for efficient and accurate simulations in multi-body gravitational problems. This paper is organized as follows. Section 2 provides the fundamental theory of SQQ based on Chebyshev interpolation method. Section 3 presents the numerical implementation, which includes three parts. The first part introduces a projection method that accelerates the computation of interpolation functions and their derivatives without requiring matrix inversion. The second part introduces the concept of the time transformation and discusses the selection of the step-size control function. The last part presents an efficient approach for solving nonlinear equations. Finally, Section 4 includes several numerical examples to demonstrate the effectiveness of the proposed method.

2 The symplectic algorithm based on Generating Function

This section gives a brief introduction to SQQ (for more details, readers are referred to Gao et al. 2012). The action of a conservative system is the integral of the Lagrangian function over time $$S={\int }_{t_a}^{t_b}[{\mathbf{p}}^T\dot{\mathbf{q}}-H(\mathbf{q},\mathbf{p})]dt,$$(1)

where q is the generalized displacement, $\dot{\mathbf{q}}$ is the generalized velocity, and p is the generalized momentum. The principle of least action is given by $$\begin{array}{rl}\delta S=& {\int }_{t_a}^{t_b}[(\delta \dot{\mathbf{q}}{)}^T\mathbf{p}-(\delta \mathbf{q}{)}^T\frac{\mathrm{\partial }H}{\mathrm{\partial }\mathbf{q}}]dt+{\int }_{t_a}^{t_b}(\delta \mathbf{p}{)}^T(\dot{\mathbf{q}}-\frac{\mathrm{\partial }H}{\mathrm{\partial }\mathbf{p}})dt\\ =& {\int }_{t_a}^{t_b}(\delta \mathbf{q}{)}^T[-\dot{\mathbf{p}}-\frac{\mathrm{\partial }H}{\mathrm{\partial }\mathbf{q}}]dt+{\int }_{t_a}^{t_b}(\delta \mathbf{p}{)}^T(\dot{\mathbf{q}}-\frac{\mathrm{\partial }H}{\mathrm{\partial }\mathbf{p}})dt\\ & +{{\mathbf{p}}^T\delta \mathbf{q}|}_{t_a}^{t_b}=0.\end{array}$$(2)

Equation (2) shows that if the Hamiltonian canonical differential equation is satisfied within the time step, the action will satisfy the following equation, $$\delta S={\mathbf{p}}^T\left(t_b\right)\mathrm{d}\mathbf{q}\left(t_b\right)-{\mathbf{p}}^T\left(t_a\right)\mathrm{d}\mathbf{q}\left(t_a\right).$$(3)

The physical meaning of Eq. (3) is that the action, S, depends only on the generalized displacement at both ends of the time interval [t_a, t_b]. In particular, S is a function of q (t_a) and q (t_b). From Eq. (3), we get $$\mathbf{p}\left(t_a\right)=-\frac{\mathrm{\partial }S(\mathbf{q}\left(t_a\right),\mathbf{q}\left(t_b\right))}{\mathrm{\partial }\mathbf{q}\left(t_a\right)},\mathbf{p}\left(t_b\right)=\frac{\mathrm{\partial }S(\mathbf{q}\left(t_a\right),\mathbf{q}\left(t_b\right))}{\mathrm{\partial }\mathbf{q}\left(t_b\right)}.$$(4)

Once the generating function is determined, the generalized displacement q (t_b) and the generalized momentum p(t_b) can be calculated according to Eq. (4). In this context, SQQ employs the Lagrange interpolation with equidistant nodes and then the Gaussian quadrature is used to calculate the approximate value of S. However, a higher order interpolation with equidistant nodes often leads to the Runge phenomenon. To avoid this phenomenon, we used the Chebyshev interpolation in this work (Cai & Song 2012; Lanczos 1938).

The generalized coordinate and momentum vectors as functions of time, t, are expressed as: $$\mathbf{q}(t)=\sum _{k=0}^m{M}_k(t){\mathbf{q}}_k,\mathbf{p}(t)=\sum _{k=0}^n{N}_k(t){\mathbf{p}}_k,$$(5)

where M_k(t) and N_k(t) are the interpolation basis functions, while q_k and p_k are defined as: $${\mathbf{q}}_k={\{q_k^1,q_k^2,\cdots ,q_k^d\}}^{\mathrm{T}},k=0,1,\cdots ,m,$$(6) $${\mathbf{p}}_k={\{p_k^1,p_k^2,\cdots ,p_k^d\}}^{\mathrm{T}},k=0,1,\cdots ,n,$$(7)

with d representing the degrees of freedom (DOF) of the system. Here, m and n denote the number of interpolation points for the generalized displacement and generalized momentum, respectively.

By substituting Eq. (5) into Eq. (1) and then combining it with Eq. (4), we obtain $$\begin{array}{rl}& \frac{\mathrm{\partial }S}{\mathrm{\partial }{\mathbf{q}}_0}+{\mathbf{p}}_0=0,\\ & \frac{\mathrm{\partial }S}{\mathrm{\partial }{\mathbf{q}}_i}=0,i=1,2,\cdots ,m-1,\\ & \frac{\mathrm{\partial }S}{\mathrm{\partial }{\mathbf{p}}_i}=0,i=0,1,\cdots ,n.\end{array}$$(8)

and $${\mathbf{p}}_n=\frac{\mathrm{\partial }S}{\mathrm{\partial }{\mathbf{q}}_m}.$$(9)

Generally, Eq. (8) is composed of nonlinear algebraic equations. Due to the introduction of the step-size control function in this paper (see Section 3.2), the analytical derivation and calculation of the Jacobian matrix become complicated and challenging. This paper employs a quasi-Newton method and accelerates the iterative computation using Broyden’s method (see Section 3.3).

3 Numerical implementation

3.1 Efficient computation of interpolation functions using the projection method

To solve Eq. (8), the action S was approximated using Gaussian integration, which requires the calculation of Gaussian points and their corresponding interpolation functions, M_k(t) and N_k(t), at each time step. Since this process involves matrix inversion, it is computationally expensive. However, because the number of interpolation points remains constant during the computation, these functions do not need to be recomputed at each step. Instead, they can be efficiently determined using a projection method. This section will illustrate this method by using N_k(t) and ${\dot{N}}_k(t)$ as examples.

Thus, we let Gaussian integral points and their weights in the interval [−1, 1] be $${\hat{\xi }}_j,{\hat{\omega }}_j,j=1,2,\cdots ,g,$$(10)

where g represents the number of Gaussian quadrature points, and $\hat{N}$ and $\hat{\dot{N}}$ are the interpolation function and its derivative, respectively. In any time interval [t_a, t_b] (0 ≤ t_a < t_b), the Gaussian integral points and their weights can be obtained as follows: $${\xi }_j=t_a+\frac{{\hat{\xi }}_j+1}{2}L,{\omega }_j=\frac{L}{2}{\hat{\omega }}_j,j=1,2,\cdots ,g,$$(11)

where L = t_b − t_a. As illustrated in Fig. 1, $\hat{\xi }$ and $\xi (\xi \ne \hat{\xi })$ denote the Gaussian quadrature points on the intervals [−1, 1] and [t_a, t_b], respectively. These points share equivalent “length coordinates” denoted by λ₁ and λ₂, which represent the proportion of the interval’s total length occupied by a given coordinate point. Consequently, the interpolation function within the interval [t_a, t_b] can be expressed as $$\mathit{N}\left({\xi }_j\right)=\hat{\mathit{N}}\left({\hat{\xi }}_j\right),j=1,2,\cdots ,g,$$(12)

Taking the derivative of both sides of the Eq. (12), we obtain $$\dot{\mathit{N}}\left({\xi }_j\right)=\frac{2}{L}\widehat{\dot{\mathit{N}}}\left({\hat{\xi }}_j\right),j=1,2,\cdots ,g.$$(13)

Clearly, it is only necessary to compute $\hat{N}$ and $\hat{\dot{N}}$, and then use the Eqs. (12) and (13) to obtain N and $\dot{\mathit{N}}$ in any time interval. The variables $\hat{\mathit{N}}$ and $\widehat{\dot{\mathit{N}}}$ can be calculated using the following equations $$\hat{\mathit{N}}(t)={\mathit{G}}^{-1}\mathit{p}(t),\dot{\mathit{N}}(t)={\mathit{G}}^{-1}\dot{\mathit{p}}(t),$$(14)

where $$\mathit{G}=\left[\begin{array}{lllll}1& t_1& t_1^2& \cdots & t_1^{n-1}\\ 1& t_2& t_2^2& \cdots & t_2^{n-1}\\ ⋮& \cdots & \cdots & \cdots & ⋮\\ ⋮& \cdots & \cdots & \cdots & ⋮\\ 1& t_n& t_n^2& \cdots & t_n^{n-1}\end{array}\right],\mathit{p}=\left[\begin{array}{l}1\\ t\\ t^2\\ ⋮\\ t^n\end{array}\right],\dot{\mathit{p}}=\left[\begin{array}{l}0\\ 1\\ 2t\\ ⋮\\ n{t}^{n-1}\end{array}\right]$$(15)

and t₁, t₂, . . . , t_n are the Chebyshev interpolation points on the interval [−1, 1].

In this paper, a warm-up procedure is used to compute $\hat{N}(\hat{\xi })$ and $\hat{\dot{N}}(\hat{\xi })$ in advance. Vectors N (ξ) and $\dot{\mathit{N}}(\xi )$ are calculated at each step by using Eqs. (12) and (13), which reduces the computational costs of the matrix inversion. Similarly, this approach is also applied to compute M(ξ) and $\dot{\mathit{M}}(\xi )$.

Fig. 1 Illustration of the length coordinates for the time intervals [−1, 1] and [ta, tb]. Although the horizontal coordinates of $\hat{\xi }$ and ξ are not equal and their length coordinates λ1 and λ2 are equal, the values of the corresponding interpolation functions are equal (Eq. (12)).

3.2 Time transformation

The time transformation t ↔ τ can be employed to implement variable step sizes. Integrating the original system with a variable step size Δt is equivalent to integrating the transformed system with a constant step size, Δτ. The relationship between the transformed time, τ, and the original time, t, is given by $$\frac{\mathrm{d}t}{\mathrm{d}\tau }=\sigma (\mathbf{p},\mathbf{q}),$$(16)

where σ > 0 is the step size control function. To ensure that the transformed system remains Hamiltonian, the Darboux–Sundman transformation (Nacozy 1977) is expressed as follows: $$K=\sigma (H-H_0).$$(17)

Here, H is the Hamiltonian of the original system, K is the new Hamiltonian system, and H₀ = H(p₀, q₀) is the initial energy of the system. For the N-body problem, the Hamiltonian system is generally separable. However, introducing the time transformation makes the system nonseparable. As a result, explicit methods become implicit, requiring the solution of nonlinear equations. Furthermore, due to the time transformation, the Jacobian matrix of the new Hamiltonian system becomes more complicated, which significantly increases the computational cost if solving nonlinear equations using Newton’s method. An efficient quasi-Newton method based on Broyden’s approach for solving nonlinear equations is presented in Section 3.3.

In some cases, appropriate step size control functions σ(p, q) are known a priori. For example, Budd & Piggott (2003) proposed an effective step size control function for the two-body problem: $${\sigma }_1(\mathbf{p},\mathbf{q})=\parallel \mathbf{q}{\parallel }^{\alpha },$$(18)

where α = 2 or 1.5. Smaller steps are taken when the bodies are close to each other. In the N-body problem, close encounters can lead to a significant increase in velocity, when the total energy of the system is primarily dominated by kinetic energy and a smaller time step should be chosen in such cases. Therefore, it is natural to use the following step size control function (Huang & Leimkuhler 1997; Hairer et al. 2006): $${\sigma }_2(\mathbf{q})={[(H_0-U(\mathbf{q}))+\mathrm{\nabla }U(\mathbf{q}{)}^T{M}^{-1}\mathrm{\nabla }U(\mathbf{q})]}^{-1/2},$$(19)

where U(q) is the potential energy. Equation (19) ensures that the step size adapts naturally to the dynamic properties of the system. To prevent the step size from becoming excessively small or large, and upper and lower bounds are introduced (Huang & Leimkuhler 1997). The step size function used in this paper is defined as: $$\sigma (p,q)=\frac{\sqrt{{\sigma }_2^2+a^2}}{\frac{1}{b}\sqrt{{\sigma }_2^2+a^2}+1},$$(20)

where 0 < a ≪ 1 ≪ b. This formulation ensures that σ(p, q) is constrained within the range defined by the lower bound a and the upper bound b. In this paper, σ(p, q) is adopted as the time step function with a = 10⁻⁶ and b = 10².

3.3 An efficient quasi-Newton method based on Broyden’s method

We let F : R^d → R^d be a vector valued function, defined by Eq. (8). Then, the Newton method for solving the nonlinear equations F(x) = 0 can be expressed as follows: $${\mathbf{x}}_{k+1}={\mathbf{x}}_k-{\mathbf{J}}_k^{-1}\mathbf{F}\left({\mathbf{x}}_k\right),$$(21)

where k is the iteration index. In each iteration of solving the nonlinear equations, the Jacobian matrix, J_k, and its inverse, ${\mathbf{J}}_k^{-1}$, are typically required. However, deriving the Jacobian analytically becomes highly complex due to the time transformation applied to the Hamiltonian. A common method to approximate the Jacobian is numerical differentiation. While effective, this approach can be computationally expensive, especially when evaluating the nonlinear function, F, involves Gaussian quadrature.

To accelerate the computation, we employ the quasi-Newton method instead of standard Newton method. Specifically, we used Broyden’s method, which iteratively updates the inverse of the Jacobian matrix without explicitly computing the Jacobian itself. This allowed us to significantly reduce the computational overhead while maintaining reasonable convergence properties.

Following by Broyden’s method, the inverse Jacobian at iteration k + 1 was updated using information from the previous iteration, k, according to the following update formula: $${\mathbf{J}}_{k+1}^{-1}={\mathbf{J}}_k^{-1}+\frac{({\mathbf{s}}_k-{\mathbf{J}}_k^{-1}{\mathbf{y}}_k){\mathbf{s}}_{\mathrm{k}}^{\mathrm{\top }}{\mathbf{J}}_k^{-1}}{{\mathbf{s}}_k^{\mathrm{\top }}{\mathbf{J}}_k^{-1}{\mathbf{y}}_k},$$(22)

where s_k = x_k − x_k−1 and y_k = F(x_k) − F(x_k−1). This formula provides an efficient method for updating ${\mathbf{J}}_k^{-1}$, without requiring a direct recalculation of the Jacobian matrix. This is especially beneficial in cases where F(x) involves computationally intensive procedures, such as Gaussian integration. The pseudocode used in this study for the quasi-Newton method based on the Broyden method, is given in Algorithm 1, presented below.

Throughout the run of Algorithm 1, the solution from the previous time step x_k is utilized as the initial guess x₀ for the subsequent time step. Similarly, the inverse Jacobian, ${\mathbf{J}}_0^{-1}$, is initialized by the inverse Jacobian, ${\mathbf{J}}_k^{-1}$, from the prior converged time step. For the first time step, the Jacobian J₀ is approximated numerically using the difference quotient method.

Nonlinear equations can also be solved by using other solvers such as the tensor method. The tensor method introduces an additional tensor term in the iteration to capture the higher-order nonlinear behavior (Schnabel & Frank 1984). Compared to quasi-Newton methods (e.g. Broyden’s method, as used in this work), the tensor method comes with higher computational costs (Nocedal & Wright 2006) and is generally more difficult to implement due to the introduction of the additional tensor term.

However, one notable problem of the quasi-Newton method is its potentially slower convergence compared to the standard Newton iteration (Broyden 1970). This issue can often be mitigated by employing smaller step sizes, which improve the stability and accuracy of the iterations while maintaining the efficiency of the quasi-Newton approach (Bulín & Hajžman 2021; Nocedal & Wright 2006).

Fig. 2 Evolution of the energy error for SQQ-P and SQQ. Here, e = 0.5, and m and n denote the number of interpolation points for the generalized displacement and generalized momentum, respectively.

4 Numerical example

To evaluate the performance of the proposed integration methods, we present three numerical examples. Before delving into the examples, we summarize the key characteristics of the integrators involved.

The original SQQ uses equidistant interpolation nodes and solves nonlinear equations with the Newton method. Five extensions of SQQ are evaluated in this Section: SQQ-P, SQQ-PN, SQQ-PQ, SQQ-PTN, and SQQ-PTQ. The suffix “T” denotes the use of the time transformation, as defined in Eq. (20), “P” refers to the projection method for efficient interpolation, while “Q” and “N” indicate the use of quasi-Newton and Newton methods for solving nonlinear equations, respectively. The characteristics of these methods are summarized in Table 1. For consistency, the convergence threshold, ε, for solving nonlinear equations is set to 1 × 10⁻¹² for all examples. All numerical computations were performed in MATLAB on a platform equipped with a 13th Generation Intel(R) Core(TM) i7-13700 processor and 32 GB of RAM.

4.1 The Kepler problem

The Kepler problem with the Hamiltonian is given by $$H=\frac{1}{2}(p_1^2+p_2^2)-\frac{1}{r},r=\sqrt{q_1^2+q_2^2},$$(23)

where the initial conditions are $$\mathbf{q}=\left[\begin{array}{c}1-e\\ 0\end{array}\right],\mathbf{p}=\left[\begin{array}{c}0\\ \sqrt{\frac{1+e}{1-e}}\end{array}\right],$$(24)

and e represents the eccentricity of the orbit. The period of this system is 2π.

To evaluate the efficiency of the projection method, e = 0.5 is selected as the test parameter. Simulations of 500 orbital periods are performed using both the SQQ and SQQ-P.

Figure 2 illustrates the evolutions of the energy error for SQQ and SQQ-P. When the number of interpolation points is relatively small (i.e., m < 15 and n < 15), the Runge phenomenon does not occur. In this example, the evolutions of the energy error for SQQ and SQQ-P are nearly identical, exhibiting oscillatory behavior within a bounded range, which is typical for symplectic methods. This demonstrates that the projection method employed in this study provides an accurate computation of the interpolation function. As the number of interpolation points increases (i.e., m = 18 and n = 18), SQQ exhibits a growth in energy error over time. This behavior arises from the excessive number of interpolation points, which triggers the Runge phenomenon. As a result, the interpolated action, S, may become inaccurate. In contrast, SQQ-P, which uses Chebyshev interpolation points, effectively suppresses the Runge phenomenon. As a result, the energy error for SQQ-P continues to oscillate within a bounded range, as shown in Fig. 2. This highlights the robustness of Chebyshev interpolation in maintaining the stability of the computed action.

Table 2 shows the time step sizes and CPU times for the various simulations. As the number of interpolation points increases, the CPU time ratio decreases (the CPU time ratio is defined as the ratio of SQQ-P’s CPU time to that of SQQ). This happens because the increase in the number of interpolation points introduces additional Gaussian nodes, significantly raising the computational cost of evaluating the interpolation function during each time step. In contrast, SQQ-P employs the projection method proposed in this study, thereby enhancing the efficiency.

The Kepler problem becomes challenging for traditional fixed-step algorithms as the eccentricity approaches to 1. We perform two numerical experiments for eccentricities of e = 0.9 and 0.99, with parameters Δτ = 0.01, m = n = 3, and a total integration time of 500 orbital periods. Fig. 3 illustrates the evolution of the energy error for SQQ-PTQ with these large eccentricities. Fig. 4 shows the corresponding phase space trajectories of q and p. It can be seen that, even over long-term integration, the phase-space structure remains nearly intact, which indicates that SQQ-PTQ effectively preserves the Hamiltonian structure.

Table 3 presents the CPU times of SQQ-PTN and SQQ-PTQ for solving the Kepler problem with high eccentricities. SQQ-PTQ exhibits an improvement in efficiency over SQQ-PTN, with the latter requiring about 2–3 times of the CPU time at e = 0.9 and e = 0.99, respectively.

Fig. 3 Evolution of the energy error for SQQ-PTQ for the Kepler problem with high eccentricities.

4.2 Three-body problem

For problems involving close encounters, any unreasonable time step will lead to increased errors or calculation failure, so an appropriate time step is necessary to accurately capture such events. The three-body problem with possible close encounters is a good option for assessing the performance of the adaptive time step strategy of integrators. Recently, Liao et al. (2022) employed machine learning techniques to discover some novel periodic solutions for the three-body problem. In this section, we use one specific periodic configuration from their study to test SQQ-PTQ’s ability to automatically select the time step. The orbital and physical parameters of the three bodies, as given in Liao et al. (2022), are as follows: $$\begin{array}{lll}{\mathbf{q}}_1={\left[\begin{array}{l}-0.2227,0\end{array}\right]}^{\mathrm{T}},& {\dot{\mathbf{q}}}_1=[0,1.7813{]}^{\mathrm{T}},& m_1=0.9,\\ {\mathbf{q}}_2=[\begin{array}{l}1,0{]}^{\mathrm{T}}\end{array},& {\dot{\mathbf{q}}}_2=[0,0.4150{]}^{\mathrm{T}},& m_2=0.85,\\ {\mathbf{q}}_3=[\begin{array}{l}0,0{]}^{\mathrm{T}}\end{array},& {\dot{\mathbf{q}}}_3=[0,-1.9559{]}^{\mathrm{T}},& m_3=1.\end{array}$$(25)

The period of the system is T = 6.3509, and the gravitational constant G is normalized to 1. We simulated 500 orbital periods of this system using the SQQ-PTQ. Figure 5 shows the trajectories of three particles, where particles 1 and 3 experience two close encounters within one period. These encounters occur at approximately 0.4T (dashed box a), and at 0.6T (dashed box b). For both encounters, the closest distance between particles 1 and 3 is approximately 0.014. Figure 6 shows the magnitude of the velocities of the three particles (during the first period). The velocities of particles 1 and 3 exhibit two peaks within one period, due to their close encounters, leading to significant velocity changes for both particles. These two close encounters cause the maximum velocities of particles 1 and 3 to be 10 and 34 times of their respective minimum velocities. The evolution of the energy error exhibits stable oscillations as shown in Fig. 7. Furthermore, the phase space trajectories shown in Fig. 8 maintain a well-preserved structure throughout long-term integration, indicating that SQQ-PTQ effectively preserves the symplectic property.

We also use the MATLAB built-in general-purpose integrator ODE45 (with relative error tolerance set to 1 × 10⁻¹⁰) and SQQ-PTN to simulate this system for 500 orbital periods. Table 4 lists the computational time and maximum energy errors for the three integrators. The only difference between SQQ-PTQ and SQQ-PTN lies in the method used to solve the nonlinear equations, which explains why their maximum energy errors are nearly identical. The CPU time of SQQ-PTQ is less than one-third of that of SQQ-PTN. This indicates that the quasi-Newton method accelerates the computation without introducing any additional errors.

Fig. 4 Phase space trajectories over 500 orbital periods for the Kepler problem with different eccentricities.

Fig. 5 Trajectories of particles in the three-body problem, with dashed boxes a and b indicating the locations of two close encounters. Box a corresponds to a close encounter occurring at approximately 0.4T, and box b corresponds to the second close encounter occurring at approximately 0.6T, where T denotes the orbital period. The two subplots below provide zoomed-in views, showing the positions and velocity directions of particles 1 and 3. For both encounters, the closest distance between particles 1 and 3 is approximately 0.014.

4.3 Outer Solar System

As the final example, we consider the outer Solar System problem, which includes the Sun, Jupiter, Saturn, Uranus, Neptune, and Pluto. The outer Solar System has multiple degrees of freedom. Its clear structure and stable orbits over long timescales make it suitable for evaluating the integrator’s robustness, stability, and computational efficiency. For a detailed description of the masses and orbital parameters of the planets, we refer to Hairer et al. (2006).

We perform simulations over 100 Jupiter orbits using three integrators: SQQ-PN, SQQ-PTN, and SQQ-PTQ. Figure 9 illustrates the evolution of energy error over time. It is worth mentioning that the evolutions of the energy error for SQQ-PTQ and SQQ-PTN are nearly indistinguishable. This is consistent with the phenomenon of the three-body problem mentioned in Section 4.2, because the only difference between SQQ-PTQ and SQQ-PTN lies in the method used to solve the nonlinear equations. Table 5 presents a summary of the CPU time and maximum energy errors for the three integrators. It can be seen that SQQ-PTN has the longest computational time, which is roughly six times that of SQQ-PTQ. This is primarily due to the time transformation, which increases the computational cost of the Jacobian matrix in each iteration. SQQ-PN, which does not use the time transformation, exhibits moderate computational efficiency. In contrast, SQQ-PTQ achieves the highest efficiency by using a quasi-Newton method based on Broyden’s method, avoiding the direct computation of the Jacobian matrix. Even compared to SQQ-PN, SQQ-PTQ requires only 50% of the computational time.

To further evaluate the long-term integration stability of SQQ-PTQ, we performed simulations over 10000 Jupiter orbital periods using both SQQ-PTQ and SQQ-PN. Figure 10 shows that the evolutions of the energy error for both SQQ-PN and SQQ-PTQ oscillate within a bounded range. Table 6 shows the CPU time and maximum energy error of the integration with different integrators. The efficiency of SQQ-PTQ is higher than that of SQQ-PN, even though SQQ-PTQ includes the time transformation. This demonstrates the effectiveness and robustness of combining the time transformation with the quasi-Newton method.

Fig. 6 Evolution of the velocity magnitudes of the three particles over one period in the three-body problem.

Fig. 7 Evolution of the energy error for the integration of the three-body problem by SQQ-PTQ.

Fig. 8 Phase space trajectories of particles in the three-body problem over 500 orbital periods.

Fig. 9 Evolutions of the energy error for the integration of the outer Solar System by SQQ-PN, SQQ-PTN, and SQQ-PTQ, with a total integration time of 100 Jupiter orbital periods.

Fig. 10 Evolutions of the energy error for the integration of the outer Solar System by SQQ-PN and SQQ-PTQ, with a total integration time of 10000 Jupiter orbital periods.

5 Conclusion

This paper introduces an adaptive symplectic integrator SQQ-PTQ, which builds on the fixed-step integrator SQQ. By combining the time transformation with a quasi-Newton method, SQQ-PTQ improves both computational efficiency and applicability.

A key feature of SQQ-PTQ is its use of time transformation to achieve adaptive step sizing. This allows SQQ-PTQ to adjust the step size in gravitational multi-body problems on different timescales. To overcome the challenges of computing Jacobian matrices, SQQ-PTQ employs a quasi-Newton method based on Broyden’s method. This approach significantly accelerates the solution of nonlinear equations, enhancing the overall computational performance.

The effectiveness and robustness of SQQ-PTQ are demonstrated through three numerical experiments. In particular, SQQ-PTQ demonstrates adaptability when handling close-encounter problems. Furthermore, over long-term integrations, SQQ-PTQ maintains energy conservation, further confirming its characteristics as a symplectic algorithm.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Nos. 12311530055 and 12472048), the Shenzhen Science and Technology Program (Grant Nos. ZDSYS20210623091808026), and by High-performance Computing Public Platform (Shenzhen Campus) of SUN YAT-SEN UNIVERSITY.

References

All Tables

All Figures

	Fig. 1 Illustration of the length coordinates for the time intervals [−1, 1] and [ta, tb]. Although the horizontal coordinates of $\hat{\xi }$ and ξ are not equal and their length coordinates λ1 and λ2 are equal, the values of the corresponding interpolation functions are equal (Eq. (12)).
In the text

	Fig. 2 Evolution of the energy error for SQQ-P and SQQ. Here, e = 0.5, and m and n denote the number of interpolation points for the generalized displacement and generalized momentum, respectively.
In the text

	Fig. 3 Evolution of the energy error for SQQ-PTQ for the Kepler problem with high eccentricities.
In the text

	Fig. 4 Phase space trajectories over 500 orbital periods for the Kepler problem with different eccentricities.
In the text

Fig. 5 Trajectories of particles in the three-body problem, with dashed boxes a and b indicating the locations of two close encounters. Box a corresponds to a close encounter occurring at approximately 0.4T, and box b corresponds to the second close encounter occurring at approximately 0.6T, where T denotes the orbital period. The two subplots below provide zoomed-in views, showing the positions and velocity directions of particles 1 and 3. For both encounters, the closest distance between particles 1 and 3 is approximately 0.014.

In the text

	Fig. 6 Evolution of the velocity magnitudes of the three particles over one period in the three-body problem.
In the text

	Fig. 7 Evolution of the energy error for the integration of the three-body problem by SQQ-PTQ.
In the text

	Fig. 8 Phase space trajectories of particles in the three-body problem over 500 orbital periods.
In the text

	Fig. 9 Evolutions of the energy error for the integration of the outer Solar System by SQQ-PN, SQQ-PTN, and SQQ-PTQ, with a total integration time of 100 Jupiter orbital periods.
In the text

	Fig. 10 Evolutions of the energy error for the integration of the outer Solar System by SQQ-PN and SQQ-PTQ, with a total integration time of 10000 Jupiter orbital periods.
In the text