Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

This question is more about how to tackle a problem numerically.

In a small project I wanted to simulate the coorbital motion of Janus and Epimetheus. This is basically a three body problem. I choose Saturn to be fixed at the origin, let $r_1$ and $r_2$ be the location vectors of janus and epimetheus, respectively. Since the effect occurs when Janus and Epimetheus are very close together I picked relative coordinates for a better resolution, i.e. $r=r_1-r_2$ and $R=r_1+r_2$. Now I get the following equations of motion:

$$ \frac {d^2}{dt^2} \binom{R}{r} = - G (m_2\pm m_1) \frac R {R^3} - 4 M G \left(\frac {r+R}{(r+R)^3} \mp \frac {r-R}{(r-R)^3}\right ) $$

where $m_i$ corresponds to the masses of the moons, $M$ is the mass of Saturn and $G$ the gravitational constant. The problem arises when I try to solve this numerically. One has to deal with values of completely different magnitudes, i.e. $M \sim e^{28}$ and $m_i \sim e^{17}$. And $r$, $R$ are in the regions of 0 to 150,000.

To be honest I am not sure if this is the place forum to discuss such numerical problems.

More Information:

Code is written in Matlab and I use a standard ODE solver to obtain the result. However this is breaking down because the step size cannot be reduced under machine precision. (I find this not to be surprising because one has to deal with the already mentioned orders of magnitude).

share|improve this question

migrated from Mar 29 '12 at 19:41

This question came from our site for active researchers, academics and students of physics.

Are you running this simulation in SI units? At the minimum, you should divide everything through by some factor of $G*m_{2}$, so that you can eliminate a few orders of magnitude. – Jerry Schirmer Nov 30 '11 at 14:58
Hi, I this, but its still not working... The same problems occur as before. :( – bios Nov 30 '11 at 16:49
You have to set your unit of mass to one of the moon's masses, and your units of length/time to set things to 1. Nothing should be smaller than 1/100 if you write it well. There is no need for an over-the-counter solver. Write code to do this yourself, where you control the stepsize. Breakdowns in stepsize with these types of potentials can occur at collisions, where the solver will try to reduce the stepsize until convergence, and at collision there is no convergence. You need to make sure that the orbits are not collinear, so you need to view the simulation. You can't get an answer as is. – Ron Maimon Nov 30 '11 at 18:00
Please avoid abbreviations in title. DGL = Differentialgleichung? – Qmechanic Jan 30 '12 at 23:01
What standard ODE solver are you using? – Geoff Oxberry Mar 29 '12 at 19:43

Your current approach ruins numerical stability; in fact you probably lose resolution this way.

Take as coordinates for each satellite its Kepler variables and the angle of the plane containing the position of the satellite, the velocity, and the origin. The differential equations in the absence of interaction between the satellites are then trivially simple, and only the interaction becomes somewhat complicated. As the interaction is tiny if the satellites are far away, the resulting dynamics should be numerically stable.

share|improve this answer

Instead of using a "classical" (stiff) ODE solver, you can use dedicated algorithms for geometric numerical integration. See for example this book and the GNI codes you can find on Ernst Hairer's website.

share|improve this answer

How about if you have three steps in your simulation:

  1. update Janus position by calculating the Janus - Saturn force.
  2. update Epimetheus position by calculating the Epimetheus - Saturn force.
  3. update Janus and Epimetheus position by calculating the Janus - Epimetheus force.

Possibly using finer timesteps for #3.

I'm not sure if this will help. I suppose the real problem is that the magnitude of the force is different in the moon - moon and the moon - Saturn case, except if the moons are close?


  1. if the moons close, compute an approximate moons - Saturn force using their center of mass vector and update both positions with the same vector.
  2. if they are far apart, update them separately.
  3. as before.

Best of luck!

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.