Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. Join them; it only takes a minute:

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

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!

share|improve this question

migrated from Jan 16 '13 at 11:23

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

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. – Wolfgang Bangerth Jan 16 '13 at 14:58
Have you tried an implicit method, such as Backward Euler? – Paul May 25 '13 at 1:52

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.

share|improve this answer

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

share|improve this answer

Your Answer


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

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