Simulating Earth's orbit using symplectic Euler in Python

Question

I'm trying to simulate the motion of earth around the sun. (This is the task I am trying to do) This is what I've come up with so far

import numpy as np
import matplotlib.pyplot as plt
#Set parameters:
N = 365                     # Earth days in a year
dt = 1.00/N                 # Time Step: Fractions of a year - 1 Earth day (i.e. 1/365)
mu = 4 * np.pi**2           # mu=4pi^2 is the Gravitational Parameter: mu = GM where G=6.67e-11 is the Universal Gravitational Constant and M is the mass of the body
rEar = 1

#Create an array, for all variables, of size N with all entries equal to zero:
xEar = np.zeros((N,))
yEar = np.zeros((N,))
vxEar = np.zeros((N,))
vyEar = np.zeros((N,))

# Initial Conditions:
xEar[0] = rEar                   # (x0 = r, y0 = 0) in AU
vyEar[0] = np.sqrt(mu/rEar)      # (vx0 = 0, v_y0 = sqrt(mu/r)) AU/yr

#Implement Verlet Algorithm:
for k in range(0,N-1):
    vxEar[k+1] = vxEar[k] - (mu*xEar[k]) / (rEar**3)*dt
    xEar [k+1] = xEar[k] + vxEar[k+1]*dt
    vyEar[k+1] = vyEar[k] - (mu*yEar[k]) / (rEar**3)*dt
    yEar [k+1] = yEar[k] + vyEar[k+1]*dt

#Plot:
plt.plot(xEar, yEar, 'go')
plt.title ('Circular Orbit of Earth')
plt.xlabel ('x')
plt.ylabel ('y')
plt.axis('equal')
plt.show()

but I don't think my implementation of verlet is quite right? With this code I think the orbit will always be circular. Any tips to improve this would be much appreciated

The implementation is correct. You can spare some memory by using a vector of 2 elements for the velocities (you need only current velocity and next velocity), in fact, you don't need them for the plot. Regarding the circular orbits, this depends on your initial speed, for instance ,if you put a 2 before the expression of VyEar[0] you get an eliptic orbit. — Fabrizio, Mar 26 '20 at 14:11

Lutz Lehmann · Answer 1 · 2020-11-19T13:45:53.243

1

Your text is slightly wrong. Symplectic Euler is an order 1 method, while Stormer-Verlet is order 2. One could bridge the gap between the methods by implementing the leap-frogging Verlet scheme where the velocities are taken at the midpoints of the time intervals. In practice that means that the initial velocity has to be shifted to the one at time t0-dt/2, everything else remains the same.

# Initial Conditions:
xEar[0] = rEar                   # (x0 = r, y0 = 0) in AU
vyEar[0] = np.sqrt(mu/rEar)      # (vx0 = 0, v_y0 = sqrt(mu/r)) AU/yr,
vxEar[0] = 0 - mu / (rEar**2)*(-0.5*dt) # corrected for time -dt/2

Your implementation of the force is wrong as it uses a constant radius and not the actual radius of the planet at the current time. Change to

#Implement Verlet Algorithm:
for k in range(0,N-1):
    rEar = (xEar[k]**2+yEar[k]**2)**0.5
    vxEar[k+1] = vxEar[k] - (mu*xEar[k]) / (rEar**3)*dt
    xEar [k+1] = xEar[k] + vxEar[k+1]*dt
    vyEar[k+1] = vyEar[k] - (mu*yEar[k]) / (rEar**3)*dt
    yEar [k+1] = yEar[k] + vyEar[k+1]*dt

edited Nov 19 '20 at 13:45

answered Apr 02 '20 at 09:15

Lutz Lehmann

25,219
2
22
51

Thanks for the help, I've been looking at some other orbital simulations and they seem to include the mass of the sun so I'm wondering if I've gone wrong by not including this? – 178971 May 04 '20 at 11:08
I have just reproduced your code. In the current simulation with one body in a central force field the mass of the sun is already included in the coefficient `mu=G*M`. The value is apparently chosen so that the period is 1 year at a distance of 1 AU. If you have a simulation with several interacting bodies you will need to explicitly scale the forces by the associated masses. – Lutz Lehmann May 04 '20 at 12:51
I am currently trying to do this by adding a calculation for acceleration for each body, using the masses and having a little difficulty. Would I put this within the for loop and create an array of zeros for acceleration also? I didn't think this would be correct as acceleration doesn't accumulate over time, but I'm struggling to see how to include it – 178971 May 04 '20 at 12:58
See https://stackoverflow.com/q/53813499/3088138 for a systematic approach, also my highest-voted answer https://stackoverflow.com/a/23671054/3088138 a gravity simulation using velocity Verlet in javascript. – Lutz Lehmann May 04 '20 at 13:16

Simulating Earth's orbit using symplectic Euler in Python

1 Answers1