and I compared it to simple verlet integration (and euler, but it has very very poor performance). It doesnt seem that RK4 with constant step is better than verlet. RK4 with adaptive step is better, but not so much. I wonder if I'm doing something wrong? Or in which sense it is said that RK4 is much superior to verlet?

The think is that Force is evaluated 4x per RK4 step, but only 1x per verlet step. So to get same performance I can set time_step 4x smaller for verlet. With 4x smaller time step verlet is more precise than RK4 with constant step and almost comparable with RK4 with addaptive step.

3 Answers
3

I don't know if I'm going to answer your specific questions, but here are my thoughts.

You have defined a very simple force model. In this case, saving some steps may not improve the performance, because calculating the new step in RK4 may take longer. If the force model is more complex, RK4 with adaptive step may save you much time. From your plot I think that the Verlet is also diverging from the correct solution, a repeating ellipse.

I know this question is quite old by now, but this really has nothing to do with the "superiority" of one of these methods over the other, or with your programming of them - they're just good at different things. (So no, this answer won't really be about code. Or even about programming. It's more about math, really...)

The Runge-Kutta family of solvers are pretty good at attacking almost any problem with quite good precision and, in case of the adaptive methods, performance. However, they are not symplectic, which shortly means that they do not conserve energy in the problem.

The Verlet method on the other hand, may require a much smaller step size than the RK methods in order to minimize oscillations in the solution, but the method is symplectic.

Your problem is energy-conserving; after an arbitrary number of revolutions, the planetary body's total energy (kinetic + potential) should be the same as it was with the initial conditions. This will be true with a Verlet integrator (at least as a time-windowed average), but it will not with an integrator from the RK family - over time, the RK solvers will build up an error that is not reduced, due to energy lost in the numerical integration.

To verify this, try saving the total energy in the system at each time step, and plot it (you might have to do much more than ten revolutions to notice the difference). The RK method will have steadily decreasing energy, while the Verlet method will give an oscillation around a constant energy.

If this is the exact problem you need to solve, I also recommend the Kepler equations, which solve this problem analytically. Even for rather complex systems with many planets, moons etc, the interplanetary interactions are so insignificant that you can use the Kepler equations for each body and it's rotational center individually without much loss of precision. However, if you're writing a game you might really be interested in some other problem, and this was just an example to learn - in that case, read up on the properties of the various solver families, and choose one that is appropriate for your problems.

OK, finaly, I used adaptive Runge-Kutta-Fehlberg (RKF45).
Interestingly, it is faster (less step is needed) when I ask for higher precission ( optimum is 1E-9 ) because at lower precision ( <1e-6 ) the solution is unstable, and lot of iterations are wasted by discarted steps ( the steps wich were to long and unprecise ). When I ask for even better precission ( 1E-12 ) it requres more steps because time steps are shorter.

For circular orbits precission can be decreesed up to ( 1e-5 ) with up to 3x speed gain, however while I need to simulated highly eccentric (eliptical orbits) I will prefer keep secure 1E-9 precission.