- #1
XanMan
- 14
- 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
Last edited by a moderator: