Computational Physics, 1B, Term 2

Worked solution for BONKERS example

pretty pictures | one-dimensional solution | two-dimensional solution
You can download all the files in a single zip file.


Some illustrations

TEN PARTICLES - equilibrium thermodynamics

ex10a In the first multiplot, the upper half shows the positions of ten particles, and the low half shows the complete state - both the positions (horizontal coordinates) and the velocities (vertical coordinates).

The blue particles on the left have masses roughly equal to 1.0. The green on the right have masses roughly equal to 4.0. One of the predictions of statistical physics is that for a large system at equilibrium, the probability of any state is just exp(-beta E)/Z, where E is the energy of the state. This prediction has astonishing consequences. First, each particle's velocity should have a Gaussian distribution. Second, all valid positions are equiprobable, and the probability of the positions and the velocities is a separable distribution -- positions and velocities are independent. All valid positions are equiprobable, whatever the masses of the particles are. So changing the particle masses may change their typical velocities but it has no effect at all on their typical positions.

ex10b This figure shows scatter plots of the velocties of particles 1 and 2 (left) and particles 9 and 10 (right).

A striking prediction of statistical physics is that both distributions should be Gaussian, and the standard deviations of v1 and v10 (velocities of particles 1 and 10) should be in the ratio 2:1, because the mass ratio is 1:4. ex10d This prediction is well borne out in the data, though the data show an intriguing dimple in P(v1) -- is this dimple because N=10 is too small for the large-N equilibrium theory to apply? Or is it because the particular random masses chosen here give a system that has a slow mixing time, so the points generated were not actually a good representation of the equilibrium distribution?

ex10c This figure shows the joint distribution of the positions of particles 1 and 2, 9 and 10.

ELEVEN PARTICLES - non-equilibrium thermodynamics

The next pin-up in my picture gallery illustrates very different ideas. In the 10-particle example above, ex11a we were assuming that equilibrium had been reached and testing the predictions of equilibrium theory. In my eleven-particle example, I have a large purple particle which I call the piston plopped in the middle. The dynamics are the same as normal, but its mass (40) is so much bigger than the others, the behaviour of them all is radically altered. In this simulation I set off all the particles with equal velocities. This means the situation is not at all typical of an equilibrium thermal distribution -- the big guy has far too much energy -- he has roughly 75% of the energy and at equilibrium he'd expect to have 9%. So what happens for the whole simulation is non-equilibrium. One way of thinking of the piston is that he's a lot like a wall, pushing down on the green community and rushing away from the blue - or vice versa. What happens to a gas that is rapidly compressed or expanded by a moving piston? Well, it depends exactly which sort of expansion we are discussing, but an adiabatic expansion is one sort. If the blue community and the green community both behave like adiabatically expanded or compressed gases, what would you predict their temperature and pressure to do as a function of 'volume' (which means the accessible volume given the piston position)? Maybe the motion of the piston can be modelled using the idea that its acceleration is the result of the pressure difference between the green gas on the right and the blue gas on the left? Under this theory, the piston should move approximately like a simple harmonic oscillator. What frequency does the adiabatic theory predict?

ex11b Here's the positions (horizontal axis) of all the particles as a function of time (vertical axis).

ex11c Here's the same thing but showing more data.

ex11d Here's the total kinetic energy of the blue guys, the green guys, and the purple guy, for the first 40 time units;

ex11e and for the first 200 time units. The way in which the excess energy of the piston is gradually dissipated and the energy sloshes between the particles deserves hours of study and thought. If the piston had very large mass, how would it ever lose energy to the gases? Perhaps this is a good test-bed for thinking about concepts of fluctuations and dissipations.


One-dimensional solution

Bonk1.cc is a worked solution for the first task:
...a single function - let's call it collide - which receives two particles with masses m1 and m2 and velocities u1 and u2, and returns the new velocities v1 and v2 of the two particles after an elastic collision between the two.
It has a structure called a particle. Each particle has a one-dimensional coordinate x, a velocity v, and an inverse-mass im. I use an inverse-mass so that infinitely-massive walls can be handled directly. Such a wall has inverse-mass zero.
// Bonk1.cc
//   program to test a two-particle 'collide' function in one dimension

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <string>

using namespace std;

#define ranf() \
  ((double)random()/(1.0+(double)RAND_MAX)) // Uniform from interval [0,1) */

struct particle {
  double x    ; // x coordinate
  double im   ; // inverse mass
  double v    ; // velocity
  double p    ; // momentum
  double T    ; // energy
} ; // Note the definition of a structure ends with a semicolon

void v2p( particle &a )
  // propagate a velocity into a momentum
  // p = m * v
{
  if( a.im > 0.0 ) 
    a.p =  a.v / a.im ;
  else // infinite mass particle can't be handled correctly
    a.p = 0.0 ;

  // update the energy too
  a.T = 0.5*a.v * a.p ; 
}

void collide ( particle & a , particle & b  )
  // two particles collide elastically
{
  // find the relative velocity 
  double velocity_along_line =  a.v - b.v;
  // find a's inverse-mass fraction
  double af = a.im / ( a.im + b.im ) ; 
  double bf = b.im / ( a.im + b.im ) ; 
  // reverse the c.o.m. velocity of each along the line of collision
  a.v -= 2.0 * af * velocity_along_line ;
  b.v += 2.0 * bf * velocity_along_line ;
  // The collision has now happened.
  // Finally update the other state information
  v2p(a); 
  v2p(b);
}

// Test the collide function using random initial conditions
int main(int argc, char* argv[])
{
  particle   a , b ;
  int Tests = 10 ; 
  // set defaults
  // read in any command-line arguments
  if(argc>1)   {
    sscanf( argv[1], "%d", &Tests ) ; // put the first command-line argument 
  }

  for ( int t = 1 ; t <= Tests ; t ++ ) {
    a.im    = ranf() + 0.1 ; // random inverse-mass
    a.v     = ranf()-0.5 ;   // random velocity
    b.im    = ranf() + 0.1 ; // random inverse-mass
    b.v     = ranf()-0.5 ;   // random velocity
    v2p(a);
    v2p(b);
    double MomentumBefore = a.p + b.p ;
    double EnergyBefore = a.T + b.T ;
    cout << "a.m:" << 1.0/a.im << "," << "b.m:" << 1.0/b.im << endl ;
    cout << a.v << "," << b.v << " -> " ;
    collide(a,b) ; 
    cout << a.v << "," << b.v << endl;
    double MomentumAfter = a.p + b.p ;
    double EnergyAfter = a.T + b.T ;
    cout << "(E,P) before: "
	 << EnergyBefore << ", "
	 << MomentumBefore << ";" << endl ; 
    cout << "(E,P) after:  "
	 << EnergyAfter << ", "
	 << MomentumAfter << "." << endl ;
    cout << "Percentage change in energy: "
	 << 100.0*(EnergyAfter-EnergyBefore)/
      (0.5*(EnergyAfter+EnergyBefore)) << endl; 
    cout << "Percentage change in momentum: "
	 << 100.0*(MomentumAfter-MomentumBefore)/
      (0.5*(MomentumAfter+MomentumBefore)) << endl;
    cout << endl; 
  }
  return 0;
}

This test program should be compiled and run in a shell. It reports for every test how big the errors in momentum and energy are.


Bonkers1.cc is a worked solution for the next task:
Write code that uses a one-dimensional collide function to simulate the motion of N particles in a one-dimensional box. Each particle can collide only with its immediate neighbours; each of the two end particles can collide with one wall too. To simulate the dynamics, you must identify the times at which collisions occur, and (assuming you want to make a realistic movie for showing in gnuplot) spit out the state of the simulation at equally spaced times.
// Bonkers1.cc
//   bonkers bonking in a one-dimensional box
//   usage:
//      Bonkers1 [N [dt [tframe(ms) [filename]]]]
//
//   by default, writes to /tmp/1
//   movie can be made by gnuplot
//         load 'gnu10'

#include <iostream>
#include <fstream> 
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <string.h>

using namespace std;
 
struct particle {
  double x    ; // (x,y) coordinates
  double p    ; // momentum
  double im   ; // inverse mass
  double v    ; // velocity
  double T    ; // kinetic energy
  double a    ; // radius of particle
} ; // Note the definition of a structure ends with a semicolon

struct control {
  int verbose    ; // program verbosity
  int printing   ; // period with which to print
  int commanding ; // period with which to issue commands
  int time_in_ms ; // time between frames (commands)
  clock_t latest_time;
  ofstream fout  ; 
};

// Implement     x += v dt for one particle
void PositionStep ( particle &a , double dt )
{ 
  a.x += dt * a.v ; 
}

void v2p( particle &a )
  // propagate the changed velocity into the momentum vector
{
  a.p =  a.v / a.im ;
}

void p2v ( particle &a ) {
  a.v  = a.p * a.im ;
}

void pv2T ( particle &a ) {
  a.T = 0.5*a.v * a.p ; 
}

void collide ( particle & a , particle & b  )
  // two particles collide elastically
{
  // find the relative velocity 
  double velocity_along_line =  a.v - b.v;
  // find a's mass fraction
  double af = a.im / ( a.im + b.im ) ; 
  double bf = b.im / ( a.im + b.im ) ; 
  // reverse the c.o.m. velocity of each along the line of collision
  a.v -= 2.0 * af * velocity_along_line ;
  b.v += 2.0 * bf * velocity_along_line ;
  v2p( a ) ; 
  v2p( b ) ; 
}

// find next collision time 
double nct( particle *a , int NN , int &whichn ) {
  double dt = 1e100 ; 
  // examine all adjacent pairs from 0,1   to  NN-2,NN-1
  for ( int n = 0 ; n < NN-1 ; n ++ ) {
    double relative_velocity = a[n].v - a[n+1].v ;
    if(relative_velocity > 0.0) {
      double collision_time =  ((a[n+1].x-a[n+1].a) - (a[n].x+a[n].a))
	                        /relative_velocity ;
      if ( collision_time < dt ) {
	dt = collision_time ;
	whichn = n ;
      }
    }
  }
  return dt ;
}

void leapForward(  particle *a , int NN , double dt ) {
  for( int n = 0 ; n < NN ; n ++ ) 
    PositionStep( a[n] , dt ) ; 
}

void showState ( particle *a , int n0 ,  int NN, ostream &fout )
{
  for( int n = n0 ; n < NN ; n ++ ) {
      fout << "\t"<<a[n].x;
      fout << "\t"<<a[n].v;
  }
  fout << endl;
}

double kineticEnergy ( particle *a , int NN )
{
  double E = 0.0 ; 
  for( int n=0 ; n < NN ; n ++ ) {
      if ( a[n].im > 0.0 ) // avoid infinitely massive objects
	E+=0.5*a[n].v*a[n].v/a[n].im;
  }
  return E ;
}

void  simulateBonkers( particle *a , int NN ,
		       double &t , double dt , double T ,
		       control &c ) {
  double next_print_dt = 0.0, next_collision_dt ;
  int whichn ;
  int we_are_printing_this_time = 1 ; 
  for(;t<=T;) {
    if( we_are_printing_this_time ) {
      c.fout << t << "\t" << kineticEnergy(a,NN) ;
      showState ( a , 0 , NN , c.fout ) ;
      next_print_dt = dt ;
    }
    // find the next event
    next_collision_dt = nct( a , NN , whichn ) ;
    // ^^ this returns the time, and sets 'whichn'
    if ( next_collision_dt < next_print_dt ) {
      // advance time to that event, have a collision
      leapForward( a ,  NN , next_collision_dt ) ;
      t += next_collision_dt ;
      next_print_dt -= next_collision_dt ;
      collide( a[whichn] , a[whichn+1] ) ;
      we_are_printing_this_time = 0 ; 
    } else {
      leapForward( a , NN , next_print_dt ) ;
      t += next_print_dt ;
      we_are_printing_this_time = 1 ; 
    }
  }
}

int main(int argc, char* argv[])
{
  int N = 4    ; // number of particles
  double T =50.0  ; // target time
  particle   *a ;
  control    c ;
  char filename[256]="/tmp/1" ;
  // set defaults
  double dt = 0.03 ;
  double t = 0.0 ;
  c.verbose = 1 ;
  c.time_in_ms = 100; // real time between plots (in ms)

  // read in any command-line arguments
  if(argc>1)   {
    sscanf( argv[1], "%d", &N ) ; // put the first command-line argument in N
  }
  if(argc>2) {
    sscanf( argv[2], "%lf", &dt ) ; // put the 2nd argument in dt
  }
  if(argc>3) {
    sscanf( argv[3], "%d", &(c.time_in_ms) ) ;
  }
  if(argc>4) {
    strncpy(filename, argv[4], 256 );
  }

  // try to write output file
  c.fout.open(filename);
  if(c.fout.good()==false){
    cerr << "can't write to file " << filename << endl;
    exit(0);
  }
  
  a = new particle[N+2] ; // 0 (wall) , 1...N , N+1 (wall)
  // two walls
  for ( int n=0 ; n <= N+1 ; n += N+1 ) {
    a[n].im = 0.0 ; 
    a[n].v  = 0.0 ;
    a[n].x  = 1.0*n ;
    a[n].a  = 0.0   ; // radius
  }
  // define some masses, positions, and velocities
  for ( int n=1 ; n <= N ; n ++ ) {
    a[n].im    = 1.0/(10.0*n + 1.23)  ; 
    a[n].v     = 1.0*n ;
    a[n].x     = 1.0*n ;
    a[n].a     = 0.1   ; // radius
    cerr<<"mass: "<<1.0/a[n].im<<endl;
  }
  for ( int n=0 ; n <= N+1 ; n ++ ) {
    v2p(a[n]);
  }
  simulateBonkers( a , N+2 , t , dt , T , c ) ; 
  
  return 0;
}

The gnuplot file gnu10
## for use with Bonkers1.cc

## Tip:
## run 
##  gnuplot -noraise
## to avoid window-raising
## or 
##  gnuplot -noraise -persist < gnu10

unset mouse; ## disables the lower-left
         ## display of (x,y) coordinates
         ## which otherwise flicker

t=0; 
T=502 ;
dt=0.1;
set nokey
set origin 0.1,0.2
set pointsize 6 ;
set autos xy
set size 0.8,0.7
print "First showing the Energy \
as a function of time..."
plot '/tmp/1' u 1:2 w l lt 7

set yrange [-5:5]
set xrange [-0.1:5.1]
set size 0.8,0.2
pause 1  ''
pause -1 'Energy - ready'

call 'gnu11'
load 'gnu10') to make a gnuplot movie of the simulation. It works by calling the gnuplot file gnu11
## called by load 'gnu10'
## displays 4 balls in a one-dimensional box

set multiplot
set origin 0.1,0.1
plot\
 '/tmp/1' every 1::t::t u 3:4 w p ps 1 pt 1,\
 '/tmp/1' every 1::t::t u 5:6 w p ps 5 pt 6,\
 '/tmp/1' every 1::t::t u 7:8 w p ps 5 pt 6,\
 '/tmp/1' every 1::t::t u 9:10 w p ps 5 pt 6,\
 '/tmp/1' every 1::t::t u 11:12 w p ps 5 pt 6,\
 '/tmp/1' every 1::t::t u 13:14 w p ps 1 pt 1

 set origin 0.1,0.5
plot\
 '/tmp/1' every 1::t::t u 3:(0.0) w p ps 1 pt 1,\
 '/tmp/1' every 1::t::t u 5:(0.0) w p ps 5 pt 6,\
 '/tmp/1' every 1::t::t u 7:(0.0) w p ps 5 pt 6,\
 '/tmp/1' every 1::t::t u 9:(0.0) w p ps 5 pt 6,\
 '/tmp/1' every 1::t::t u 11:(0.0) w p ps 5 pt 6,\
 '/tmp/1' every 1::t::t u 13:(0.0) w p ps 1 pt 1
 unset multiplot
 
t=t+1;
pause dt;
if(t<T) call 'gnu11'
which then calls itself recursively.


Two-dimensional solution

Bonk.cc is a worked solution for the first task:
...a single function - let's call it collide - which receives two particles with masses m1 and m2 and velocities u1 and u2, and returns the new velocities v1 and v2 of the two particles after an elastic collision between the two.

It is a general worked solution for velocities u1 u2 in any number of dimensions D.
// Bonk.cc -- demonstrates a structure called 'particle'
//            and a function called 'collide', which collides two particles
//            in D dimensions.

#include <iostream>
#include <cmath>
using namespace std;

#define D  2  // number of dimensions
struct particle {
  double x[D] ; // (x,y) coordinates
  double p[D] ; // momentum
  double im   ; // inverse mass
  double v[D] ; // velocity
  double T    ; // kinetic energy
} ; // Note the definition of a structure ends with a semicolon


// Here's a different way passing a structure by reference
void showVelocityC ( particle & a ) { // the "&" indicates 'this is a pass by
                                      //  reference. 
  for ( int i = 0 ; i < D ; i++ ) { 
    cout << "v["<<i<<"] = "<<a.v[i]<<endl; 
  }
  //  cout << "Kinetic energy = " << a.T << endl ; 
}

void makeVelocityC ( particle & a ) {
  a.T=0.0;
  for ( int i = 0 ; i < D ; i++ ) { 
    a.v[i]  = a.p[i] * a.im ;
    a.T    += 0.5*a.v[i] * a.p[i] ; 
  }
}

void v2p( particle & a )
  // propagate a change in velocity into the momentum vector
{
  for ( int i = 0 ; i < D ; i++ ) {
    a.p[i] = (1.0/a.im) * a.v[i] ;
  }
}


void collide ( particle & a , particle & b , double d[D] )
  // two particles collide elastically; the line joining
  // their centres is in direction d, which should be normalized
{
  // find the relative velocity 
  double velocity_along_line = 0.0 ; 
  for ( int i = 0 ; i < D ; i++ ) {
    velocity_along_line += ( a.v[i] - b.v[i] ) * d[i] ;
  }
  if(velocity_along_line<0.0) {
    cerr << "warning: particle collided but velocity "
	 << "along alleged line joining centres "
	 << "is negative: " << velocity_along_line
	 << endl ;
  }
  // find a's mass fraction
  double af = a.im / ( a.im + b.im ) ; 
  double bf = b.im / ( a.im + b.im ) ; 
  // reverse the c.o.m. velocity of each along the line of collision
  for ( int i = 0 ; i < D ; i++ ) {
    a.v[i] -= 2.0 * af * velocity_along_line * d[i] ;
    b.v[i] += 2.0 * bf * velocity_along_line * d[i] ;
  }
  v2p( a ) ; 
  v2p( b ) ; 
}

void checkCollide(  particle & a , particle & b , double d[D] )
{
  // check normalization of d
  double mag = 0.0;
  for ( int i = 0 ; i < D ; i++ ) {
    mag += d[i] * d[i] ;
  }
  if ( mag < 0.99 || mag > 1.01 ) {
    cerr << "Warning, collision direction not normalized\n" ;
    if (mag>0.0) {
      mag = 1.0 / sqrt(mag) ;
      for ( int i = 0 ; i < D ; i++ ) {
	d[i] *= mag ;
      }
    } else {
      cerr << "WARNING, collision direction too small for progress\n" ; 
    }
  }
  // assume v is right and make sure p is consistent
  v2p( a ) ; 
  v2p( b ) ; 

  // check energy and momentum before
  double momentum_before[D], momentum_after[D] , diff[D] ;
  double energy_before=0.0 , energy_after =0.0; 
  for ( int i = 0 ; i < D ; i++ ) {
    momentum_before[i] = a.p[i]+b.p[i] ;
    energy_before += 0.5*a.p[i]*a.v[i] +  0.5*b.p[i]*b.v[i] ; 
  }
  
  collide( a , b , d );

  for ( int i = 0 ; i < D ; i++ ) {
    momentum_after[i] =  a.p[i]+b.p[i] ;
    energy_after += 0.5*a.p[i]*a.v[i] +  0.5*b.p[i]*b.v[i] ;
  }
  if( fabs(energy_after-energy_before) > 1e-4 ) {
    cerr << "Warning - energy mismatch?\n" ; 
  }
  cerr << "\tBefore\tAfter\n" ;
  cerr << "Energy\t" << energy_before 
       << "\t" << energy_after << endl;
  int  mismatch = 0 ;
  cerr << "\tBefore\t\tAfter\n" ;
  for ( int i = 0 ; i < D ; i++ ) {
    diff[i] = momentum_before[i] - momentum_after[i] ;
    if (( diff[i]*diff[i] / ( momentum_before[i] + momentum_after[i] + 1e-5 )
	  / ( momentum_before[i] + momentum_after[i] + 1e-5 ) )>1e-10 ) {
      mismatch ++ ; 
    }
    cerr << "Mom["<<i<<"]\t"<< momentum_before[i] << "\t\t"
	 << momentum_after[i]  << endl ; 
  }
  if ( mismatch ) 
    cerr << "WARNING -  momentum mismatch?\n" ; 
  
}
  

int main()
{
  particle   a ;
  particle   b ; 

  a.im    = 1.0/10 ;
  a.v[0] = 1.0;
  a.v[1] = 0.0;

  b.im    = 1  ;
  b.v[0] = 0.0;
  b.v[1] = 0.0;

   double d[2] = { 0.7 , 0.7 } ;
  //    double d[2] = { 1 , 0 } ;
  cout << "a's velocity before\n" ; 
  showVelocityC( a ) ; 
  cout << "b's velocity before\n" ; 
  showVelocityC( b ) ; 
  checkCollide( a , b , d ) ;
  cout << "a's velocity after\n" ; 
  showVelocityC( a ) ; 
  cout << "b's velocity after\n" ; 
  showVelocityC( b ) ; 

  a.im    = 1.0/10 ;
  a.v[0] = 1.0;
  a.v[1] = 0.0;

  b.im    = 1.0/10 ;
  b.v[0] = 0.0;
  b.v[1] = 0.0;

  cout << "a's velocity before\n" ; 
  showVelocityC( a ) ; 
  cout << "b's velocity before\n" ; 
  showVelocityC( b ) ; 
  checkCollide( a , b , d ) ;
  cout << "a's velocity after\n" ; 
  showVelocityC( a ) ; 
  cout << "b's velocity after\n" ; 
  showVelocityC( b ) ; 
  
  return 0;
}


David MacKay
Last modified: Wed Nov 14 10:11:32 2007