Previous ToC Up Next

9. Energy Conservation

9.1. Estimating Numerical Errors

Bob: That was quite a long session yesterday, when we were chasing those bugs while we were trying to implement a better Body class!

Alice: Yes, but we learned a lot about Ruby.

Bob: And even a thing or two about debugging.

Alice: Perhaps this is a good time to return to our real project: to build progressively better integrators. But before we move on to second-order integrators, as well as introducing graphics and other complications, I would like to do one more test. Even though we now know that our forward Euler scheme converges, we don't know yet how well it conserves energy. Shouldn't we test that explicitly?

Bob: Great minds think alike. I had come to the same conclusion. For the two-body problem, of course you can check all kind of things. You can measure how accurately a bound orbit remains the same ellipse, as it should be, rather than slowly drifting and changing shape, and whether the motion traverses the ellipse at the right speed, and so on. The reason is that we have an analytic solution for the two-body problem.

For the many-body problem in general, however, we are a lot less lucky. All we know is that the total energy and total angular momentum of the system are conserved. That gives us one scalar and one vector, or in three dimensions it gives us four conserved quantities. In addition, there are the position and velocity of the center-of-mass motion, which give us another six conserved quantities, ten in total. So we cannot compare all the details of the motion of the particles with theory, not by a long shot. An N-body system in three dimensions has degrees of freedom, far more than the 10 handles we have in our 2-body case.

However, things are not as bad as they may seem. If you make errors in your integration, and you always make errors at some level, those errors will let the simulated orbit drift away from the true orbit in a random fashion. It would require a very clever conspiracy for those errors to work together in such a way as to keep the energy conserved to a high degree of accuracy, while still introducing much larger errors elsewhere. In practice, therefore, checking energy conservation in an N-body system has become the standard test to see how accurate the integration of an N-body system is.

Alice: Another way of saying this is to picture the phase space of the whole system.

Bob: You mean the six-dimensional space in which you plot for each particle its position and velocity?

Alice: I was thinking about the alternative way of picturing the whole system as one point in a -dimensional system. The evolution of an N-body system can then be viewed as the complex motion of this one point through this huge space. Now in this space you can define layers of subspaces, on each of which the total energy of the system is constant. Once the particle starts on such a hypersurface, it should stay there, because energy is conserved. But if we make arbitrary errors in computing its motion, our simulated particle will show some numerical drift in all dimensions available, including the direction perpendicular to the energy hypersurface.

Bob: Which reaches the same conclusion as what I just summarized.

Alice: Yes. Since the dynamics of this one point in dimensions is incredibly complicated, the distance away from the energy hypersurface will typically be of the same order of magnitude as the total length of the drift of that master particle that symbolizes our whole N-body system.

Bob: More mathematically minded students might prefer your explanation; I think I'll just stick to my more lowbrow explanation. But however we tell the story, we should of course warn them that this argument fails when we are dealing with an integrator that has a built-in type of energy conservation.

Alice: Yes, indeed, that is very important. It would be easy to construct an integrator that projects the orbit of the master particle in dimensions back onto the original constant-energy hypersurface. But that would be cheating. It would get rid of only one degree of freedom of the total error, and leave the other degrees of freedom uncorrected. It would look great when you test for energy conservation, but at the cost of having lost your handle on checking what is going on.

Bob: Or the way I would put it, you could just wiggle one particle far out, perhaps a particle that has already escaped from the system. By changing the velocity of that particle a little, you will change its kinetic energy, and so you can accurately balance the numerical errors that occur in the rest of the system. But of course that does not make the rest of the system behave in a more accurate way.

Alice: Of course, the two examples that you and I just cooked up are extreme. It is not at all a bad idea to conserve energy, if your approach is part of a systematic way to make the whole integration more accurate. Symplectic integrators are an example.

Bob: I have heard of the name, but I must admit I don't know what they are.

Alice: They are a lot of fun, from a mathematical point of view. How useful they may be for real N-body simulations, that is not too clear. They are more accurate, for sure, but they are hard to generalize to individual time steps. And for our game, we have no other choice but to go to individual time steps. Perhaps we can come back to this topic, some day, but not any time soon. But just to give you one example of the simplest symplectic integrator: the good old leapfrog scheme, also known in some areas of physics as the Verlet algorithm.

Bob: Which is probably the next algorithm we want to implement.

Alice: Yes. So shall we first test energy conservation of our forward Euler code?

9.2. A Driver

Bob: Well, I have a surprise for you. Last night I continued to tinker with Ruby, and I had so much fun doing so that I stayed at it long enough to wind up with a version of our code that produces energy conservation diagnostics.

Alice: How nice! Can I see it?

Bob: First of all, here is the driver, in file euler.rb, the only part of the program that we will change, while exploring the effect of different time step sizes and integration times: 

 require "body.rb"
 
 include Math
 
 dt = 0.01           # time step
 dt_dia = 1          # diagnostics printing interval
 dt_out = 1          # output interval
 dt_end = 1          # duration of the integration
 
 STDERR.print "dt = ", dt, "\n",
       "dt_dia = ", dt_dia, "\n",
       "dt_out = ", dt_out, "\n",
       "dt_end = ", dt_end, "\n"
 
 b = Body.new
 b.simple_read
 b.evolve(dt, dt_dia, dt_out, dt_end)

Alice: So at the bottom you create a new particle, which will represent the relative motion between the two bodies in our system, with the call to Body.new, then you read in the initial conditions, as we did before, and then you tell the particle to evolve itself.

Bob: Yes, with evolve I mean that we invoke here a method that computes the evolution of the system in time. I have made this a method in the Body class, so we can invoke it with b.evolve. The only extra thing needed is to tell the particle how long to integrate, with what time step, and how frequently to report the degree of energy conservation and the actual position and velocity of the relative motion.

Alice: So the actual program is a mere three lines, at the bottom, with the four lines of assignments at the top just the way to specify each of those four parameters you just mentioned.

Bob: Exactly. Ruby is very economical! And just as a courtesy, I decided to let the program print out the values of the four parameters. We can then change those parameters in the current file to our hearts' content, and we can always check from earlier outputs what parameters were used for a given integration.

Alice: Good idea. Clean and careful. Of course, some day soon we will want to give those values as command line options, rather than always editing the file.

Bob: I suppose so, although we'll have to think about that. It is of course possible to give Unix-style command line options, but there are several other ways as well; in fact Ruby has a nice built-in way to do that.

Alice: Later. Can I have a look at what the Body class now looks like? It must have grown quite a bit, since it now contains not only the forward Euler integrator, but also your whole evolve method, as well as whatever you wrote to get all the energy diagnostics done.

Bob: It's not that bad, actually. Ruby really is a compact language. Here, let me show you the whole file body.rb first, and then we can go through each of the methods used.

9.3. The Body Class

 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(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
       evolve_step(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 evolve_step(dt)
     r2 = @pos*@pos
     r3 = r2*sqrt(r2)
     acc = @pos*(-@mass/r3)
     @pos += @vel*dt
     @vel += acc*dt
   end
 
   def ekin                     # kinetic energy
     0.5*(@vel*@vel)
   end
 
   def epot                     # potential energy
     -1/sqrt(@pos*@pos)
   end
 
   def e_init                   # initial total energy
     @e0 = ekin + epot
   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

Alice: Indeed, not as long as I would have guessed. Can you give me a guided tour?

Bob: My pleasure!
Previous ToC Up Next