# Understanding the Two-Body Problem with Python Simulation

Written on

## Introduction

In an **inertial frame**, two masses experience a **gravitational attraction** towards each other. The goal of the two-body problem is to analyze how these two objects move in response to their mutual gravitational force. For an overview, refer to the **animation** above, showcasing the simulation developed in this article. To delve deeper into the underlying theory, consult Example 2.2 in the latest edition of *Orbital Mechanics for Engineering Students*.

## Understanding the System

The diagram illustrates the **two-body problem**, featuring two masses, **m1** and **m2**, within an **inertial reference frame**—a frame that is not accelerating. Here, the only force acting on the system is gravity, with no external influences.

Key parameters depicted in the diagram include:

**R1**: the position vector from the origin of the inertial frame to m1, as expressed in Equation 1.**R2**: the position vector from the origin of the inertial frame to m2, defined in Equation 2.**r**: the position vector from m1 to m2, calculated using Equation 3.**G**: the center of mass (COM) of the system. If**m1 > m2**, the COM is closer to m1.**RG**: the position vector from the inertial frame's origin to the COM, expressed in Equation 4.

## Deriving the Equations of Motion

By applying **Newton’s Law of Universal Gravitation** alongside **Newton’s Second Law**, we derive the equations of motion (EOM). These equations, combined with the initial position and velocity vectors, allow us to predict the state of the system over time.

Equation 5 represents the gravitational force, **F**, acting between the two masses.

Equation 6 defines **Newton’s Second Law**, which relates an object's acceleration to the sum of forces acting upon it.

In our two-body system, the only force at play is gravity. Hence, the **force F12** acting from m1 to m2 combines into the two-body EOM represented in Equation 7.

According to Newton's third law, the force from m2 towards m1 is equal and opposite to **F12**. Equation 8 shows the unit vector from object 1 to object 2.

Substituting the unit vector into Equation 7 leads to the final motion equations.

Equation 7 provides a vector representation of acceleration with components in the X, Y, and Z directions. Decomposing Equations 9 and 10 into their components results in Equations 11–16.

## Numerical Integration

To compute the position and velocity of the masses over time, we numerically integrate the equations of motion. The odeint function from the scipy library effectively serves as an integrator. The solver, such as **Runge-Kutta**, requires **first-order ordinary differential equations (ODEs)**. Thus, various **auxiliary variables** are introduced for system modification, as shown in Equations 17–29.

The second-order system of six equations is transformed into **twelve first-order ODEs** using the defined variables, represented in Equations 30–41.

Initial conditions specify the positions and velocities of both masses at the beginning of the simulation. Equations 42–45 outline these conditions as presented in the orbital mechanics textbook utilized in the simulation.

## Python Implementation

Gist 1 showcases the Python code that establishes the **simulation parameters**, including the masses of the bodies and their initial positions and velocities within the **state vector** **y**.

The first-order ODEs of the two-body system are represented as a **Numpy** array in Python, as indicated in Gist 2. Ensure these equations align with Equations 11 through 16.

The propagation of initial conditions over time to solve the two-body problem can be performed with a single line of Python code, as demonstrated below.

# solve two-body problem y = odeint(two_body_eqm, y0, time, args=(G, m1, m2))

The variable y contains the historical state vectors, including the positions and velocities of m1 and m2, from time **t0 = 0** to **t = 480** seconds.

- y[0:3]: X, Y, and Z components of
**R1**at time**t0**. - y[3:6]: Position of m2 in relation to the inertial frame.

Gist 3 provides the Python code for calculating additional position vectors relative to significant reference points, enabling motion visualization from various perspectives, such as from the center of mass or one of the bodies.

## Simulation Results

Figure 2(a) presents the animation of the simulation from two viewpoints.

- The paths of m1 and the center of mass as observed from body 1. Both the COM and body two appear to trace an
**ellipse**. - The orbits of both bodies visible to an observer positioned at the COM, with both orbits being elliptical.

Figure 2(b) illustrates the movement of the two bodies concerning inertial space.

> The two-body system exhibits a periodic spiral motion around the straight-line trajectory of the center of mass through the inertial frame.¹

## Conclusion

Within the inertial frame, the **center of mass follows a straight-line path** with a **constant velocity**, confirming that the **system is not subject to external forces**. This motion continues perpetually.

This article illustrates a technique for deriving and numerically integrating the equations of motion governing two bodies influenced solely by the gravitational force between them. The approach to solving the two-body problem can often be simplified to a single-body scenario. Refer to the **references** below for further analysis.

## References

*Orbital Mechanics for Engineering Students*- Simplifying the Two-Body Problem to a Single Body Problem