N-body simulation - straight line orbits

Join the discussion
Ask a follow-up here, or get your own question answered by working scientists, mathematicians and engineers — people, not an autocomplete.
Real named experts · corrections over time · the nuance an AI answer skips
2 replies · 3K views
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: 541
Last edited by a moderator:
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: 512
  • Like
Likes   Reactions: mfb