Previous ToC Up Next

1. Leapfrog

1.1. The Basic Equations

Alice: Hi, Bob! Any luck in getting a second order integrator to work?

Bob: No problem, it was easy! Actually, I got two different ones to work, and a fourth order integrator as well.

Alice: Wow, that was more than I expected!

Bob: Let's start with the second order leapfrog integrator.

Alice: Wait, I know what a leapfrog is, but we'd better make some notes about how to present the idea to your students. How can we do that quickly?

Bob: Let me give it a try. The name leapfrog comes from one of the ways to write this algorithm, where positions and velocities `leap over' each other. Positions are defined at times , spaced at constant intervals , while the velocities are defined at times halfway in between, indicated by , where . The leapfrog integration scheme then reads:

Note that the accelerations are defined only on integer times, just like the positions, while the velocities are defined only on half-integer times. This makes sense, given that : the acceleration on one particle depends only on its position with respect to all other particles, and not on its or their velocities. Only at the beginning of the integration do we have to set up the velocity at its first half-integer time step. Starting with initial conditions and , we take the first term in the Taylor series expansion to compute the first leap value for :

We are then ready to apply the first equation above to compute the new position , using the first leap value for . Next we compute the acceleration , which enables us to compute the second leap value, , using the second equation above, and then we just march on.

A second way to write the leapfrog looks quite different at first sight. Defining all quantities only at integer times, we can write:

This is still the same leapfrog scheme, although represented in a different way. Notice that the increment in is given by the time step multiplied by , effectively equal to . Similarly, the increment in is given by the time step multiplied by , effectively equal to the intermediate value . In conclusion, although both positions and velocities are defined at integer times, their increments are governed by quantities approximately defined at half-integer values of time.

1.2. Time reversibility

Alice: A great summary! I can see that you have taught numerical integration classes before. At this point an attentive student might be surprised by the difference between the two descriptions, which you claim to describe the same algorithm. They may doubt that they really address the same scheme.

Bob: How would you convince them?

Alice: An interesting way to see the equivalence of the two descriptions is to note the fact that the first two equations are explicitly time-reversible, while it is not at all obvious whether the last two equations are time-reversible. For the two systems to be equivalent, they'd better share this property. Let us inspect.

Starting with the first set of equations, even though it may be obvious, let us write out the time reversibility. We will take one step forward, taking a time step , to evolve to , and then we will take one step backward, using the same scheme, taking a time step . Clearly, the time will return to the same value since , but we have to inspect where the final positions and velocities are indeed equal to their initial values . Here is the calculation, resulting from applying the first set of equations twice:

In an almost trivial way, we can see clearly that time reversal causes both positions and velocities to return to their old values, not only in an approximate way, but exactly. In a computer application, this means that we can evolve forward a thousand time steps and then evolve backward for the same length of time. Although we will make integration errors (remember, leapfrog is only second-order, and thus not very precise), those errors will exactly cancel each other, apart from possible round-off effects, due to limited machine accuracy.

Now the real fun comes in, when we inspect the equal-time version, the second set of equations you presented:

In this case, too, we have exact time reversibility. Even though not immediately obvious from an inspection of your second set of equations, as soon as we write out the effects of stepping forward and backward, the cancellations become manifest.

Bob: A good exercise to give them. I'll add that to my notes.

1.3. A New Driver

Alice: Now show me your leapfrog, I'm curious.

Bob: I wrote two new drivers, each with its own extended Body class. Let me show you the simplest one first. Here is the leapfrog driver integrator_driver1.rb: 

 require "lbody.rb"
 
 include Math
 
 dt = 0.001            # time step
 dt_dia = 10           # diagnostics printing interval
 dt_out = 10           # output interval
 dt_end = 10           # duration of the integration
 ##method = "forward"    # integration method                               
 method = "leapfrog"   # integration method                                
 
 STDERR.print "dt = ", dt, "\n",
       "dt_dia = ", dt_dia, "\n",
       "dt_out = ", dt_out, "\n",
       "dt_end = ", dt_end, "\n",
       "method = ", method, "\n"
 
 b = Body.new
 b.simple_read
 b.evolve(method, dt, dt_dia, dt_out, dt_end)

Same as before, except that now you can choose your integrator. The method evolve, at the end, now has an extra parameter, namely the integrator.

Alice: And you can supply that parameter as a string, either "forward" for our old forward Euler method, or "leapfrog" for your new leapfrog integrator. That is very nice, that you can treat that choice on the same level as the other choices you have to make when integrating, such as time step size, and so on.

Bob: And it makes it really easy to change integration method: one moment you calculate with one method, the next moment with another. You don't even have to type in the name of the method: I have written it so that you can switch from leapfrog back to forward Euler with two key strokes: you uncomment the line 

 ##method = "forward"    # integration method

and comment out the line 

 method = "leapfrog"   # integration method

Alice: It is easy to switch lines in the driver, but I'm really curious to see how you let Ruby actually make that switch in executing the code differently, replacing one integrator by another!

1.4. An Extended Body Class

Bob: Here is the new version of our Body class, in a file called lbody.rb with l for leapfrog. It is not much longer than the previous file body.rb, so let me show it again in full: 

 require "vector.rb"
 
 class Body
 
   attr_accessor :mass, :pos, :vel
 
   def initialize(mass = 0, pos = Vector[0,0,0], vel = Vector[0,0,0])
     @mass, @pos, @vel = mass, pos, vel
   end
 
   def evolve(integration_method, dt, dt_dia, dt_out, dt_end)
     time = 0
     nsteps = 0
     e_init
     write_diagnostics(nsteps, time)
 
     t_dia = dt_dia - 0.5*dt
     t_out = dt_out - 0.5*dt
     t_end = dt_end - 0.5*dt
 
     while time < t_end
       send(integration_method,dt)
       time += dt
       nsteps += 1
       if time >= t_dia
         write_diagnostics(nsteps, time)
         t_dia += dt_dia
       end
       if time >= t_out
         simple_print
         t_out += dt_out
       end
     end
   end
 
   def acc
     r2 = @pos*@pos
     r3 = r2*sqrt(r2)
     @pos*(-@mass/r3)
   end    
 
   def forward(dt)
     old_acc = acc
     @pos += @vel*dt
     @vel += old_acc*dt
   end
 
   def leapfrog(dt)
     @vel += acc*0.5*dt
     @pos += @vel*dt
     @vel += acc*0.5*dt
   end
 
   def ekin                        # kinetic energy
     0.5*(@vel*@vel)               # per unit of reduced mass
   end
 
   def epot                        # potential energy
     -@mass/sqrt(@pos*@pos)        # per unit of reduced mass
   end
 
   def e_init                      # initial total energy
     @e0 = ekin + epot             # per unit of reduced mass
   end
 
   def write_diagnostics(nsteps, time)
     etot = ekin + epot
     STDERR.print <<END
 at time t = #{sprintf("%g", time)}, after #{nsteps} steps :
   E_kin = #{sprintf("%.3g", ekin)} ,\
  E_pot =  #{sprintf("%.3g", epot)} ,\
  E_tot = #{sprintf("%.3g", etot)}
              E_tot - E_init = #{sprintf("%.3g", etot-@e0)}
   (E_tot - E_init) / E_init =#{sprintf("%.3g", (etot - @e0) / @e0 )}
 END
   end
 
   def to_s
     "  mass = " + @mass.to_s + "\n" +
     "   pos = " + @pos.join(", ") + "\n" +
     "   vel = " + @vel.join(", ") + "\n"
   end
 
   def pp               # pretty print
     print to_s
   end
 
   def simple_print
     printf("%24.16e\n", @mass)
     @pos.each{|x| printf("%24.16e", x)}; print "\n"
     @vel.each{|x| printf("%24.16e", x)}; print "\n"
   end
 
   def simple_read
     @mass = gets.to_f
     @pos = gets.split.map{|x| x.to_f}.to_v
     @vel = gets.split.map{|x| x.to_f}.to_v
   end
 
 end

1.5. Minimal Changes

Alice: Before you explain to me the details, remember that I challenged you to write a new code while changing or adding at most a dozen lines? How did you fare?

Bob: I forgot all about that. It seemed so unrealistic at the time. But let us check. I first corrected this mistake about the mass factor that I had left out in the file body.rb. Then I wrote this new file lbody.rb. Let's do a diff: 

    |gravity> diff body.rb lbody.rb
    11c11
    <   def evolve(dt, dt_dia, dt_out, dt_end)
    - -
    >   def evolve(integration_method, dt, dt_dia, dt_out, dt_end)
    22c22
    <       evolve_step(dt)
    - -
    >       send(integration_method,dt)
    36c36
    <   def evolve_step(dt)
    - -
    >   def acc
    39c39,49
    <     acc = @pos*(-@mass/r3)
    - -
    >     @pos*(-@mass/r3)
    >   end    
    > 
    >   def forward(dt)
    >     old_acc = acc
    >     @pos += @vel*dt
    >     @vel += old_acc*dt
    >   end
    > 
    >   def leapfrog(dt)
    >     @vel += acc*0.5*dt
    41c51
    <     @vel += acc*dt
    - -
    >     @vel += acc*0.5*dt
To wit: four lines from the old code have been left out, and twelve new lines appeared.

Alice: Only twelve! You did it, Bob, exactly one dozen, indeed.

Bob: I had not realized that the changes were so minimal. While I was playing around, I first added a whole lot more, but when I started to clean up the code, after it worked, I realized that most of my changes could be expressed in a more compact way.

Alice: Clearly Ruby is very compact.

Bob: Let me step through the code.
Previous ToC Up Next