8

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!

J. M.
  • 3,155
  • 28
  • 37
Fin_de_Siecle
  • 81
  • 2
  • 3
    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 01:52

2 Answers2

7

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.

Jed Brown
  • 25,650
  • 3
  • 72
  • 130
0

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.

Bill Greene
  • 6,064
  • 1
  • 16
  • 24