Runge–Kutta RK4 not better than Verlet?

Question

I'm just testing several integration schemes for orbital dynamics in game. I took RK4 with constant and adaptive step here http://www.physics.buffalo.edu/phy410-505/2011/topic2/app1/index.html

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.

See the image: https://lh4.googleusercontent.com/-I4wWQYV6o4g/UW5pK93WPVI/AAAAAAAAA7I/PHSsp2nEjx0/s800/kepler.png

10T means 10 orbital periods, the following number 48968,7920,48966 is number of force evaluations needed

the python code (using pylab) is following:

from pylab import * 
import math

G_m1_plus_m2 = 4 * math.pi**2

ForceEvals = 0

def getForce(x,y):
    global ForceEvals
    ForceEvals += 1
    r = math.sqrt( x**2 + y**2 )
    A = - G_m1_plus_m2 / r**3
    return x*A,y*A

def equations(trv):
    x  = trv[0]; y  = trv[1]; vx = trv[2]; vy = trv[3];
    ax,ay = getForce(x,y)
    flow = array([ vx, vy, ax, ay ])
    return flow

def SimpleStep( x, dt, flow ):
    x += dt*flow(x)

def verletStep1( x, dt, flow ):
    ax,ay = getForce(x[0],x[1])
    vx   = x[2] + dt*ax; vy   = x[3] + dt*ay; 
    x[0]+= vx*dt;        x[1]+= vy*dt;
    x[2] = vx;        x[3] = vy;

def RK4_step(x, dt, flow):    # replaces x(t) by x(t + dt)
    k1 = dt * flow(x);     
    x_temp = x + k1 / 2;   k2 = dt * flow(x_temp)
    x_temp = x + k2 / 2;   k3 = dt * flow(x_temp)
    x_temp = x + k3    ;   k4 = dt * flow(x_temp)
    x += (k1 + 2*k2 + 2*k3 + k4) / 6

def RK4_adaptive_step(x, dt, flow, accuracy=1e-6):  # from Numerical Recipes
    SAFETY = 0.9; PGROW = -0.2; PSHRINK = -0.25;  ERRCON = 1.89E-4; TINY = 1.0E-30
    scale = flow(x)
    scale = abs(x) + abs(scale * dt) + TINY
    while True:
        x_half = x.copy();  RK4_step(x_half, dt/2, flow); RK4_step(x_half, dt/2, flow)
        x_full = x.copy();  RK4_step(x_full, dt  , flow)
        Delta = (x_half - x_full)
        error = max( abs(Delta[:] / scale[:]) ) / accuracy
        if error <= 1:
            break;
        dt_temp = SAFETY * dt * error**PSHRINK
        if dt >= 0:
            dt = max(dt_temp, 0.1 * dt)
        else:
            dt = min(dt_temp, 0.1 * dt)
        if abs(dt) == 0.0:
            raise OverflowError("step size underflow")
    if error > ERRCON:
        dt *= SAFETY * error**PGROW
    else:
        dt *= 5
    x[:] = x_half[:] + Delta[:] / 15
    return dt    

def integrate( trv0, dt, F, t_max, method='RK4', accuracy=1e-6 ):
    global ForceEvals
    ForceEvals = 0
    trv = trv0.copy()
    step = 0
    t = 0
    print "integrating with method: ",method," ... "
    while True:
        if method=='RK4adaptive':
            dt = RK4_adaptive_step(trv, dt, equations, accuracy)
        elif method=='RK4':
            RK4_step(trv, dt, equations)
        elif method=='Euler':
            SimpleStep(trv, dt, equations)
        elif method=='Verlet':
            verletStep1(trv, dt, equations)
        step += 1
        t+=dt
        F[:,step] = trv[:]
        if t > t_max:
            break
    print " step = ", step


# ============ MAIN PROGRAM BODY =========================

r_aphelion   = 1
eccentricity = 0.95
a = r_aphelion / (1 + eccentricity)
T = a**1.5
vy0 = math.sqrt(G_m1_plus_m2 * (2 / r_aphelion - 1 / a))
print " Semimajor axis a = ", a, " AU"
print " Period T = ", T, " yr"
print " v_y(0) = ", vy0, " AU/yr"
dt       = 0.0003
accuracy = 0.0001

#                 x        y     vx  vy
trv0 = array([ r_aphelion, 0,    0, vy0 ])             

def testMethod( trv0, dt, fT, n, method, style ):
    print " "
    F = zeros((4,n));
    integrate(trv0, dt, F, T*fT, method, accuracy);
    print "Periods ",fT," ForceEvals ",  ForceEvals
    plot(F[0],F[1], style ,label=method+" "+str(fT)+"T "+str(  ForceEvals ) ); 

testMethod( trv0, dt, 10, 20000  , 'RK4', '-' )
testMethod( trv0, dt, 10, 10000  , 'RK4adaptive', 'o-' )
testMethod( trv0, dt/4, 10, 100000, 'Verlet', '-' )
#testMethod( trv0, dt/160, 2, 1000000, 'Euler', '-' )

legend();
axis("equal")
savefig("kepler.png")
show();

Answer 1

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.

For orbital mechanics you can also try the RK7(8) adaptive integrator, the Adams-Bashforth multistep or the Gauss Jackson method. Here is a paper showing some of these methods: http://drum.lib.umd.edu/bitstream/1903/2202/7/2004-berry-healy-jas.pdf

And finally, if your force model will be always a simple central force, like in this example, take a look to the Kepler equation. Solving it is precise, fast and you can jump to an arbitrary time.

Answer 2

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.

Answer 3

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.

There is the code if somebody will solve similar problem http://www.openprocessing.org/sketch/96977 it also shows how many force evaluations it requires to simulate one time unit of trajectrory length