Comparative Study of Runge-Kutta (RK4) vs. Adams-Bashforth-Moulton for Solving the Three-Body Problem

Introduction Study Background

The three-body problem describes the motion of three-point masses, and according to Newton's law of universal gravitation, these objects interact through gravitational forces. Many scientists have been interested in this problem, such as Newton, who laid the foundation, and Leonhard Euler and Joseph-Louis Lagrange, who discovered special equilibrium solutions, and it has also attracted the attention of Henri Poincaré [1] . In his work, for which he won the King Oscar II prize, Poincaré demonstrated that the three-body problem lacks a general analytical solution, and proved that it exhibits sensitive dependence on initial conditions (SDIC), as this concept has become one of the main pillars of chaos theory [2] . This discovery has contributed to making numerical methods the primary tool for understanding the dynamics of astronomical systems. The greatest difficulty in solving the differential equations governing the three-body problem is their non-integrability, which means the lack of sufficient first integrals to reduce the degrees of freedom precludes simple and periodic solutions in the general case, and thus approximating the solutions with numerical precision, which is a process highly sensitive to the integration algorithm used, remains the only feasible option [1] .

Using Numerical Methods in Solving the Three-Body Problem

The use of numerical integration has become essential for studying the temporal evolution of astronomical systems, as closed-form analytical solutions do not exist except in highly specialized cases such as equilibrium configurations (Lagrange points) or particular mass ratios, and short-term approximate solutions cannot generally be relied upon for applications such as computing spacecraft trajectories, assessing the stability of planetary systems, or modeling stellar dynamics in open clusters [3] . Press et al. (2007) [3] indicate that employing an algorithm with insufficient accuracy or inadequate stability may lead to catastrophically erroneous physical predictions, such as spurious orbital ejection of a planet or non-physical collision events, due to the accumulation of truncation errors over extended integration intervals. In Hamiltonian systems, this issue is significantly exacerbated: while energy conservation is theoretically preserved by the exact flow, in the absence of symplectic integrators, one must resort to comparing numerical integrators that exhibit bounded energy drift or superior long-term stability properties.

Methods of Solving Numerical Differential Equations (ODEs)

Numerical methods for solving ordinary differential equations (ODEs) are generally divided into two main categories: Runge-Kutta methods and Adams-Bashforth-Moulton methods.

1.2.1

RK Methods Named after the German mathematicians Carl Runge and Wilhelm Kutta, these methods are the most widely used class of one-step integrators. The RK4 method is characterized by the fact that it does not require solution values from previous steps, as it relies directly on the current state of the system to compute the next step, which makes it flexible and easy to apply, particularly when adjusting the step size. The method evaluates the derivative function at each time step, which increases the computational cost, especially in complex cases [4] .

1.2.2

Adams-Bassforth-Moulton Methods (ABM Methods) These methods were developed by John Couch Adams and Francis Bashforth in the 19th century, and were later extended by Forest Ray Moulton. ABM belongs to the category of multistep methods. These methods rely on several previous steps in order to calculate the solution at the current step [5] . The method proposed in this research consists of a predictor step and an implicit corrector step, as this method is characterized by high efficiency. It requires only two function evaluations per step after the initialization phase, instead of four as in the case of RK4, which makes it computationally faster and reduces the computational cost in many cases. However, the main disadvantage of this method is that it requires historical solution data and necessitates a separate startup algorithm to initialize the integration [4] .

Research Objectives and Research Questions

The research aims to conduct an in-depth applied comparison between the RK4 and ABM methods in the context of solving the three-body problem, and to clarify the behavior of each of the two methods under chaotic conditions and close approaches. The research attempts to answer the following questions: 1 How does the method structure (one-step vs. multistep) affect energy drift in the short term (10 integration periods) compared to the long term (1000 integration periods)? 2 What is the relationship between accuracy and computational cost?

3 How does each method behave when confronted with singularities in the Pythagorean problem, where bodies approach close encounter regimes with r ij → 0?

Research Structure

This paper is divided into six sections; in addition to this first section, which provides a general introduction to the research and explains its objectives, numerical tools, and methodology, the second section is devoted to the theoretical framework, where the equations of motion are derived and conservation laws are explained. The third section explains the mathematical details of the RK4 and ABM methods, including the derivation of their coefficients. The fourth section describes the experimental methodology and test cases. The fifth section presents the results and error analysis.

The sixth section discusses the results in the context of chaos and numerical stability.

Look at the issue of the three Assams 2.1 Equations of motion in Cartesian coordinates

Assume a system of three-point masses moving in threedimensional space. Let m 1 , m 2 , m 3 the masses, and let r i be the position vector of each body within an inertial reference frame. According to Newton's law of universal gravitation, the force exerted on body i by body j s given by the relation [1] :

where G is the gravitational constant. Applying Newton's second law of motion (F = ma), the equation of motion for each body can be written as follows [1] :

where r̈i represents the acceleration vector.

When the distance between two bodies approaches zero (∥ r j -r i ∥→ 0), a singularity occurs, which corresponds to an infinite force. This system of second-order differential equations must be converted into a system of first-order equations in order to obtain a numerical solution. The velocity vector is defined by the relation [5] :

This yields a system of first-order ordinary differential equations, in which each body is described by three position variables and three velocity components in three-dimensional space. The complete state of the system can be expressed in the state vector y [5] :

This is written in the standard form dy dt = f(t, y), which is required for the application of the Runge-Kutta and Adams-Bashforth-Moulton methods.

Conserved Quantities and Accuracy Standards

Conserved physical quantities serve as a primary metric for assessing the quality of a numerical integrator in an isolated system, since a high-quality integrator exhibits minimal drift in these quantities over time [3] . Total energy is defined as the sum of the kinetic energy and gravitational potential energy of the bodies, and is expressed mathematically as follows [4] :

his quantity is conserved in the exact solution, with dE dt = 0 theoretically. However, in numerical integration, a relative energy error emerges, which serves as a metric for integrator quality. This error is monitored at each integration step; if it grows linearly with the number of steps, the integrator is unsuitable for long-term orbital simulations ΔE/E 0 [4] . Total angular momentum represents the conserved rotational quantity of the system and is expressed by the following relation [3] :

This is also a conserved quantity; deviations in its magnitude ∥ L ∥ indicate structural deficiencies or non-symplectic behavior in the integrator.

Chaotic Nature and Susceptibility to Errors

The three-body problem is characterized by the presence of positive Lyapunov exponents, so any error, no matter how small, whether it is a measurement error or an approximation error in one of the six state variables describing the system, will increase exponentially with time, as in the following relation [1] :

where λ is the Lyapunov exponent. Equation (7) shows that after a certain period of time, the numerical solutions will begin to diverge significantly from each other, whether two separate integrators are used or the same integrator is employed but with different step sizes [1] . As a result of equation ( 7), studies of chaotic systems focus on the physical and statistical properties of the system under investigation, such as energy conservation. The ABM method usually exhibits better stability in long-term energy conservation compared to the RK4 method [5] .

3. Mathematical Derivation of Numerical Methods

Classic Fourth-Order Runge-Kutta Method (RK4)

The RK4 method is a standard method in the numerical solution of ordinary differential equations. The basic idea of this method is to approximate the derivative over the interval [t n , t n+1 ] y a weighted mean of four slope evaluations at different points within that interval [4] . The general form of a first-order ordinary differential equation is first considered:

The subsequent value y n+1 is calculated by the relationship:

The stage values k i are calculated as follows:  k 4 = f(t n + h, y n + hk 3 ): Slope at the end of the interval, using the estimate from k 3 .

These expressions were derived by matching the Taylor series expansion of the exact solution y(t + h) up to the h 4 term. The local truncation error is of order O(h 5 ), while the global error is of order O(h 4 ). The absolute stability region is one of the most important features of the RK4 method, as it covers a portion of the negative real axis in the complex plane, which makes it very suitable for non-stiff problems [4] . On the other hand, the energy error accumulates linearly over time in periodic orbits, making it unsuitable for long-term integrations [3] .

Fourth-Order Adams-Bashforth-Moulton (ABM) Method

These methods are based on previous values and results, i.e., the stored function values at previous grid points, as they are derived from the integration of the following differential equation:

Within the integration, f(t, y(t)) is replaced by an interpolating polynomial that passes through the previous points P(t)(t n , f n ), (t n-1 , f n-1 ), … A third-order polynomial P 3 (t) that passes through the points f n , f n-1 , f n-2 , f n-3 to form the explicit Adams-Bashforth predictor is used, yielding the following relationship [5] :

Equation ( 11 into the function to compute an initial estimate f n+1 (p) , and then applying the following corrector formula [5] :

The local truncation error constants of the fourth-order Adams-Bashforth method (251/720) and the implicit Adams-Moulton method (19/720) are derived from the truncation error analysis of the interpolating polynomial, which are standard values documented in the classical literature of numerical analysis [5] . When comparing the explicit formula in equation ( 11) with the implicit formula in equation ( 12), a difference in coefficients is observed, because the implicit formula possesses a smaller error constant than the explicit formula, making it more accurate and stable [5] . The complete ABM algorithm [4] proceeds as follows:

1 Retrieval of pre-stored values: Only one new function evaluation is required per step (after initialization). Compared to RK4, which requires four function evaluations, ABM is theoretically significantly faster. However, it is important to consider the computational overhead of storing previous step values and the complexity of stepsize adaptation in the ABM method [4] .

Startup Procedures

The ABM method requires the solution values at four previous grid points to initiate the multistep integration, but at the initial time 𝑡 0 ,only the initial conditions are available. In this study, for 𝑛 = 0,1,2,3) at 𝑡 = 0 with initial value 𝑦 0 , the RK4 method was employed to generate the first three solution values 𝑦 1 , 𝑦 2 , 𝑦 3 )with high accuracy, using the same step size ℎ to ensure consistency of the local truncation errors. This initialization procedure incurs a small computational overhead at the beginning of the simulation (a negligible cost for longterm integrations), but is necessary to leverage the computational efficiency of the ABM method in subsequent steps [3] .

Methodology

he experiments were designed to ensure a comprehensive assessment by covering two contrasting dynamical regimes: the first being a stable periodic configuration (for quantifying energy drift), and the second being a chaotic configuration (for assessing numerical stability and stiffness handling). This design establishes a practical framework for testing and evaluating the numerical methods described above.

Selection of Test Cases

The Figure-8 periodic orbit and Burrau's problem were selected as benchmark problems commonly used to compare numerical integrators in the scientific literature.

The Figure-8 Periodic Orbit

This periodic solution was discovered by Chenciner and Montgomery (2000) [6] . This configuration assumes three point masses with equal mass m 1 = m 2 = m 3 = 1 moving along the same closed trajectory that forms a figure-eight shape, as shown in Figure 1. This orbit is considered linearly stable, as it exhibits robustness to small perturbations and does not immediately destabilize under minor disturbances. Initial conditions (in dimensionless units) were defined as follows [6] :

Since the orbit is strictly periodic with period T ≈ 6.3259, any numerical error will manifest clearly as a deviation from the closed trajectory. This benchmark is well-suited for assessing long-term energy conservation properties [3] .

Burrau's Problem

This problem was first investigated by Burrau in 1913, involving three-point masses with mass ratios (m 1 = 3, m 2 = 4, m 3 = 5), initially positioned at the vertices of a rightangled triangle with sides 3, 4, and 5, and with zero initial velocities (v = 0) [5] , as shown in Figure 2. The initial conditions are defined as follows [5]: Mass m 1 = 3 at (1 , 3), Mass m 2 = 4 at (-2, -1), Mass m 3 = 5 at (1, -1), with all initial velocity components equal to zero. This system is characterized by its rapid chaotic evolution, in which the bodies undergo close encounters to form a tight binary pair, while the third body is ejected via gravitational scattering. This scenario involves extremely close approaches where gravitational forces become singular, representing a stringent test of numerical stability and stiffness handling in integration algorithms (Close-Encounter Robustness Test) [5] .

Numerical Implementation and Experimental Setup

The simulations were implemented using the MATLAB programming environment (R2025b), with standardized test conditions as follows: Total Integration Time: yielding T_total ≈ 632.59 (dimensionless time units).  Pythagorean problem: Integrated from t = 0 to t = 20, which is sufficient for binary formation and ejection to occur.

Step Sizes (h): Several values were tested: h = 0.001,0.005,0.02,0.08, in order to observe the behavior of the global error as the step size varies. Units: Canonical (dimensionless) units with (G = 1) were employed to avoid numerical errors associated with extremely large or small physical quantities.

Evaluation Criteria

He following quantitative metrics were adopted to assess numerical performance [3] , [5] :

Relative energy error δE(t) =|

, where E 0 denotes the total energy at t = 0. he quantity log 10 (δE) is plotted against time to monitor the rate of error growth.

Computational efficiency: Measured via CPU time using MATLAB's tic and toc functions, together with the total number of function evaluations (i.e., acceleration vector computations f(t, y)).

Results and analysis

Numerical simulations were performed using both the fourthorder Runge-Kutta method (RK4) and the Adams-Bashforth-Moulton method (ABM) for the two test cases. The results are summarized in Tables 1 and2, while Figures 3 to 6 illustrate the geometric behavior of the trajectories and the temporal evolution of the energy error. The trajectories in the standard configuration with step size h = 0.01 are illustrated as follows:

Stable Orbit Energy Conservation Analysis (Figure-8)

The Figure-8 orbit is a periodic and linearly stable configuration, making it an ideal benchmark for quantifying energy drift over medium and long integration intervals. The results in Table 1 show a clear difference in the performance of the two methods. The higher accuracy of RK4 is evident when using a step size h = 0.01 and an integration interval t = 20, with RK4 recording a maximum energy error of 8.33 × 10⁻⁹, while ABM recorded an error of 8.97 × 10⁻⁷.

That is, RK4 was nearly two orders of magnitude more accurate in conserving the total energy of the system. Temporal error evolution (Figure 4): The logarithmic plot of the energy error shown in Figure 4 indicates that both methods are non-symplectic, resulting in a gradual drift of error over time. The RK4 curve (solid blue) remains below the ABM curve (dashed red) by at least one order of magnitude throughout the simulation period (t = 200), confirming its superior numerical stability in periodic systems.

The literature suggests that multistep methods may achieve slower error growth over very long integration intervals (T ≫ 10³) , but the current result (T ≤ 200) show that ABM exhibits a larger local truncation error constant during the initialization phase, which increases the initial magnitude of energy fluctuations. RK4 relies on function evaluations at the current step only, which reduces the accumulation of historical truncation errors and achieves higher instantaneous accuracy within this time frame.

Computational efficiency and the effect of step size

Efficiency was measured via CPU time recorded in Tables 1 and2, where ABM was 40% to 65% faster compared to RK4 across all step sizes and integration intervals. For example, in the case of the Figure-8 orbit at h = 0.01 and t = 200, ABM required 0.1217s versus 0.2150s for RK4. This behavior corresponds to theoretical computational complexity, as RK4 requires four function evaluations of f(t, y) per time step, while ABM (after the startup phase) requires only two evaluations (one for prediction and one for correction). The reduction in function evaluations clearly outweighs the computational overhead resulting from storing and managing previous step values. The two tables show that reducing h lowers the energy error of both methods at a rate consistent with fourth-order convergence, confirming the correctness of the numerical implementation and the convergence of both algorithms (O(h⁴)) Any decrease in h increases the computation time approximately linearly, highlighting the classic trade-off between accuracy and computational cost.

Dynamic Response in the Chaotic Order (The Pythagorean Question)

The dynamics in the Pythagorean problem are fundamentally different due to the presence of close encounters and abrupt variations in gravitational acceleration, which serves as a stringent test of algorithmic robustness. Table 2 shows that the energy errors increased by several orders of magnitude (from 10¹ to ³ 10 ) compared to the stable configuration. At h = 0.01, the RK4 error was approximately 42.6, while the ABM error increased to 1.29 × ³ 10 . This substantial difference reflects the high sensitivity of chaotic systems to initial conditions and the nature of the numerical integration. Figures 3 and4 show that RK4 handles the close-encounter phase with more controlled trajectory divergence, while ABM exhibits faster and less uniform error growth, suggesting a more rapid loss of numerical stability during abrupt variations in the derivative function. The logarithmic plot (Fig. 6) shows that the energy error does not grow indefinitely, but rather saturates at characteristic levels (~⁰ 10 for RK4 and ~³ 10 for ABM). This behavior is expected in Lyapunov-based chaotic systems, where nearby trajectories diverge exponentially, making precise energy conservation impossible when using fixed-step, nonsymplectic integrators. RK4 exhibits greater robustness in handling sharp variations in the derivative function during close encounters, because it relies solely on the current state without accumulating truncation errors from previous steps, which provides more precise control over error growth near singularities.

Reconciling Numerical Results with the Theoretical Framework

The experimental results confirm the validity of the theoretical hypotheses presented in Sections 2 and 3, while clarifying some nuances. ABM excels computationally by reducing function evaluations, while RK4 excels in accuracy and stability, especially in extreme conditions. This is consistent with the analysis of Press et al. (2007) [3] regarding the computational cost of explicit versus multistep algorithms. The apparent discrepancy between the long-term theoretical superiority of ABM and its weaker performance in these simulations is due to the fact that the local truncation error constant of the Adams-Bashforth predictor (251/720) is greater than that of the Adams-Moulton corrector (19/720), raising the initial error amplitude. In intermediate-time simulations (t ≤ 200) of chaotic systems, this disadvantage overshadows the potential long-term drift advantage. The local truncation error constant of the corrected Adams-Moulton method (251/720 ≈ 0.026), although smaller than the Adams-Bashforth coefficient (251/720 ≈ 0.348), remains greater than the effective error constant of RK4 in non-stiff systems, which explains the slower accumulation of error in the single-step method under the same step size h [4] . RK4 is recommended when dealing with chaotic systems, close encounters, or when the accuracy of conserving physical invariants is a primary priority, even at the expense of computational time. ABM is recommended for very long-term simulations of smooth systems, where its computational efficiency outweighs the allowable energy error tolerance, provided that the step size is carefully controlled.

Discussion 6.1 Trade-offs

This study illustrates the classic trade-off in numerical analysis. The RK4 method represents a general-purpose choice, owing to its simplicity, ease of implementation, minimal memory requirements for storing solution history, and favorable stability region for non-stiff problems. ABM, by contrast, is a specialized choice for long-term simulations of smooth systems. The primary trade-off lies between implementation complexity and computational efficiency: while RK4 can be implemented in a few lines of code, ABM requires careful management of time arrays and historical solution values. According to Press et al. (2007) [3] , investing in ABM implementation yields significant benefits in astrophysical simulations that span extended integration intervals (e.g., modeling solar system formation), where saving 50% of computational time becomes crucial, particularly when balanced against short-term energy drift behavior.

The Impact of Chaos on the Reliability of Results

Chaotic systems necessitate a redefinition of the concept of numerical precision, as RK4 and ABM will produce substantially divergent trajectories after a relatively short integration time, due to positive Lyapunov exponents. Consequently, relying on exact geometric position agreement after long integration periods as the sole measure of validity is insufficient.

The conservation of the system's structural invariants serves as the most reliable indicator of integrator quality, such as the stability of total energy and the magnitude of angular momentum. If an algorithm maintains these conserved quantities within narrow bounds, it exhibits shadowing of a physically plausible solution to the problem, even if the apparent trajectory diverges [5] . From this theoretical perspective, ABM is expected to provide more stable verylong-term energy behavior (T ≫ 10 3 ) compared to the cumulative linear drift characteristic of RK4, although the numerical results in this study (for t ≤ 200) demonstrated the superiority of RK4 in instantaneous accuracy, reflecting the importance of selecting an appropriate evaluation metric aligned with the target simulation interval.

Study Limitations

This study identifies a set of constraints that delineate the scope of applicability of its findings, the most significant of which are: Fixed Step Size: A fixed step size was adopted in all simulations. Advanced applications may employ adaptive step-size control techniques, which could reduce the performance gap between the two methods, particularly during close-encounter phases in chaotic systems. Comparison Limited to Fourth Order: Only fourth-order methods were compared. Higher-order methods (e.g., RK8 or Dormand-Prince 5(4)) or symplectic integrators specifically designed for Hamiltonian systems may offer superior performance for the three-body problem in terms of long-term preservation of phase-space structure. Both the RK4 and ABM methods employed in this study are classified as non-symplectic integrators, which accounts for the gradual energy drift observed in long-term simulations, particularly in conservative Hamiltonian systems. In recent numerical mechanics literature, symplectic integrators or regularized methods are preferred for the long-term integration of celestial mechanical systems, as they preserve the geometric structure of phase space and ensure enhanced stability of motion invariants over extended integration intervals [5] . The incorporation of adaptive step-size control mechanisms can improve the efficiency of both methods in handling close encounters without compromising numerical stability. The integration of these advanced methodologies (high-order symplectic methods with adaptive control) represents a promising direction for future research to enhance the accuracy and efficiency of simulations of chaotic multi-body systems.

Conclusion

Theoretical analyses and numerical experiments (based on the Figure-8 and Pythagorean test cases) have demonstrated that the fourth-order Adams-Bashforth-Moulton (ABM) method significantly outperforms the Runge-Kutta method (RK4) in terms of conserving the total energy of the system and computational efficiency in long-term simulations. This superiority is attributed to the efficiency of the corrector step and the reduction in the number of function evaluations required per integration step. The RK4 method remains a powerful and reliable choice for applications involving rapid variations in velocity or close encounters, owing to its robustness as a single-step method and its ease of implementation. Based on these findings, the study recommends the use of the ABM method in astronomical simulations that require high accuracy in preserving conserved quantities (energy and angular momentum) over extended integration intervals (such as assessing the stability of planetary systems), while recommending the use of RK4 in the study of short-term dynamics, during preliminary investigations or algorithm prototyping, or in cases involving chaotic dynamics and close encounters between bodies.

References

1. Liao S, Li X. On the periodic solutions of the three-body problem. Natl Sci Rev. 2019;6(6):1070-1071. DOI:10.1093/nsr/nwz102.
2. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes: the art of scientific computing. 3rd ed. Cambridge: Cambridge University Press; 2007.
3. Butcher JC. Numerical methods for ordinary differential equations in the 20th century. J Comput Appl Math. 2000;125(1-2):1-29. DOI:10.1016/S03770427(00)00455-6.
4. Hairer E, Nørsett SP, Wanner G. Solving ordinary differential equations I. 2nd ed. Berlin, Heidelberg: Springer; 1993. DOI:10.1007/978-3-540-78862-1.
5. Chenciner A, Montgomery R. A remarkable periodic solution of the three-body problem in the case of equal masses. Ann Math. 2000;152(3):881-901. DOI:10.2307/2661357.