I've come across a stiff equation in solving the Circular Restricted Three Body Problem. [An object is moving considering the effect of the gravitational forces caused by two gravitational sources fixed in a 2D Space.]

The equations are these:

$x''=-\frac{GM_1 (x-x_1)}{\sqrt{(x-x_1)^2+y^2}^3}--\frac{GM_2 (x-x_2)}{\sqrt{(x-x_2)^2+y^2}^3}$

$y''=-\frac{GM_1 y}{\sqrt{(x-x_1)^2+y^2}^3}--\frac{GM_2 y}{\sqrt{(x-x_2)^2+y^2}^3}$

Neither Euler Method or Runge Kutta will work as the property near $(x_1, 0)$ or $(x_1, 0)$ is not good. The derivatives change too fast. The simulation can't be solved out right. The object is too easy to hit on the gravitational source.

How can I fix this?

Thank you!

  • 3
    $\begingroup$ What do you mean by "Neither Euler Method or Runge Kutta will work"? Are these your own implementations of these methods? Are you using a fixed time step size? I believe that both of these methods should work just fine even close to the singularities if you choose an adaptive time step size. $\endgroup$ Commented Jan 16, 2013 at 14:58
  • $\begingroup$ Have you tried an implicit method, such as Backward Euler? $\endgroup$
    – Paul
    Commented May 25, 2013 at 1:52

2 Answers 2


Conventional integrators do not preserve the "shape" of phase space, leading to systematic energy gain or loss, thus you should consider "sympelectic" integrators. For close encounters, you should consider the method of Chambers (1999) A hybrid symplectic integrator that permits close encounters between massive bodies.


You can find a very readable introduction to modern ODE solvers in chapter 7 of Cleve Moler's book, Numerical Computing With MATLAB, available online here: http://www.mathworks.com/moler/odes.pdf Among the topics he discusses are stability, how to obtain a prescribed accuracy with variable time step algorithms, and application to the two-body problem.


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.