4
\$\begingroup\$

I am looking at various integration methods for my n-body simulation and I'm slightly confused about actual implementation of leapfrog integration.

According to the wikipedia page leapfrog method is defined like this:

leapfrog forumla from wikipedia

When I imagine how to code it, i come up with something like this:

while true:
    dt = time since last frame

    for every object:
        object.position += object.velocity * dt

    for every object:
        acceleration = object.calculate_acceleration()
        object.velocity += acceleration * dt

    render()

But this is Euler integration, right?

So moving the initial velocity back by one half step makes a first order method into a second order method? But for me a small change of initial velocity is not important, it will be set as a random value anyway...

Why I won't just use velocity verlet: This is for a game where tens of thousands particles fly around affected by forces roughly resembling gravity. I'd like to keep amount of data stored in every partile as low as possible to improve cache locality, which is why I don't want to use a method that requires storing acceleration (and recalculating it would be too slow).

\$\endgroup\$
1
  • 1
    \$\begingroup\$ Indeed, what you wrote is called semi-implicit or symplectic Euler. Both Verlet/Leap-frog and the Semi-implicit Euler methods are symplectic (which is a special class of integrators that preserve the energy of the system if conservatory forces are used - e.g. without friction). Here's a quick and dirty (not mathematically ladden) reference codeflow.org/entries/2010/aug/28/… \$\endgroup\$
    – teodron
    Dec 16, 2015 at 8:29

1 Answer 1

7
\$\begingroup\$

The leapfrog method works by using a variable at a non-integer time step to determine the new value of a different variable.

Image demonstrating the leapfrog method
(source: drexel.edu)

So the new position values are calculated using the velocity half a time step ahead of the position. This is done for the same reason that a similar method is used in the mid-point method (however the leapfrog method has some advantages i.e. it is reversible and hence useful for oscillating systems).

If you want to run the leapfrog method one possible way is using,

while true
    dt = time since last frame

    for every object:
        // If first time object is moved determine velocity at midpoint
        if object.firstmove
            object.velocity += acceleration * dt / 2 
        object.position += object.velocity * dt


    for every object:
        acceleration = object.calculate_acceleration()
        object.velocity += acceleration * dt

    render()

So the objects store the velocity not at the given time step but at time step + 1/2.

Potential problems

Since you are not storing the velocity at a particular time step but instead at time step + 1/2 any processes which use this velocity for calculations may not be correct.

This is also true for the acceleration, the leapfrog method assumes that the acceleration does not depend on the velocity. So the acceleration only depends on the position, hence the velocity 'leaps' over the position.

However, if your object.calculate_acceleration() function uses the velocity then it will not be a true leapfrog method as the velocity will be leaping over itself and you will be using information both from the current time step and the time step + 1/2.

Example; friction in the form of drag is often dependent on velocity^2 (as one would expect the faster the object the larger the drag).

Other notes

The added line in my code snippet is simply changing the initial velocity to the new time step. If this is not important then you simply return to the normal Euler integration

On the order of the method

According to wikipedia the leapfrog method can also be written as,

enter image description here

which is a second order method. These equations can also be used to perform your integration and may be easier to understand. Again note that here the acceleration cannot depend on the velocity otherwise the second equation becomes much more complicated.

\$\endgroup\$
2
  • \$\begingroup\$ This is a great answer, one thing I would propose adding is a small example that calculate_acceleration with friction generally requires the velocity. \$\endgroup\$
    – Elva
    Dec 16, 2015 at 16:59
  • \$\begingroup\$ Ah yes I'll mention that \$\endgroup\$
    – Malrig
    Dec 16, 2015 at 17:00

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .