Three body movement simulator
Simulates the chaotic gravitational motion of three bodies by numerically integrating Newtonian equations using RK4 for accurate trajectory prediction.
Like a lobster shell, security has layers — review code before you run it.
License
SKILL.md
Skill: 3-Body Movement Simulator
Name: 3-body-simulator
Description: A comprehensive guide to architecting and implementing a numerical simulator for the three-body problem, covering the physics, mathematical formulation, and computational strategies required to model chaotic gravitational systems.
Keywords: ["three-body problem", "physics simulation", "numerical integration", "chaos theory", "python", "astrophysics", "computational physics", "newtonian gravity"]
3-Body Movement Simulator
Objective To guide users through the process of building a functional and accurate simulator for the three-body problem, transforming theoretical physics into a visual, interactive computational model.
Core Concept: The Dance of Chaos
The three-body problem is a classic problem in physics and celestial mechanics that seeks to predict the motion of three celestial bodies interacting through their mutual gravitational attraction. Unlike the two-body problem, which has a clean, analytical solution (elliptical orbits), the three-body problem is famously chaotic. This means that even minuscule differences in the initial positions and velocities of the bodies can lead to vastly different outcomes, making long-term prediction impossible.
- The Goal: To create a computer program that calculates the trajectories of three bodies over time by numerically solving the equations of motion.
- Key Challenge: The system is highly sensitive to initial conditions. Our simulator must use precise mathematical methods to ensure the results are physically plausible and stable over the simulation's duration.
The Physics: Newton's Law of Universal Gravitation
The foundation of our simulator is Newton's law of gravitation. The force exerted on body i by body j is given by:
Fᵢⱼ = (G * mᵢ * mⱼ / |rᵢⱼ|³) * rᵢⱼ
Where:
- Fᵢⱼ is the force vector on body
ifrom bodyj. - G is the gravitational constant.
- mᵢ and mⱼ are the masses of the bodies.
- rᵢⱼ is the vector pointing from body
ito bodyj(i.e., rⱼ - rᵢ). - |rᵢⱼ| is the distance between the two bodies.
To find the total force on any single body, we must sum the gravitational forces exerted on it by the other two bodies.
Computational Strategy: Numerical Integration
Since we cannot solve the equations of motion analytically, we must use numerical integration. This involves breaking time into tiny steps (dt) and calculating the new position and velocity of each body at each step.
- The Naive Approach (Euler Method): This is the simplest method but is highly inaccurate and unstable for orbital mechanics. It tends to add energy to the system, causing orbits to spiral outwards. We will avoid this method.
- The Recommended Approach (Runge-Kutta 4th Order - RK4): RK4 is a robust and widely used method that offers a good balance between accuracy and computational cost. It calculates the slope of the solution at multiple points within a time step to produce a much more accurate result than the Euler method.
Step-by-Step Implementation Guide
Step 1: Define the System State
We need a way to represent the state of our entire system at any given moment. The state consists of the position and velocity of all three bodies. We can represent this as a single state vector S.
S = [x₁, y₁, z₁, vx₁, vy₁, vz₁, x₂, y₂, z₂, vx₂, vy₂, vz₂, x₃, y₃, z₃, vx₃, vy₃, vz₃]
This vector contains 18 elements: 3 position coordinates and 3 velocity components for each of the 3 bodies.
Step 2: Implement the Derivative Function
This is the core physics engine. This function takes the current state vector S and time t as input and returns the derivative of the state vector, dS/dt. The derivative of position is velocity, and the derivative of velocity is acceleration.
import numpy as np
G = 6.67430e-11 # Gravitational constant
def calculate_derivatives(t, state, masses):
"""
Calculates the time derivatives of the state vector.
Args:
t (float): Current time (not used directly as gravity is time-independent)
state (np.array): The current state vector [x1, y1, z1, vx1, ...]
masses (list): List of masses [m1, m2, m3]
Returns:
np.array: The derivatives of the state vector [vx1, vy1, vz1, ax1, ...]
"""
# Unpack the state vector for clarity
r1 = state[0:3]
v1 = state[3:6]
r2 = state[6:9]
v2 = state[9:12]
r3 = state[12:15]
v3 = state[15:18]
m1, m2, m3 = masses
# Calculate the force on body 1 from bodies 2 and 3
r12 = r2 - r1
r13 = r3 - r1
F1 = (G * m1 * m2 / np.linalg.norm(r12)**3) * r12 + \
(G * m1 * m3 / np.linalg.norm(r13)**3) * r13
# Calculate accelerations (F = ma -> a = F/m)
a1 = F1 / m1
# Repeat for body 2 (forces from 1 and 3)
r21 = r1 - r2
r23 = r3 - r2
F2 = (G * m2 * m1 / np.linalg.norm(r21)**3) * r21 + \
(G * m2 * m3 / np.linalg.norm(r23)**3) * r23
a2 = F2 / m2
# Repeat for body 3 (forces from 1 and 2)
r31 = r1 - r3
r32 = r2 - r3
F3 = (G * m3 * m1 / np.linalg.norm(r31)**3) * r31 + \
(G * m3 * m2 / np.linalg.norm(r32)**3) * r32
a3 = F3 / m3
# Construct and return the derivative vector
# d(position)/dt = velocity
# d(velocity)/dt = acceleration
dstate_dt = np.array([*v1, *a1, *v2, *a2, *v3, *a3])
return dstate_dt
Step 3: Implement the RK4 Integrator
This function will use the derivative function to advance the state of the system by one time step dt.
def rk4_step(t, state, dt, masses):
"""
Performs a single Runge-Kutta 4th order integration step.
"""
k1 = calculate_derivatives(t, state, masses)
k2 = calculate_derivatives(t + dt/2, state + dt/2 * k1, masses)
k3 = calculate_derivatives(t + dt/2, state + dt/2 * k2, masses)
k4 = calculate_derivatives(t + dt, state + dt * k3, masses)
new_state = state + (dt / 6) * (k1 + 2*k2 + 2*k3 + k4)
return new_state
Step 4: The Main Simulation Loop Now we bring it all together. We set the initial conditions, define the time step, and run the loop.
def run_simulation():
# 1. Set Initial Conditions
# Example: A figure-8 orbit (a known stable solution)
# These are dimensionless units for simplicity
masses = [1.0, 1.0, 1.0]
initial_state = np.array([
0.97000436, -0.24308753, 0, 0.466203685, 0.43236573, 0, # Body 1
-0.97000436, 0.24308753, 0, 0.466203685, 0.43236573, 0, # Body 2
0.0, 0.0, 0, -0.93240737, -0.86473146, 0 # Body 3
])
# 2. Simulation Parameters
t = 0.0
dt = 0.001 # Time step
total_time = 10.0
# 3. Run the Loop
current_state = initial_state
trajectory = [current_state.copy()] # Store states for visualization
while t < total_time:
current_state = rk4_step(t, current_state, dt, masses)
trajectory.append(current_state.copy())
t += dt
return np.array(trajectory)
# To visualize, you would extract the positions from the trajectory array
# and plot them using a library like Matplotlib.
Best Practices and Considerations
- Choose an Appropriate Time Step (
dt): A smallerdtincreases accuracy but makes the simulation slower. A largerdtis faster but can lead to instability and energy drift. You must find a balance. - Monitor Energy Conservation: In a closed gravitational system, the total energy (kinetic + potential) should remain nearly constant. Calculating the total energy at each step is a great way to check the accuracy and stability of your integrator. If the energy is steadily increasing or decreasing, your
dtis likely too large. - Start with Known Solutions: Test your simulator with known stable or periodic solutions, like the figure-8 orbit, to verify that your implementation is correct before experimenting with random initial conditions.
- Visualization: The results are just arrays of numbers. Use a plotting library like Matplotlib (for Python) to create 2D or 3D plots of the trajectories to see the beautiful and chaotic dance of the three bodies.
Files
1 totalComments
Loading comments…
