Earth and Moon computer simulation

Question

Edit: I made a test model with 2 bodies. The sun and earth or w/e. But I now have the problem of zero divisions every time that one of the coordinates is zero.

import numpy as np
from numpy import exp, sqrt, pi, linspace
from time import clock,time
import pylab as pylab

S1 = 1E3                   #Mass Star
M1 = 1                     #Mass Planet


Rs = np.array([0, 0])      #position
R1 = np.array([0,200])

Vs = np.array([0,0])       #velocity
V1 = np.array([10,0])

As = np.array([0,0])       #acceleration
A1 = np.array([0,0])


dt = 0.1                   #timestep
N = 1000                   #steps
G = 1                       #G const

for i in range(1,N):

    deltaR = Rs - R1
    Rsq = deltaR**2
    F12 = -G*M1*S1/(Rsq) 
    F21 = -F12

    As = F21/S1
    A1 = F12/M1

    Vs += As*dt
    V1 += A1*dt

    Rs += Vs*dt + 0.5*As*dt**2
    R1 += V1*dt + 0.5*A1*dt**2

    x = [ Rs[0], R1[0]]
    y = [ Rs[0], R1[0]]
    pylab.ion()
    pylab.scatter( x, y )
    pylab.axis([-200, 200, -200, 200])
    pylab.plot
    pylab.clf()

Any tips?

So I want to simulate the solar system but want to start simple with one orbiting body. However, I never did anything like this before and was wondering if anyone here could give me some hints.

Clearly the first thing we need to use is the law of gravitation and Kepler. Then I'd like to throw everything into an array and just calculate the timesteps. However, I am not sure on how to get the position of the planet/moon as a function of time. (Can I just integrate the acceleration via the law of gravitation?)

Look at answers to this question as you will need a somewhat specialized time advancement technique for this problem. — Godric Seer, May 15 '13 at 15:24
To summarize: you're dealing with the $n$-body problem, $n=2$. The DEs are in that article. As @Godric hints, you'll need a symplectic integrator for your Hamiltonian system. — J. M., May 15 '13 at 15:36

score 3 · Answer 1 · edited May 16 '13 at 15:42

3

For this type of problem, the leapfrog time-integration method is a good choice, see details here.

edited May 16 '13 at 15:42

J. M.

3,155
28
37

answered May 16 '13 at 10:56

Maxim Umansky

2,525
1
12
15

Max Hutchinson · Answer 2 · 2013-05-16T15:44:49.203

1

If you want to start simple, you can just use simple kinematics: \begin{align*} a(t) &= \text{Gravity}[x(t)] \\ v(t+ \Delta t) &= v(t) + a(t) \Delta t \\ x(t+\Delta t) &= x(t) + v(t) \Delta t + a(t) (\Delta t)^2/2 \\ \end{align*} That should give you something that looks OK for small time-steps, $\Delta t$, and short duration. I also think it helps demonstrate the connection between continuous physics and explicit methods.

For all kinds of reasons, simple methods like these don't work for long durations. For a quantitatively valid scheme, check out the question's comments or a different answer.

edited May 16 '13 at 15:44

answered May 15 '13 at 16:59

Max Hutchinson

3,051
16
29

3

This simplest of all possible quadratures fails to conserve energy and won't generate closed orbits. The class of integrators that work properly are said to be symplectic. – dmckee --- ex-moderator kitten May 17 '13 at 03:00

Earth and Moon computer simulation

2 Answers2