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.scatter( x, y )
    pylab.axis([-200, 200, -200, 200])

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?)

  • $\begingroup$ Look at answers to this question as you will need a somewhat specialized time advancement technique for this problem. $\endgroup$ May 15, 2013 at 15:24
  • 3
    $\begingroup$ 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. $\endgroup$
    – J. M.
    May 15, 2013 at 15:36

2 Answers 2


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


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.

  • 3
    $\begingroup$ 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. $\endgroup$ May 17, 2013 at 3:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.