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: 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](http://en.wikipedia.org/wiki/Symplectic_integrator), 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.

