Quantum
Quest

Simulation of the solar system dynamics using Newtonian and Yukawa potentials

In my latest project, I explore the dynamic complexities of the solar system, exploring the gravitational interactions between the Sun and its orbiting planets. By employing both Newtonian gravity and the Yukawa potential as models for these forces, I’ve constructed a simulation that provides a rich understanding of celestial mechanics.

The gravitational interactions within our solar system can be modeled accurately using Newton’s law of universal gravitation, which asserts:

\mathbf{F} = -G \frac{m_1 m_2}{R^3} (\mathbf{r}_2 - \mathbf{r}_1)

Here, G is the gravitational constant, m_1 and m_2 are the masses of two celestial bodies, and \mathbf{r}_1 and \mathbf{r}_2 their respective position vectors. The distance R between these bodies directly influences the magnitude of the gravitational force, which is always directed along the line connecting the two masses.

To introduce modifications to the classical understanding of gravity, I also considered the Yukawa potential, which modifies the gravitational potential by a decaying exponential term:

U(r) = \frac{K_1}{R} \exp\left(-\frac{\left|\mathbf{r}_2 - \mathbf{r}_1\right|}{K_2}\right)

where K_1 and K_2 are constants with K_1 = GmM and K_2 a characteristic length. This results in a force that not only depends on the inverse square of the distance but also includes an exponential factor that diminishes the force over a characteristic range:

\mathbf{F}_{\text{Yukawa}} = \mathbf{F}_{\text{Newton}} \left(1 + \frac{R}{K_2}\right) \exp\left(-\frac{R}{K_2}\right)

The simulation is built on the framework of the n-body problem, where the positions and velocities of celestial bodies are computed using numerical integration. For this purpose, the adaptive Runge-Kutta method (RK45) is employed to handle the differential equations derived from the force equations. Here’s a snippet of the code used to calculate the gravitational force in Python:

def gravitational_force_newton(m1, m2, r_ij):
    r = np.linalg.norm(r_ij)
    return cst.G * m1 * m2 * r_ij / r**3

def gravitational_force_yukawa(m1, m2, r_ij):
    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)

Through these calculations, I assess the system’s kinetic energy, potential energy, and total energy, where:

1. Kinetic Energy (T) is given by:

T = \sum_{i=1}^n \frac{1}{2} m_i v_i^2

2. Potential Energy (U) involving interactions such as:

U = \sum_{i=1}^n \sum_{j=i+1}^n U_{ij}

This detailed approach allows me to explore the stable and chaotic regions of the solar system’s orbital dynamics. The variations in the force laws provide intriguing insights into how minor changes in the force model could affect the stability and trajectories of planetary orbits.

For more insights into this topic, you can find the details here

For access to the complete simulation code, please visit the GitHub repository here.