15

I'm creating a simple astronomy simulator that should use Newtonian physics to simulate movement of planets in a system (or any objects, for that matter). All the bodies are circles in an Euclidean plane, that have properties such as position, velocity, mass, radius and the resultant force.

I want to update the universe in small time steps, usually a few milliseconds, but I'm not sure how to correctly calculate the changes in position.

The force is simple: fr = sum(G * body.m * bodyi.m / dist(body, bodyi)^2).

But how do I go on from there?

I could do this:

a = Fr/body.m
v += a*dt
position += v*dt

But that would, of course, be false. Maybe if I added 0.5 as a factor in position calculation?

jcora
  • 251
  • 1
  • 5

1 Answers1

17

You've essentially got the answer - no need for the factor of 0.5.

Essentially you have a two-dimensional system of first-order ODEs: \begin{align} \dot{x} & = v \\ \dot{v} & = \frac{F}{m}, \end{align} where everything is a function of time except presumably $m$, and dots denote time-derivatives. If you do a simple, forward-Euler-esque first-order differencing of these, you find \begin{align} \frac{x_{n+1}-x_n}{\Delta t} & = v_n \\ \frac{v_{n+1}-v_n}{\Delta t} & = \frac{F_n}{m}, \end{align} or \begin{align} x_{n+1} & = x_n + \Delta t \cdot v_n \\ v_{n+1} & = v_n + \Delta t \frac{F_n}{m}. \end{align} Here I'm indexing the timestep with $n$.

However, forward-Euler is inherently unstable. Fortunately, there's a symplectic method right around the corner. (That linked article is more of a stub, but it might contain some useful links.) The key is to advance positions from $t_n$ to $t_{n+1}$ using velocities at $t_{n+1/2}$. That is, suppose you were given $x_0$ and $v_{1/2}$ for each particle. Then you could use \begin{align} x_{n+1} & = x_n + \Delta t \cdot v_{n+1/2} \\ v_{n+1/2} & = v_{n-1/2} + \Delta t \frac{F_n}{m} \end{align} to integrate forward in time. This is known as the leapfrog method. With this, your system conserves an energy of a sort, and orbits are less likely to fly off to infinity or some such thing due to exponential growth of roundoff error.

The only catch is how to get $v_{1/2}$, since presumably you start with $v_0$. There you should just use as accurate a scheme as you are up to writing, fourth-order Runge-Kutta being a popular choice. It may be long-term unstable, but there's only so much error you'll introduce in half a timestep, and that error will be kept small afterward by the leapfrog scheme.

Finally, this answer applies to any general Newtonian gravity sim. If you really want perfect circles, as mentioned in passing in the question, then you won't get those except in an idealized system in which planets do not interact with each other and the initial conditions are chosen just right. If that's the case, then, you don't need to integrate at all, since the angular velocity (radians per unit time) of such an object is simply $$ \omega = \sqrt{\frac{GM}{r^3}}, $$ where $M$ is the mass of the central object and $r$ is the radius of the orbit. This can be used to test the accuracy of your simulation.

  • Hey, can you explain why don't I need the 0.5 factor? It seems to be doing the same thing as taking the speed n-1/2dt seconds ago, which is what it seems you are suggesting. – jcora Oct 27 '12 at 23:16
  • That would only be the same if the speed at $(n-1)$ was 0. What you want is a first-order estimate of the mean of $v_n$ and $v_{n+1}$ (the latter you don't know), not the mean of $v_n$ and $0$. –  Oct 28 '12 at 03:02