1

When i set my timestep to dt = 3600*24 (1 day in secs) and run my simulation for 365 steps,I would expect Earth to make one full rotation around the sun. However the result is not even close to a quarter orbit. The number of steps I need to make one Earth orbit is more than I expect for the dt value I select, I could make dt bigger but it doesn't make sense in the units I am using.

class Planet:
    G = 6.67430e-11
    def __init__(self, name, position, velocity, mass):
        self.name = name
        self.position = position
        self.velocity = velocity
        self.mass = mass
        self.dr = np.array([0,0,0])
        self.dv = np.array([0,0,0])
        self.pos_data = [self.position.copy()]
        self.vel_data = [self.position.copy()]

    def acceleration(self, pos, planets):
        acceleration = np.array([0,0,0])
        for planet in planets:
            if planet != self:
                r = planet.position - pos
                r_mag = np.linalg.norm(r)
                force_magnitude = (self.G * self.mass * planet.mass) / (r_mag ** 2)
                accel_magnitude = force_magnitude / self.mass
                acceleration = acceleration + accel_magnitude * (r / r_mag)
        return acceleration

    def rk4(self, dt, planets):
        k1v = dt * self.acceleration(self.position, planets)
        k1r = dt * self.velocity        
        intermediate_position = self.position + 0.5 * k1r
        k2v = dt * self.acceleration(intermediate_position, planets)
        k2r = dt * (0.5 * k1v)        
        intermediate_position = self.position + 0.5 * k2r
        k3v = dt * self.acceleration(intermediate_position, planets)
        k3r = dt * (0.5 * k2v)        
        intermediate_position = self.position + k3r
        k4v = dt * self.acceleration(intermediate_position, planets)
        k4r = dt * (k3v)
        self.dp = (k1r + 2 * k2r + 2 * k3r + k4r) / 6
        self.dv = (k1v + 2 * k2v + 2 * k3v + k4v) / 6
    
    def update_state(self):
        self.position = self.position + self.dp
        self.velocity = self.velocity + self.dv
        self.pos_data.append(self.position)
        self.vel_data.append(self.velocity)

My usage of the class

sol     = Planet("Sun"    , np.array([0.0       , 0.0, 0.0]), np.array([0.0, 0.0    , 0.0]), 1.989e30)
mercury = Planet('Mercury', np.array([0.39e12   , 0.0, 0.0]), np.array([0.0, 47000.0, 0.0]), 3.285e23)
venus   = Planet('Venus'  , np.array([0.72e12   , 0.0, 0.0]), np.array([0.0, 35000.0, 0.0]), 4.867e24)
earth   = Planet('Earth'  , np.array([1.0e12    , 0.0, 0.0]), np.array([0.0, 29783.0, 0.0]), 5.972e24)
moon    = Planet('Moon'   , np.array([1.08e12    , 0.0, 0.0]), np.array([0.0, 30805.0, 0.0]), 7.348e22)

planets = [sol, earth, venus, moon, mercury]

dt = 3600*24 # 1 day
for _ in range(365*10):
    for planet in planets:
        planet.rk4(dt, planets)
    for planet in planets:
        planet.update_state()

plt.clf()
for planet in planets:
    x = [pos[0] for pos in planet.pos_data]
    y = [pos[1] for pos in planet.pos_data]
    z = [pos[2] for pos in planet.pos_data]
    plt.scatter(x[-1], y[-1])
    plt.plot(x,y)
plt.scatter(0,0)
plt.show()

enter image description here

Maybe I'm wrong and shouldn't expect it to be exact, but I thought it at least should be close. Any helpp would be great. Thanks.

This is my first ever question, sorry if I format the question wrong.

Fergie
  • 11
  • 3
  • I would think there would be a much better approximation in [Kepler's laws](https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion), which are specifically observed planetary motion? You are going to run into problems with the [precession of Mercury](https://en.wikipedia.org/wiki/Tests_of_general_relativity). – Neil Aug 11 '23 at 16:04
  • 2
    Isn't the precession of Mercury an issue for all non-relativistic approaches, including Kepler's laws? – slothrop Aug 11 '23 at 16:06
  • 2
    The distances look off in the initial setup. The Earth-Sun distance is more like 1.5e11 metres than 1.0e12. – slothrop Aug 11 '23 at 16:11
  • @slothrop Exactly; you would have to correct for gr. This is a very famous problem that led to _eg_ [Vulcan](https://en.wikipedia.org/wiki/Vulcan_(hypothetical_planet)) being proposed. – Neil Aug 11 '23 at 16:26

1 Answers1

1

I'm not an expert on this, so it's possible the code can be improved further. But the main issues are:

  1. The calculation of k2r, k3r and k4r should incorporate the updated velocity, not the size of the update to the velocity.

So instead of:

k2r = dt * (0.5 * k1v) 

k3r = dt * (0.5 * k2v) 

k4r = dt * (k3v)

it should be:

k2r = dt * (self.velocity + 0.5 * k1v)

k3r = dt * (self.velocity + 0.5 * k2v)

k4r = dt * (self.velocity + k3v)
  1. The initial distances are off. The Earth-Sun distance (i.e. 1 astronomical unit) is around 1.50e11 metres. The 1e12 in the current code is off by a factor of 6.67 from that. The other planetary distances seem to be off by the same proportional amount, while the Earth-Moon distance is considerably off (it should be about 0.00266 AU or 3.84e8 metres).

Correcting the setup code to:

au = 1.50e11

sol     = Planet("Sun"    , np.array([0.0       , 0.0, 0.0]), np.array([0.0, 0.0    , 0.0]), 1.989e30)
mercury = Planet('Mercury', np.array([0.39*au   , 0.0, 0.0]), np.array([0.0, 47000.0, 0.0]), 3.285e23)
venus   = Planet('Venus'  , np.array([0.72*au   , 0.0, 0.0]), np.array([0.0, 35000.0, 0.0]), 4.867e24)
earth   = Planet('Earth'  , np.array([1.0*au    , 0.0, 0.0]), np.array([0.0, 29783.0, 0.0]), 5.972e24)
moon    = Planet('Moon'   , np.array([1.00266*au , 0.0, 0.0]), np.array([0.0, 30805.0, 0.0]), 7.348e22)

gives much more realistic-looking results (though I'm not sure about the behaviour of the Moon - you might need a shorter timestep to get that to work well).

slothrop
  • 3,218
  • 1
  • 18
  • 11
  • k2r = dt * (self.velocity + 0.5 * k1v) – Fergie Aug 16 '23 at 09:37
  • `k2r = dt * (self.velocity + 0.5 * k1v)` is indeed correct, as said in the answer. The problem is that your current code instead has `k2r = dt * (0.5 * k1v)`. – slothrop Aug 16 '23 at 09:39