This simulation models the solar system’s dynamics by calculating the gravitational interactions between the Sun and the planets. The system begins with the Sun positioned at the origin (0, 0, 0) and uses real-world data for each planet’s mass, position, and velocity. The gravitational forces can be modeled using either Newtonian gravity or the Yukawa potential, allowing for the exploration of different force laws. The simulation employs an adaptive Runge-Kutta method (RK45) for numerical integration, which efficiently computes the evolving positions and velocities of the celestial bodies over time. This setup provides a comprehensive and flexible framework for studying the intricate gravitational interactions in our solar system.
The n-body problem involves predicting the individual motions of a group of celestial bodies interacting with each other through forces. Each body experiences a force that depends on the positions, velocities, and possibly other properties of the other bodies. The force acting on the i^{th} body, denoted as \mathbf{F}_i, is the sum of the forces exerted by all other bodies.
Mathematically, the motion of each body can be described by Newton’s second law:
m_i \frac{\mathrm d^2 \mathbf{r}_i}{\mathrm dt^2} = \mathbf{F}_i
where m_i is the mass of the i^{th} body, \mathbf{r}_i is its position vector, and \mathbf{F}_i represents the total force acting on it. The form of \mathbf{F}_i can vary depending on the nature of the interactions between the bodies, such as gravitational forces, electrostatic forces, or other potentials.
For n = 2, there is an exact analytical solution describing the motion of the bodies. However, for n > 2, the problem does not have a general exact solution and can only be solved numerically, due to the complex interactions between the bodies. Numerical methods and simulations are used to approximate the positions and velocities over time.
In classical mechanics, the gravitational force between two point masses m_1 and m_2, separated by a distance R = \left\|\mathbf{r}_2 - \mathbf{r}_1\right\|, is given by Newton’s law of universal gravitation:
F = G \frac{m_1 m_2}{R^2}
To express this force as a vector, pointing from the first body to the second, we use the position vectors \mathbf{r}_1 and \mathbf{r}_2. The vector form of the gravitational force \mathbf{F} exerted on the first body by the second is:
\mathbf{F} = -G \frac{m_1 m_2}{\left\|\mathbf{r}_2 - \mathbf{r}_1\right\|^2} \frac{\mathbf{r}_2 - \mathbf{r}_1}{\left\|\mathbf{r}_2 - \mathbf{r}_1\right\|} = -G \frac{m_1 m_2}{R^3} (\mathbf{r}_2 - \mathbf{r}_1)
Here, the term \mathbf{r}_2 - \mathbf{r}_1 represents the displacement vector from the first body to the second. The factor \frac{1}{\left\|\mathbf{r}_2 - \mathbf{r}_1\right\|^3} ensures the correct magnitude and direction, reducing the force with the square of the distance and maintaining the vector nature of the interaction. This formulation accurately describes the attractive gravitational force, which always acts along the line connecting the two bodies.
Besides the classical Newtonian gravitational force, other central force laws can be applied to the n-body problem. One such example is the Yukawa force, which has a potential of the form:
U(r) = \frac{K_1}{R} \exp\left(-\frac{\left\|\mathbf{r}_2 - \mathbf{r}_1\right\|}{K_2}\right)
where K_1 = GmM is a constant proportional to the product of the gravitational constant G and the masses m and M of the interacting bodies, and K_2 is a constant with dimensions of length. For large values of K_2, the Yukawa potential approximates Newtonian gravity, as K_2 \to \infty leads to the familiar 1/R dependence. The Yukawa potential can give rise to bounded orbits (circular and elliptical) when K_2 is sufficiently large.
The force corresponding to this potential, \mathbf{F}_{\text{Yukawa}}, is given by the negative gradient of the potential:
\mathbf{F}_{\text{Yukawa}} = -\nabla U(r)
Calculating this, we get:
\mathbf{F}_{\text{Yukawa}} = -K_1 \left(\frac{1}{R} + \frac{1}{K_2}\right) \frac{\exp\left(-\frac{R}{K_2}\right)}{R^2} \frac{\mathbf{r}_2 - \mathbf{r}_1}{R}
This force can be expressed in terms of the Newtonian gravitational force, \mathbf{F}_{\text{Newton}}, as:
\mathbf{F}_{\text{Yukawa}} = \mathbf{F}_{\text{Newton}} \left(1 + \frac{R}{K_2}\right) \exp\left(-\frac{R}{K_2}\right) This expression shows how the Yukawa force modifies the Newtonian force by an additional factor dependent on the distance R and the range parameter K_2.
On each body within the n-body system, the total forces acting are the vector sum of all the individual forces exerted by the other bodies. If we denote the force on the i-th body due to the j-th body as \mathbf{F}_{ij}, the total force \mathbf{F}_i on the i-th body is given by:
\mathbf{F}_i = \sum_{\substack{j=1 \\ j \neq i}}^n \mathbf{F}_{ij}
This total force \mathbf{F}_i leads to an acceleration \mathbf{a}_i of the i-th body according to Newton’s second law:
\mathbf{a}_i = \frac{\mathbf{F}_i}{m_i}
where m_i is the mass of the i-th body. The specific form of \mathbf{F}_{ij} depends on the nature of the force between the bodies, such as Newtonian gravity, Yukawa forces, or other potential interactions. The resulting accelerations determine the trajectories of the bodies over time, which can be computed using numerical methods.
The total energy of an n-body system consists of the kinetic energy of all the bodies and the potential energy due to the interactions between them.
T_i = \frac{1}{2} m_i v_i^2
where m_i is the mass and v_i is the velocity magnitude of the i-th body. The total kinetic energy T of the system is the sum of the kinetic energies of all the bodies:
T = \sum_{i=1}^n \frac{1}{2} m_i v_i^2
U_{ij} = -G \frac{m_i m_j}{R_{ij}}
where R_{ij} = ||\mathbf{r}_i - \mathbf{r}_j||. The total potential energy U of the system is the sum of the potential energies of all unique pairs of bodies:
U = \sum_{i=1}^n \sum_{j=i+1}^n U_{ij}
For other forces, like the Yukawa force, the potential energy expression will differ accordingly.
The total energy E of the system is the sum of its kinetic and potential energies:
E = T + U
This total energy is a crucial quantity in the study of dynamical systems, as it is conserved in the absence of external forces, providing insights into the stability and behavior of the system over time.
To numerically integrate the equations of motion in an n-body system, we employ a matrix description of the system’s state. This approach transforms the higher-order differential equations governing the system’s dynamics into a set of first-order differential equations, facilitating their numerical solution.
Given a system described by an nth-order differential equation, it can be recast as a system of n first-order equations. Let \mathbf{X} represent the state vector containing the positions and velocities (or other relevant variables) of all bodies. The derivatives of \mathbf{X} with respect to time are defined by:
\frac{\mathrm d^n \mathbf{X}}{\mathrm dt^n} = \mathbf{Y}_n
where \mathbf{Y}_n is a vector representing the nth derivative of \mathbf{X} with respect to time. By introducing auxiliary variables for each derivative of \mathbf{X}, we can express the system as:
\frac{\mathrm d\mathbf{X}_1}{\mathrm dt} = \mathbf{X}_2, \quad \frac{\mathrm d\mathbf{X}_2}{\mathrm dt} = \mathbf{X}_3, \quad \dots, \quad \frac{\mathrm d\mathbf{X}_{n-1}}{\mathrm dt} = \mathbf{X}_n, \quad \frac{\mathrm d\mathbf{X}_n}{\mathrm dt} = \mathbf{Y}_n
where \mathbf{X}_1 = \mathbf{X} and subsequent \mathbf{X}_i are introduced to reduce the order of differentiation.
In the context of the n-body problem, \mathbf{X} might include the positions and velocities of all bodies, while \mathbf{Y}_n includes the accelerations derived from the forces acting on each body. By setting up the system in this way, we can use standard numerical integration techniques, such as the Runge-Kutta methods or symplectic integrators, to advance the system state in time.
This matrix formulation allows for a unified and systematic treatment of the system’s dynamics, ensuring accurate and efficient computation of the trajectories of the bodies under various force laws.
For a second-order differential equation, such as \mathbf F = m \mathbf a, the system can be reduced to two first-order equations. Using matrix notation, the state vector consists of the position x and velocity v. The system can then be expressed as:
\frac{\mathrm d}{\mathrm dt} \begin{bmatrix} \mathbf x \\ \mathbf v \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} \mathbf x \\ \mathbf v \end{bmatrix} + \begin{bmatrix} 0 \\ \frac{\mathbf F}{m} \end{bmatrix}
Here, the equation \dot{\mathbf x} = \mathbf v is represented by the first row, where the derivative of position with respect to time gives velocity. The second row corresponds to m\dot{\mathbf v} =\mathbf F, representing the acceleration. The matrix on the right-hand side contains the coefficients that relate the state variables \mathbf x and \mathbf v to their derivatives. The force \mathbf F is typically a function of position and velocity, and it is scaled by the inverse of the mass m in the system.
This matrix form is suitable for numerical integration methods, as it clearly delineates the dependencies between position, velocity, and acceleration.
The simulation of the solar system considers the Sun at the position (0, 0, 0) and retrieves the mass, distance, and instantaneous velocity of each planet. The initial data for the celestial bodies are represented in JSON format, as shown below:
{
"sun": {
"label": "Sun",
"mass": 1.989e30,
"scale": 1.5,
"color": [1, 1, 0],
"radius": 6.96342e8,
"position": [0.0, 0.0, 0.0],
"velocity": [0.0, 0.0, 0.0]
},
"mercury": {
"label": "Mercury",
"mass": 3.301e23,
"scale": 2.7,
"color": [0.662, 0.662, 0.662],
"radius": 2.4397e6,
"position": [5.791e10, 0.0, 0.0],
"velocity": [0.0, 4.74e4, 0.0]
},
"venus": {
"label": "Venus",
"mass": 4.867e24,
"scale": 2.7,
"color": [1.0, 0.8, 0.6],
"radius": 6.0518e6,
"position": [1.082e11, 0.0, 0.0],
"velocity": [0.0, 3.502e4, 0.0]
},
"earth": {
"label": "Earth",
"mass": 5.972e24,
"scale": 2.7,
"color": [0, 0, 1],
"radius": 6.371e6,
"position": [1.496e11, 0.0, 0.0],
"velocity": [0.0, 2.978e4, 0.0]
}
}
In this approximation, the planets are initially positioned along the x-axis with their respective velocities perpendicular to the radius vector from the Sun, lying in the same plane. The JSON structure contains relevant physical parameters, such as mass, radius, position, and velocity, which are essential for simulating the dynamics of the solar system under gravitational interactions. The initial setup provides a basis for implementing the numerical integration of the equations of motion to simulate the orbital trajectories of the planets.
The integration routine utilizes an adaptive Runge-Kutta method of 4th/5th order (RK45) to solve the n-body problem. This method dynamically adjusts the step size to maintain accuracy while efficiently computing the trajectories over time. The routine considers whether to use Newtonian or Yukawa gravitational forces, depending on the configuration.
def gravitational_force_newton(m1, m2, r_ij):
"""Compute the gravitational force between two masses using Newton's gravitational law."""
r = np.linalg.norm(r_ij)
return cst.G * m1 * m2 * r_ij / r**3
def gravitational_force_yukawa(m1, m2, r_ij):
"""Compute the gravitational force between two masses using the Yukawa gravitational law."""
r = np.linalg.norm(r_ij)
return cst.G * m1 * m2 * r_ij / r**3 * (1 + r / cfg.yukawa_coeff) * np.exp(-r / cfg.yukawa_coeff)
The gravitational_force function is selected based on the configuration setting cfg.use_yukawa. If cfg.use_yukawa is true, the Yukawa force is used; otherwise, the Newtonian force is applied.
def n_body_problem(t, y, masses):
num_bodies = len(masses)
positions = y[:3 * num_bodies].reshape((num_bodies, 3))
velocities = y[3 * num_bodies:].reshape((num_bodies, 3))
# Initialize accelerations array
accelerations = np.zeros_like(positions)
# Calculate accelerations due to gravity
for i in range(num_bodies):
for j in range(num_bodies):
if i != j:
# Vector from body i to body j
r_ij = positions[j] - positions[i]
# Add contribution of the force to the acceleration
accelerations[i] += gravitational_force(masses[i], masses[j], r_ij) / masses[i]
# Return the derivative of the state vector
derivatives = np.concatenate((velocities.flatten(), accelerations.flatten()))
return derivatives
This function computes the positions and velocities of all bodies and calculates the accelerations due to gravitational forces. The resulting derivatives of the state vector are returned for integration.
def compute(self):
t_eval = np.arange(0, cfg.time_span[1], cst.day_in_seconds)
# Extract masses
masses = np.array([body.mass for body in self._bodies])
# Flatten initial conditions
initial_conditions = np.concatenate([np.concatenate([body.position for body in self._bodies]
+ [body.velocity for body in self._bodies])])
num_bodies = len(masses)
# max_step is set to ensure convergence of RK45
solution = solve_ivp(n_body_problem, cfg.time_span, initial_conditions, args=(masses,),
method='RK45', t_eval=t_eval, max_step=10 * cst.day_in_seconds)
self._positions = solution.y[:num_bodies * 3].reshape((num_bodies, 3, -1))
self._velocities = solution.y[num_bodies * 3:].reshape((num_bodies, 3, -1))
The compute method integrates the system over the specified time span, cfg.time_span. The solve_ivp function from SciPy is used with the RK45 method, providing accurate solutions for the positions and velocities of the bodies. The time points for evaluation (t_eval) are defined, and the initial conditions are prepared from the mass and state vectors of the bodies. The solution is then reshaped and stored for further analysis and visualization.
The use of an adaptive Runge-Kutta method allows for efficient and accurate simulations of complex gravitational interactions in the n-body system, taking into account both traditional and modified gravitational forces.
To enhance the realism of the simulation, I used data from NASA’s JPL HORIZONS system. This system provides customizable production of accurate ephemerides for observers, mission-planners, researchers, and the public by numerically characterizing the location, motion, and observability of solar system objects over time. The data was retrieved using a script to ensure precise and up-to-date information on the planets’ positions and velocities.
Additionally, the simulation includes an option to visualize the orbits of the planets, along with their projections on a 2D plane. This feature allows for a clear and detailed representation of the solar system’s dynamics.
The final results, including the simulation code and data, are available for download on my GitHub repository here.