Runge Kutta RK4 is not better than Verlet?

Question

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

and I compared it to the simple integration of verbal elements (and the Euler, but it has very low performance). RK4 with a constant pitch seems to be better than willow. RK4 with adaptive pitch is better, but not so much. I wonder what I'm doing something wrong? Or in what sense is it said that the RK4 is far superior to the camel?

Think of Force being evaluated 4 times per RK4 step, but only 1 time per step to the top. Thus, to get the same performance, I can set time_step 4x less for willow. With a 4-fold time step, the camel is more accurate than RK4 with a constant step and is almost comparable to RK4 with an additional step.

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

10T means 10 periods of the orbit, the next number 48968.7920.48966 is the number of required strength estimates

python code (using pylab) is as follows:

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 # xy 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 do not know if I will answer your specific questions, but here are my thoughts.

You have defined a very simple model of strength. In this case, saving some steps may not improve performance, because it may take longer to calculate a new step in RK4. If the strength model is more complex, the adaptive pitch RK4 can save you a lot of time. From your plot, I think Verlet also deviates from the correct solution, a repeating ellipse.

For orbital mechanics, you can also try the adaptive integrator RK7 (8), the multi-stage Adams-Bashfort method, or the Gauss Jackson method. Here is a document 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 is always a simple central force, as in this example, take a look at the Kepler equation. The solution is accurate, fast, and you can go at any time.

Answer 2

I know this question is quite old, but it really has nothing to do with the "superiority" of one of these methods over another or with their programming on them - they are just good in different things. (So ​​no, this answer will not really be about code or even about programming. It's more about math, really ...)

The Runge-Kutta family of solvers pretty well attacks almost any problem with pretty good accuracy and, in the case of adaptive methods, performance. However, they are not symplectic , which soon means that they do not save energy in this problem.

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

Your problem is energy conservation; after an arbitrary number of revolutions, the total energy of the planetary body (kinetic + potential) should be the same as under the initial conditions. This will be true with the Verlet integrator (at least, as a time average), but not with the integrator from the RK family - over time, RK developers will create an error that will not be reduced due to the energy lost during numerical integration.

To make sure of this, try to save the full energy in the system at each time step and draw it (you may need much more than ten revolutions to notice the difference). The RK method will steadily reduce energy, and the Verlet method will oscillate around constant energy.

If this is the exact problem you need to solve, I also recommend the Kepler equations, which analytically solve this problem. Even for fairly complex systems with many planets, moons, etc. Interplanetary interactions are so negligible that you can use the Kepler equations for each body and its rotational center individually without significant loss of accuracy. However, if you are writing a game, you might be interested in some other problem, and this was just an example to study - in this case, get acquainted with the characteristics of different families of solvers and choose the one that suits your problems.

Answer 3

OK, finally, I used the adaptive Runge-Kutta-Felberg (RKF45). It is interesting that it is faster (a smaller step is required) when I ask for a higher dose (optimal - 1E-9), because with less accuracy (<1e-6) the solution is unstable, and many iterations are lost as a result of the steps dropped (steps that should have lasted and inaccurate). When I request an even better precession (1E-12), it asks for more steps because the time steps are shorter.

For circular orbits, the precession can be reduced to (1e-5) with a gain of up to 3x, however, although I need to model very eccentric (elliptical orbits), I prefer to keep the safe 1E-9 precipitation.

There is a code if someone solves a similar problem http://www.openprocessing.org/sketch/96977 it also shows how many power estimates are required to simulate a unit unit of the path length