Previous ToC Up Next

7. Energy Checks

7.1. More Modular

Alice: Hi Bob! Did you recover from our mathematical adventures?

Bob: Well, I must say, I'm glad to be back in the normal world. Let me show you how to deal with a real example of Plummer's model, rather than with a mathematical abstraction.

Alice: You mean the code you showed me, which triggered our lengthy discussion?

Bob: Yes, but I made one modification, following your suggestion that I could factor out the spherical angle treatment, which was so similar for creating the position and velocity vectors. I've named it mkplummer2.rb. There is now one new method: 

 def spherical(r)
   vector = Vector.new
   theta = acos(frand(-1, 1))
   phi = frand(0, 2*PI)
   vector[0] = r * sin( theta ) * cos( phi )
   vector[1] = r * sin( theta ) * sin( phi )
   vector[2] = r * cos( theta )
   vector
 end

It takes only one argument, the absolute value of the vector that will be returned by this method with random values for the two spherical angles and . You can see how I invoke this method from within the inner loop of the mkplummer method: 

   nb.body.each do |b|                                                        
     b.mass = 1.0/n                                                           
     radius = 1.0 / sqrt( rand ** (-2.0/3.0) - 1.0)                           
     b.pos = spherical(radius)                                                
     x = 0.0                                                                  
     y = 0.1                                                                  
     while y > x*x*(1.0-x*x)**3.5                                             
       x = frand(0,1)                                                         
       y = frand(0,0.1)                                                       
     end                                                                      
     velocity = x * sqrt(2.0) * ( 1.0 + radius*radius)**(-0.25)               
     b.vel = spherical(velocity)                                              
   end

Alice: I like this better, yes, it is more modular.

Bob: And a lot shorter, so you can see more clearly what is happening now.

7.2. Repeatability

Alice: just to see the raw output, can you show me what a 3-body realization looks like, of Plummer's model?

Bob: Sure. Here is one: 

 |gravity> kali mkplummer2.rb -n 3
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 84602776384174400121148759980534896287
 ACS
   NBody 
     Array body
       Body body[0]
         Float mass
              3.3333333333333331e-01
         Vector pos
              1.1107754173294841e+00   8.7598177115218653e-02   1.1024225203258864e+00
         Vector vel
             -2.8803050799931779e-01   4.2827424016947752e-01  -2.4596499224386348e-01
       Body body[1]
         Float mass
              3.3333333333333331e-01
         Vector pos
             -2.5020638780233173e-01   4.5646218917292558e-01  -4.7024123561676934e-01
         Vector vel
             -5.4336097846255682e-01  -4.8179216868648034e-01   6.1323589122991451e-01
       Body body[2]
         Float mass
              3.3333333333333331e-01
         Vector pos
              7.0982750725799604e-01   4.5418343016141471e-01  -9.8938292838395003e-01
         Vector vel
              3.5888119559927067e-01   1.3607894924065120e-01  -3.6610065188231811e-01
 SCA
Alice: I see that you are echoing the seed from random number generator, as you explained before. Let us check whether we indeed can repeat the same model generation. I'd be happy to just compare the last particle. What value should we take for the seed of the random number generator?

Bob: 42.

Alice: Of course, that is the answer. But what a mundane question it was!

Bob: Different questions can have the same answer. Here is the first attempt: 

 |gravity> kali mkplummer2.rb -n 3 -s 42 | tail -2
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 42
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 42
              1.0707835396652617e-01   3.2726929775243468e-01   2.3909866955874151e-01
 SCA
And now let me repeat the command with the same seed: 

 |gravity> kali mkplummer2.rb -n 3 -s 42 | tail -2
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 42
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 42
              1.0707835396652617e-01   3.2726929775243468e-01   2.3909866955874151e-01
 SCA
And now with a different seed: 

 |gravity> kali mkplummer2.rb -n 3 -s 43 | tail -2
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 43
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 43
              1.0397916156880645e-01  -1.4513970203331489e-01   1.6854199026403588e-01
 SCA
Alice: Good! And the full values from the first 3-body version look plausible. But of course we cannot check, just by staring at them, whether they really correspond to Plummer's model. We'll have to come up with some checks.

7.3. Energy Diagnostics

Bob: I thought about this checking question, so I cobbled together a tool to do a quick energy check. Here it is: 

 #!/usr/local/bin/ruby -w
 
 require "acs"
 
 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 ekin                        # kinetic energy
     0.5*@mass*(@vel*@vel)
   end
 
   def epot(body_array)                  # potential energy
     p = 0
     body_array.each do |b|
       unless b == self
         r = b.pos - @pos
         p += -@mass*b.mass/sqrt(r*r)
       end
     end
     p
   end
 
 end
 
 class NBody
 
   attr_accessor :time, :body
 
   def initialize
     @body = []
   end
 
   def ekin                        # kinetic energy
     e = 0
     @body.each{|b| e += b.ekin}
     e
   end
 
   def epot                        # potential energy
     e = 0
     @body.each{|b| e += b.epot(@body)}
     e/2                           # pairwise potentials were counted twice
   end
 
   def write_diagnostics
     etot = ekin + epot
     print <<END
   E_kin = #{sprintf("%.3g", ekin)} ,\
  E_pot =  #{sprintf("%.3g", epot)} ,\
  E_tot = #{sprintf("%.3g", etot)}
 END
   end
 
 end
 
 include Math
 
 nb = ACS_IO.acs_read(NBody)
 nb.write_diagnostics

This is all very similar to what we did within our N-body integrators.

7.4. Measurements

Alice: Let's start with a few small realizations, to see how large the fluctuations are from case to case, in the kinetic and potential energies, as well as the total energy. 

 |gravity> kali mkplummer2.rb -n 10 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 10
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 314084175335806434408394307870262799293
   E_kin = 0.108 , E_pot =  -0.283 , E_tot = -0.175
 |gravity> kali mkplummer2.rb -n 10 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 10
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 269773554575491416225742644056394692430
   E_kin = 0.133 , E_pot =  -0.227 , E_tot = -0.0938
 |gravity> kali mkplummer2.rb -n 10 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 10
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 237787308285128236747451513485022117506
   E_kin = 0.199 , E_pot =  -0.371 , E_tot = -0.172
Bob: Better to increase the number of particles, to see whether we actually converge to something reasonable. An particle system should have deviations around three percent or so. 

 |gravity> kali mkplummer2.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 120675979305433593447309484955153184810
   E_kin = 0.147 , E_pot =  -0.279 , E_tot = -0.131
 |gravity> kali mkplummer2.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 229279117753660015758434507167504350638
   E_kin = 0.139 , E_pot =  -0.286 , E_tot = -0.147
 |gravity> kali mkplummer2.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 219014694910382831760580860722482279459
   E_kin = 0.14 , E_pot =  -0.284 , E_tot = -0.144
Bob: Good enough! Also note that the virial theorem is obeyed quiet well: the total energy has a magnitude comparable to the kinetic energy, and the potential energy has a magnitude that is twice as large.

Alice: Yes, that is encouraging. The question is: are these magnitudes right?

Bob: Well, we can check by looking up what the kinetic and potential energy should be for a Plummer model. But something tells me that you'd rather compute it for yourself.

Alice: How did you guess! But we have to compute only one quantity. Since we start with a model that is in equilibrium, we can assume the virial theorem to hold, so as soon as we know one of the three energy diagnostics that you are printing out here, we can determine the other two immediately.

Bob: In that case, the total potential energy of the system would seem to be the best candidate.

Alice: I agree. We already know the density and the potential, so it is just a matter of integration.

7.5. Two Roads

Bob: Let's see. The total potential energy must be the sum of the potential energy that each mass element feels, with respect to the rest of the star cluster:

Alice: Ah, but now you're double counting all pairwise interactions. You have to divide this by a factor two:

Bob: I guess that is right, or is it? The more I think about it, I wonder. You are talking about pairwise interactions, but we are letting each mass element feel the rest of the cluster through the whole potential . If I were to think about pairwise interactions, I would immediately think about how the cluster can be build up by bringing each new mass element in from an infinite distance. At first there would be little mass, and later, when more and more mass is accumulated, each new mass element feels more attraction, and gains more energy falling in. Yet each mass element attracts, and is in turn attracted by the earlier mass that already has fallen in. So . . . I admit, I'm a bit confused.

Alice: I'm sure it is wrong to leave out the factor two. in the above equation. But I like the notion of building a cluster up, bit by bit. Given that we have spherical symmetry, we can make life easiest by letting each radial shell fall in, as if it was prefabricated. Like building a pre-fab house: you order the parts and put it together.

Bob: So every shell with destination radius and destination thickness then has a mass . Even if the shell has a much lower density at first, and perhaps a much larger thickness, and certainly a much larger radius, when it settles in it must have that amount of mass, and since mass is conserved, this must be the mass that the pre-fab shell had before it got compressed into place. The term pre-fab is not a very good metaphor, perhaps, when we have to crunch the components, but let's not worry about words for now.

At the moment that the shell reaches its proper place, the material inside that shell has already been put there, while the material outside that shell is still waiting to fall in. So during the trip from an infinite distance with zero binding energy, down to its final destination, the infalling shell acquires an amount of potential energy of

Now if we look at it this way, and repeat the procedure for all mass shells from the inside outward, we must have for the total potential energy:

And I'm sure that there is no factor two here, no matter what you may argue about pairwise interactions.

Alice: Why do you call this energy top rather than pot?

Bob: Because I want to stress the difference with the expression you wrote down above, and you already claimed the label pot. So I just reversed the letters. Besides, my way of constructing the star cluster by letting mass rain down from infinity is surely a top-down method.

Alice: As you wish. And I agree that there should not be a factor two in your case, probably because on average, each mass element sees only half of the rest of the mass. But I must admit, I'm not totally sure that both expressions are correct. Why don't we work them both out, and see whether they boil down to the same result.

7.6. One Destination

Bob: You go first, since you wrote down the first expression. I'll hand you the expressions for the density and the potential, from way back when, eq. (12). Here they are again:

Alice: Okay. Eq. (84) then becomes:

Before working this out further, let's see whether you get at least the same expression. Your turn!

Bob: I'll try. I first need to recover the expression for the accumulative mass, eq. (26):

where I added the right-hand term to make the expression more similar to the expressions that give us density and potential. My top expression then becomes:

Bob: Hmm. That doesn't look like your expression.

Alice: It is more complicated than my expression, alright, but an integration by parts may help out. May I try?

where the first term in the next-to-last line vanishes, because the contributions at both zero and infinity are zero.

Bob: Thanks! Now it is clear that and are identically the same.

7.7. Potential Energy

Alice: All that is left to do is to solve it. Let me try another integration by parts:

where again the first term in the next-to-last line vanishes for the same reason as before.

Now how shall we tackle this integral? I wonder how we can simplify this further.

Bob: Now this is really simple enough for my taste. Yesterday I let you get away with going back to square one, but if we keep doing that, we'll never build an N-body environment.

We should be able to find this integral easily in a book with a table of integrals, or by using a symbolic package. Ah, you see, here it is, in good old Abramowitz and Stegun, right at the beginning of their list of most common integrals:

Alice: Okay, I give in.

Bob: Substituting eq. (92) into eq. (91), we get:

Alice: It looks plausible: at least it has the right dimensions and the right dependency on mass and structural parameter.

7.8. Validation

Bob: I'm really curious to see whether my code gives the potential energy that we have now calculated, at least within the statistical noise. In my code, following Aarseth et al., I have , so I should find a value of

For the kinetic and total energy, we should similarly expect:

Alice: And look, those are indeed the type of numbers we got. Why don't you create a few more 1000-particle realizations.

Bob: Here are a few more. 

 |gravity> kali mkplummer2.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 302041592721379705309481006390526637954
   E_kin = 0.142 , E_pot =  -0.275 , E_tot = -0.132
 |gravity> kali mkplummer2.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 8717335803285052460172785737126512016
   E_kin = 0.145 , E_pot =  -0.287 , E_tot = -0.142
 |gravity> kali mkplummer2.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
              actual seed used  : 283452213331896683662695500815454999402
   E_kin = 0.144 , E_pot =  -0.298 , E_tot = -0.153
You're right: the numbers work out beautifully.

Alice: Congratulations!
Previous ToC Up Next