9
$\begingroup$

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

Improve this question
$\endgroup$
5
  • 2
    $\begingroup$ 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. $\endgroup$
    – Jerry Schirmer
    Nov 30, 2011 at 14:58
  • $\begingroup$ Hi, I this, but its still not working... The same problems occur as before. :( $\endgroup$
    – bios
    Nov 30, 2011 at 16:49
  • $\begingroup$ 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. $\endgroup$
    – Ron Maimon
    Nov 30, 2011 at 18:00
  • 1
    $\begingroup$ Please avoid abbreviations in title. DGL = Differentialgleichung? $\endgroup$
    – Qmechanic
    Jan 30, 2012 at 23:01
  • $\begingroup$ What standard ODE solver are you using? $\endgroup$
    – Geoff Oxberry
    Mar 29, 2012 at 19:43

3 Answers 3

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
2
$\begingroup$

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.

Improve this answer
$\endgroup$
0
$\begingroup$

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?

Alternatively:

  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!

Improve this answer
$\endgroup$

Your Answer

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