Softening

Close Encounters

Alice: I’m very glad to see that we can integrate eight bodies in a cold collapse system. This is quite a bit more demanding than integrating a handful of bodies in a virialized system. However, in both cases, sooner or later there will be close encounters between two or more of the particles. Our code will never be able to handle all of those close encounters. No matter how small a time step we give it, sooner or later there will be particles that approach each other closely enough to have a near miss that takes less time than the time step size. This will necessarily lead to large numerical errors.

Bob: This is of course why people have introduced variable time steps, as well as a whole order set of algorithmic tools to tame the unruly behavior of particles that get too close to the inverse square singularities of Newtonian gravity.

Alice: Soon we will introduce those extensions in our codes, but for now, there are more urgent things on our agenda. I guess we just have to live with it, and make sure the students realize that this first N-body tool is not to be trusted under all circumstances.

Bob: Hmm. I don’t much like the idea of giving someone a tool that cannot be trusted. How about adding softening, as an option?

Alice: You mean to soften the potential, from an inverse square law to a form that remains finite in the center?

Bob: Indeed. We start from the singular Newtonian potential energy between two particles with positions and and masses :

The standard softening approach is to replace this by a regular variant, simply by adding the square of a small quantity :

When you differentiate this modified potential with respect to the position of a particle, you obtain a modified acceleration:

And of course, in the limit that , this last equation again returns to the Newtonian gravitational acceleration.

Fuzzy-Point Particles

Alice: Yes, this is what is often used in collisionless stellar dynamics, to suppress the effect of close encounters. I can’t say I’m very happy with this softening approach, since it’s not the real thing. It is purely a mathematical trick, to avoid numerical problems.

Bob: Well, you can give it a physical interpretation. Instead of using point particles, which are not very physical in the first place, each particle gets a more extended mass distribution. In fact, you can easily show that a softened potential corresponds to a mass distribution given by a polytrope of index five, better known as a Plummer mass distribution:

Alice: But look, your mass distribution stretches all the way to infinity! Even though most of the mass in concentrated in a small region, with a radius of order the softening length . You solution works, in the sense of avoiding singularities, and it gives a roughly reasonable answer, but it does come at the cost of smearing each particle all over space.

Bob: It would be quite easy to use a different mass distribution, corresponding a finite support. This is what people so who work with SPH particles, for example. However, for or current purpose, the main thing is to provide a tool that works, and we can worry later about aesthetic details.

Alice: Okay. Even though I can’t say I’m very happy with it, I see your point, and it is certainly safer to give the students a tool that is guaranteed to give finite answer.

Bob: It should be easy to add softening to our code. Time to create another version for our N-body code! So we will call this new file rknbody9.rb. Well, this will take me a while.

Alice: Okay, I’m way behind in reading the astro-ph abstracts. This will give me a chance to catch up. I’ll come back when I’ve gone through them.

A New Driver

Bob: Here it is, the new version of our N-body code, now with softening build in. It was quite straightforward to make the changes. First of all, here is the new driver:

 require "rknbody9.rb"

 include Math

 eps = 0
 dt = 0.001           # time step
 dt_dia = 2.1088      # diagnostics printing interval
 dt_out = 2.1088      # output interval
 dt_end = 2.1088      # duration of the integration
 init_out = false     # initial output requested ?
 x_flag = false       # extra diagnostics requested ?
 #method = "forward"  # integration method
 #method = "leapfrog" # integration method
 #method = "rk2"      # integration method
 method = "rk4"       # integration method

 STDERR.print "eps = ", eps, "\n",
       "dt = ", dt, "\n",
       "dt_dia = ", dt_dia, "\n",
       "dt_out = ", dt_out, "\n",
       "dt_end = ", dt_end, "\n",
       "init_out = ", init_out, "\n",
       "x_flag = ", x_flag, "\n",
       "method = ", method, "\n"

 nb = Nbody.new
 nb.simple_read
 nb.evolve(method, eps, dt, dt_dia, dt_out, dt_end, init_out, x_flag)

As you can see, minimal differences, contained in three lines. The method evolve has an extra parameter, eps, the softening length. The default value is zero, which means no softening at all. The third new line is where the value of eps is echoed on the standard error stream.

Alice: So now evolve has eight parameters. At some point we may want to think about grouping them together, perhaps creating a class for them, since there is clear substructure: two flags controlling the amount of output, three variables giving intervals between output times, and three other variables.

Bob: But not now.

Alice: Not now, no. Can you show me the code itself?

A Code with Softening

Bob: Here it is. Almost all changes speak for themselves.

 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 calc(softening_parameter, body_array, time_step, s)
     ba  = body_array
     dt = time_step
     eps = softening_parameter
     eval(s)
   end

   def acc(body_array, eps)
     a = @pos*0                              # null vector of the correct length
     body_array.each do |b|
       unless b == self
         r = b.pos - @pos
         r2 = r*r + eps*eps
         r3 = r2*sqrt(r2)
         a += r*(b.mass/r3)
       end
     end
     a
   end

   def ekin                         # kinetic energy
     0.5*@mass*(@vel*@vel)
   end

   def epot(body_array, eps)        # potential energy
     p = 0
     body_array.each do |b|
       unless b == self
         r = b.pos - @pos
         p += -@mass*b.mass/sqrt(r*r + eps*eps)
       end
     end
     p
   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 ppx(body_array, eps)         # pretty print, with extra information (acc)
     STDERR.print to_s + "   acc = " + acc(body_array, eps).join(", ") + "\n"
   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

 class Nbody

   attr_accessor :time, :body

   def initialize
     @body = []
   end

   def evolve(integration_method, eps, dt, dt_dia, dt_out, dt_end,
              init_out, x_flag)
     @dt = dt
     @eps = eps
     nsteps = 0
     e_init
     write_diagnostics(nsteps, x_flag)
     t_dia = dt_dia - 0.5*dt
     t_out = dt_out - 0.5*dt
     t_end = dt_end - 0.5*dt

     simple_print if init_out

     while @time < t_end
       send(integration_method)
       @time += dt
       nsteps += 1
       if @time >= t_dia
         write_diagnostics(nsteps, x_flag)
         t_dia += dt_dia
       end
       if @time >= t_out
         simple_print
         t_out += dt_out
       end
     end
   end

   def calc(s)
     @body.each{|b| b.calc(@eps, @body, @dt, s)}
   end

   def forward
     calc(" @old_acc = acc(ba,eps) ")
     calc(" @pos += @vel*dt ")
     calc(" @vel += @old_acc*dt ")
   end

   def leapfrog
     calc(" @vel += acc(ba,eps)*0.5*dt ")
     calc(" @pos += @vel*dt ")
     calc(" @vel += acc(ba,eps)*0.5*dt ")
   end

   def rk2
     calc(" @old_pos = @pos ")
     calc(" @half_vel = @vel + acc(ba,eps)*0.5*dt ")
     calc(" @pos += @vel*0.5*dt ")
     calc(" @vel += acc(ba,eps)*dt ")
     calc(" @pos = @old_pos + @half_vel*dt ")
   end

   def rk4
     calc(" @old_pos = @pos ")
     calc(" @a0 = acc(ba,eps) ")
     calc(" @pos = @old_pos + @vel*0.5*dt + @a0*0.125*dt*dt ")
     calc(" @a1 = acc(ba,eps) ")
     calc(" @pos = @old_pos + @vel*dt + @a1*0.5*dt*dt ")
     calc(" @a2 = acc(ba,eps) ")
     calc(" @pos = @old_pos + @vel*dt + (@a0+@a1*2)*(1/6.0)*dt*dt ")
     calc(" @vel += (@a0+@a1*4+@a2)*(1/6.0)*dt ")
   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, @eps)}
     e/2                           # pairwise potentials were counted twice
   end

   def e_init                      # initial total energy
     @e0 = ekin + epot
   end

   def write_diagnostics(nsteps, x_flag)
     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
     if x_flag
       STDERR.print "  for debugging purposes, here is the internal data ",
                    "representation:\n"
       ppx
     end
   end

   def pp                           # pretty print
     print "     N = ", @body.size, "\n"
     print "  time = ", @time, "\n"
     @body.each{|b| b.pp}
   end

   def ppx                          # pretty print, with extra information (acc)
     print "     N = ", @body.size, "\n"
     print "  time = ", @time, "\n"
     @body.each{|b| b.ppx(@body, @eps)}
   end

   def simple_print
     print @body.size, "\n"
     printf("%24.16e\n", @time)
     @body.each{|b| b.simple_print}
   end

   def simple_read
     n = gets.to_i
     @time = gets.to_f
     for i in 0...n
       @body[i] = Body.new
       @body[i].simple_read
     end
   end

 end

Details

Alice: Even though the changes may speak for themselves, I have some questions. First of all, the value of eps has to be passed on from the driver, where it is defined, through evolve into Nbody and then down to the methods within Body that do all the hard work.

Bob: First of all, I gave the Nbody class an extra instance variable, @eps, which stores the value of the softening. As soon as evolve is executed, within the Nbody class, the first thing it does is assign the proper value to @eps, as well as to @dt, as was done already in our previous version:

     @dt = dt
     @eps = eps

Alice: I see. In that way, you don’t have to give an extra argument to the integration methods, for example: they can just pick up the value of the softening length from @eps, to which they have automatic access, as Nbody class methods. But of course they do have to pass that value down to the particles, which are realized as instances of the Body class, since otherwise the particles would not know what softening to use.

Bob: The one thing I didn’t like very much is that the lines in the integration methods have become somewhat longer. Forward Euler, for example, has grown now from:

   def forward
     calc(" @old_acc = acc(ba) ")
     calc(" @pos += @vel*dt ")
     calc(" @vel += @old_acc*dt ")
   end

to:

   def forward
     calc(" @old_acc = acc(ba,eps) ")
     calc(" @pos += @vel*dt ")
     calc(" @vel += @old_acc*dt ")
   end

I’m not too happy with the fact that acc now has to get a second argument. But that’s the way it is.

If I really would want to make the lines shorter, I could use shorter variables than ba and eps, for example a and e. But let us not spend more time on such niceties, which give us diminishing return in clarity at the cost of making the code more complex and hence less clear.

Alice: I fully agree. Instead, let’s see how your new code performance in our cold collapse experiment.