14
$\begingroup$

I would like to use Runge-Kutta 8th order method (89) in a celestial mechanics / astrodynamics application, written in C++, using a Windows machine. Therefore I wonder if anyone knows a good library / implementation that is documented and free to use ? It is ok if it is written in C, as long as there aren't any compilation problems to be expected.

So far I have found this library ( mymathlib ). The code seems ok, but I haven't found any information on licensing.

Can you help me by revealing some of the alternatives you might know and would suit my problem ?

EDIT :
I see that there aren't really that many C/C++ source codes available as I expected. Therefore a Matlab/Octave version would be ok too (still has to be free to use).

Improve this question
$\endgroup$

4 Answers 4

Reset to default
10
$\begingroup$

Both GNU Scientific Library (GSL) (C) and Boost Odeint (C++) feature 8th order Runge-Kutta methods.

Both are opensource, and under linux and mac they should be directly available from the package manager. Under windows, it will probably be easier for you to use Boost rather than GSL.

GSL is published under the GPL license, and Boost Odeint under the Boost license.

Edit: Ok, Boost Odeint does NOT have the Runge-Kutta 89 method, only the 78, but it does provide a recipe for making arbitrary Runge-Kutta steppers.

8th order methods are quite high however, and most likely overkill for your problem.

Prince-Dormand refers to a specific kind of Runge-Kutta, and is not directly related to the order but the most common is 45. Matlabs ode45, which is their recommended ODE algorithm is a Prince-Dormand 45 implementation. This is the same algorithm as implemented in Boost Odeint Runge_Kutta_Dopri5.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Thank you for answer. OK that is embarrassing now, I have taken a look at Boost Odeint even before asking here, and only found "runge_kutta_fehlberg78". Is this the right thing ? Actually I don't know the differences between methonds when used in practice, but I was looking for a RK89 ( called also Dormand-Prince as I search the internet ). Can you comment or expand your answer regarding this matter please ? Thank you. $\endgroup$
    – James C
    Commented Aug 17, 2014 at 13:58
  • $\begingroup$ Updated post to answer your questions. Prince-Dormand 45 will most likely solve your problems nicely. $\endgroup$
    – LKlevin
    Commented Aug 17, 2014 at 18:22
17
$\begingroup$

If you're doing celestial mechanics over long time scales, using a classical Runge-Kutta integrator will not preserve energy. In that case, using a symplectic integrator would probably be better. Boost.odeint also implements a 4th-order symplectic Runge-Kutta scheme that would work better for long time intervals. GSL does not implement any symplectic methods, as far as I can tell.

Improve this answer
$\endgroup$
7
  • $\begingroup$ Thank you for answer. Would a 4th-order symplectic Runge-Kutta give better results than RKF78, if used with Earth satellites ( low orbit as well as deeper space orbit ), perhaps over a period of 1-3 orbits ? $\endgroup$
    – James C
    Commented Aug 18, 2014 at 11:30
  • $\begingroup$ @JamesC Yes. In a long period, symplectic method is much better. $\endgroup$
    – eccstartup
    Commented Aug 18, 2014 at 13:27
  • $\begingroup$ @eccstartup - What would you consider a long period here ? Because it could be as long as one orbit of a planet around Sun, or a few orbits of a weather satellite around Earth etc.. $\endgroup$
    – James C
    Commented Aug 18, 2014 at 13:29
  • $\begingroup$ @JamesC I have not observed that big problem. But for my model problems, with many orbits calculated, symplectic methods give very perfect orbits. $\endgroup$
    – eccstartup
    Commented Aug 18, 2014 at 13:35
  • $\begingroup$ So, it is an advice to program you own a version of implicit Runge-Kutta method, which includes a lot of symplectic methods with as higher order as you want. $\endgroup$
    – eccstartup
    Commented Aug 18, 2014 at 13:36
6
$\begingroup$

summarizing some points:

  1. If it's a long-term integration of a non-dissapative model, a symplectic integrator is what you're looking for.
  2. Otherwise, since it's an equation of motion, Runge-Kutta Nystrom methods will be more efficient than a transformation to a first order system. There are high order RKN methods due to DP. There are some implementations, like here in Julia they are documented and here's a MATLAB one.
  3. High order Runge-Kutta methods are only needed if you want a high accuracy solution. If it's lower tolerances then a 5th order RK will likely be faster (for the same error). The best thing to do if you need to solve this often is to test a bunch of different methods. In this set of benchmarks on 3-body problems we see that (for the same error) high order RK methods are only a marginal improvement in speed, though as error -> 0 you can see that the improvement already goes to >5x against Dormand-Prince 45 (DP5) when you're looking at 4 digits of accuracy (tolerances are wayyyy lower for this though. Tolerances are only a ballpark in any problem). As you pull the tolerances even lower the improvement from a high order RK method grows, but you may need to start using higher precision numbers.
  4. Dormand-Prince order 7/8 algorithm has a different 8th order tableau than the DP853 method of Hairer's dop853 and DifferentialEquations.jl's DP8 (which are the same). The latter 853 method cannot be implemented in standard tableau version of an Runge-Kutta method since it's error estimator is non-standard. But this method is much more efficient and I wouldn't recommend even using the older Fehlberg 7/8 or DP 7/8 methods.
  5. For high order RK methods, the Verner "Efficient" methods are the gold standard. That shows up in the benchmarks I linked. You can code those into Boost yourself, or use one of the 2 packages which implement them if you want those easier (Mathematica or DifferentialEquations.jl).
Improve this answer
$\endgroup$
3
$\begingroup$

I would like to add that while what Geoff Oxberry suggests for long-term integration (using symplectic integrators) is true, in some cases it won't work. More specifically, if you have dissipative forces, your system does not preserve the energy anymore, and therefore you cannot resort to symplectic integrators in that case. The person asking the question was talking about low Earth orbits, and such orbits exhibit a big amount of atmospheric drag, that is a dissipative force that precludes using such symplectic integrators.

In that specific case (and for cases where you cannot use /do not have access to/do not want to use symplectic integrators), I would recommend the use of Bulirsch-Stoer integrator if you need precision and efficiency over long timescales. It works well by experience, and is also recommended by the Numerical Recipes (Press et al., 2007).

Improve this answer
$\endgroup$
6
  • $\begingroup$ No, do not recommend Numerical Recipes. Especially, in most cases Burlirsch-Stoer shouldn't be recommended. This is a well-known problem with the book. See a bunch of rebuttles from top researchers in the field here: uwyo.edu/buerkle/misc/wnotnr.html. If you want benchmarks on this, see Hairer's first book where you'll see BS almost never does well. Higher order is only more efficient when errors are low enough, and we (and others) have done benchmarking to show quite consistently that it's only efficient for sub-floating point precision. $\endgroup$
    – Chris Rackauckas
    Commented Dec 7, 2017 at 15:42
  • $\begingroup$ I can't speak for NR too much since I used it mostly for ODEs, but it seems to me that the complaints on the page you link to are old and have been adressed by NR's authors in their response (end of the page), but this is off topic. Concerning the long-term integration of orbits with high precision (say, 13-14 digits) which is what I was mentioning in my answer, it is proven since long that extrapolation methods work well (see Montenbruck & Gill's chapter on Numerical Integration). More recent papers use it too, and it has proven to me and others a reliable and efficient method. $\endgroup$
    – viiv
    Commented Dec 7, 2017 at 17:39
  • $\begingroup$ M&G only test it against dop853 and more modern high order RK methods, like those due to Verner, are much more efficient. M&G also only seems to measure using function evaluations, which are a weak indicator of timing. It also doesn't time against Runge-Kutta Nystrom methods which are specifically for 2nd order ODEs and are more efficient than the first-order RK methods applied to second order by quite a bit. At 13-14 digits BS is probably competitive on most problems, but it's far from the obvious choice and I haven't seen a work-precision diagram with recent methods disagree with that. $\endgroup$
    – Chris Rackauckas
    Commented Dec 7, 2017 at 18:39
  • $\begingroup$ M&G do test RKN against RK, and BS and others against RKN (pages 123-132 and 151-154) and say they are the most efficient of the RK methods (not including Verner even though they cite him). BS is showed to be efficient at 13-14 digits, which was my claim, I have seen it tested against dop853, ABM(12), Taylor, and standard RK8 and it performs well. I have to admit I haven't seen it tested againt RKN but from what I can see from M&G it is not far from FILG11 for instance. I'm genuinely interested in Verner's RK and will look at your links above. Do you have a paper that tests all of them to see? $\endgroup$
    – viiv
    Commented Dec 7, 2017 at 19:32
  • $\begingroup$ I went back and re-ran a bunch of benchmarks at DiffEqBenchmarks.jl and odex doesn't tend to fair well. So at least for 1st order ODEs and for tolerances >=1e-13, extrapolation doesn't seem to do well and it's usually not even close. This is in line with the claim above. $\endgroup$
    – Chris Rackauckas
    Commented Dec 9, 2017 at 14:45

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.