N-body simulation - straight line orbits

Click For Summary
SUMMARY

The discussion focuses on implementing an N-body simulation using Python and Numpy, specifically employing the Euler-Cromer method to simulate the gravitational interactions between the Sun and the first four planets of the solar system. The user encountered issues with straight-line orbits tending towards infinity, attributed to incorrect initial velocities that were too high, resulting in negligible gravitational influence from the Sun. The equations used for acceleration, velocity, and position updates were correctly stated, but the user needed to adjust the initial conditions to achieve realistic orbital behavior.

PREREQUISITES
  • Understanding of N-body simulations
  • Familiarity with Python programming
  • Knowledge of Numpy for numerical computations
  • Basic principles of gravitational physics
NEXT STEPS
  • Review the Euler-Cromer method for numerical integration in physics simulations
  • Learn about gravitational forces and their calculations in N-body systems
  • Investigate the impact of initial conditions on orbital mechanics
  • Explore visualization techniques using Matplotlib for simulation results
USEFUL FOR

Students and developers working on physics simulations, particularly those interested in celestial mechanics and numerical methods for solving differential equations.

XanMan
Messages
14
Reaction score
1

Homework Statement



The problem is your typical N-body simulation, implemented using Python and Numpy. The implementation specifically calls for using the Euler-Cromer method. For this particular case I used the Sun and the first 4 planets of the solar system.

Essentially the problem is I'm getting straight orbits, tending towards infinity. Attached is an image of the output.

Homework Equations


[/B]
The following are the equations used to evaluate the system :
$$\mathbf{a_i} = \sum^N_{j \neq i} \cfrac{Gm_j}{|\mathbf{R}_{ij}|^3}\mathbf{R}_{ij}$$
$$\mathbf{v_i}(t + \Delta t) = \mathbf{v_i}(t) + \mathbf{a_i}\Delta t$$
$$\mathbf{x_i}(t + \Delta t) = \mathbf{x_i}(t) + \mathbf{v_i}\Delta t$$

The Attempt at a Solution


[/B]
Below is my code (commented as best as I can) - I tried various slight alterations, such as calculating the position first, the the velocity, rather than the other way round - but to no avail!

My suspicion is something related to the acceleration - I simply cannot figure out what. I noticed (after printing 10,000 values!) that as the x-coordinate of a body approaches zero, the y-coordinate simply increases more and more, and never goes below it's previous value. In order words I noticed the following issue with regards to the y-positions (which I cannot understand why or if its actually an issue) :
$$y_n < y_{n+1} \ \text{(always)}$$

Appreciate any help - Cheers!

Python:
import numpy as np

from matplotlib import pyplot as plt

# SECTION I - VARS, DATA STRUCTURES

# Defining Constants
AU = 149.6*(10**9)
G = 6.6743*(10**(-11))
deltaT = 24*3600

# Initializing Planetary Data
# planetary_data = [['name', mass, initial position, initial velocity]]
planetary_data = [['Sun',1.989*(10**30),0*AU,0],
                  ['Mercury',0.330*(10**24),0.387*AU,47.36*(10**6)],
                  ['Venus',4.869*(10**24),0.723*AU,35.02*(10**7)],
                  ['Earth',5.974*(10**24),1*AU,29.78*(10**7)],
                  ['Mars',0.647*(10**24),1.524*AU,24.07*(10**7)]]

# Structures to store vel(ocity) vectors and pos(ition) vectors after one full iteration on all bodies
vel_vectors = np.array([[0,0],
                        [0,0],
                        [0,0],
                        [0,0],
                        [0,0]])

pos_vectors = np.array([[0,0],
                        [0,0],
                        [0,0],
                        [0,0],
                        [0,0]])

# Structures to store temp_vel(ocity) vectors and temp_pos(ition) vectors during calculations
temp_vel_vectors = np.array([[0,0],
                             [0,0],
                             [0,0],
                             [0,0],
                             [0,0]])

temp_pos_vectors = np.array([[0,0],
                             [0,0],
                             [0,0],
                             [0,0],
                             [0,0]])
 

Attachments

  • figure_1.png
    figure_1.png
    14.6 KB · Views: 509
Last edited by a moderator:
Physics news on Phys.org
Assuming you use SI base units, your initial velocities are a large fraction of the speed of light (or even above that for Venus). At that speed you won't see a relevant gravitational influence of the Sun.
 
  • Like
Likes   Reactions: XanMan
mfb said:
Assuming you use SI base units, your initial velocities are a large fraction of the speed of light (or even above that for Venus). At that speed you won't see a relevant gravitational influence of the Sun.

Much better now - I feel silly for copying those off my assignment sheet without noticing the mistakes! Cheers!
 

Attachments

  • figure_1.png
    figure_1.png
    27.3 KB · Views: 487
  • Like
Likes   Reactions: mfb

Similar threads

  • · Replies 6 ·
Replies
6
Views
2K
  • · Replies 3 ·
Replies
3
Views
2K
  • · Replies 6 ·
Replies
6
Views
2K
Replies
3
Views
1K
  • · Replies 19 ·
Replies
19
Views
4K
  • · Replies 15 ·
Replies
15
Views
3K
  • · Replies 4 ·
Replies
4
Views
5K
  • · Replies 7 ·
Replies
7
Views
3K
  • · Replies 3 ·
Replies
3
Views
3K
Replies
8
Views
15K