First part:

Write a function

function [ v1, v2 ] = collide(  u1, u2 , im1, im2  ) 

Work out what this function should be. Think about special cases, then generalize. Check your answer on the special cases. Check that your general answer does indeed conserve momentum and energy.

Solution: collide.m, and a testing script: CHECKcollide.m.

function [ v1, v2 ] = collide(  u1, u2 , im1, im2  ) 
  ## returns the new velocities of two particles 
  ## with inverse-masses im1 and im2 
  ## after an elastic collision. 
  ##
  ## Dimensions: u1,u2:   [L/T]
  ##             im1,im2: [1/M]
  ##             v1,v2:   [L/T]

  ## compute centre of mass velocity (see proof below)
  c         = ( im2 * u1 + im1 * u2  ) / ( im1 + im2 ) ;
  ##  [L/T] = ( [L/T/M]  +  [L/T/M]  ) / ( [1/M]+[1/M] )  

  ## in the CM frame, velocities reverse
  v1       = c - ( u1 - c ) ;
  v2       = c - ( u2 - c ) ;
  ## [L/T] = [L/T] - [L/T] + [L/T] 
  
endfunction

## See CHECKcollide for a script that checks various special cases.

## centre of mass velocity:
## m_TOT c = m1 u1 + m2 u2
##  c = ( m1 u1 + m2 u2 ) / ( m1 + m2 )
##    =  (im2 u1 + im1 u2) / (m1+m2)/(m1m2)
##    =  (im2 u1 + im1 u2) / (im1 + im2 ) 

Second part:

Simulate a collection of N particles with different masses, bouncing off each other, and moving in one dimension. For the final part of this project, it's important that the particles should all have different masses (for example, ranging between 1 and 9 mass units). But initially when developing your program it might be wise to check what happens in the special case where they are equal.

Collect data at every time Tc, 2Tc, 3Tc, ... If you are making a movie, you will want to have Tc a little smaller than the typical time between collisions. If you are just gathering statistics and not making a movie, you can use a bigger Tc -- perhaps N times bigger is reasonable, where N is the number of particles.

Run the simulation for a long time, and then plot histograms of the velocities and positions of all the different masses. Pick two particles, having the biggest and smallest masses, and look at the histograms of their velocities. Make a theory to describe these histograms.

Here are histograms of velocities using 400 samples. The histogram with velocities ranging from -8 to 8 is for a particle with mass 9. The histogram with velocities ranging from -24 to 24 is for a particle with mass 1.

To make good histograms you will want 1000s of samples. When making histograms, you may find the octave function histo useful. However, in my experience, this function does not do what I really want. Its bin sizes are too big. I therefore offer to you my histogram routine, djcmhisto.m, which also returns a Gaussian Fit to the histogram, and tells you the mean and standard deviation of each column.

You can see a simple example use of this histogram program here (TESTdjcmhisto.m).

Solution

A solution is in bounce1.m.

Here are examples of the types of graph that are interesting to plot. In order to not give away the results that you will get for a simulation with lots of masses, with masses ranging from 1 to 9, I show here some small histograms and scatterplots for a system with just two balls.

(1) (most important): Histogram of velocity of each particle

(2) scatter plot of positions of particles, or histograms of positions.

(3) scatter plot of position versus velocity.

bounce1.m is invoked by RUNbounce0.m

## Given a state (x,v) of length N and inverse-masses im(1..N)
## where x(1) = left-hand wall  and im(1) = 0 
## and   x(N) = right hand wall and im(N) = 0 

## determine which of the inter-particle collisions will take place next
## -- or is the next event one of the data-collection times {Tc, 2 Tc, 3Tc...}? --
## and simulate the motion until that event, and implement the collision. 

time_to_next_clock = Tc ;

j = 0 ; ## counter: will run up to Number_samples_required 
logbook = zeros(Number_samples_required , 2*N ) ; 

########################################################################
## set up graph
figure(1);
y = zeros(size(x)) ;
## look at the current velocities to guess a good range for the velocity graph
vmax=max(v);
vmin=max(-v) ;
vmax=max(vmax,vmin);
axis([x(1)-1 x(N)+1 -vmax*5.8 vmax*5.8]); 
gset border 15 
gset pointsize 3
########################################################################

################################################################################
## Main loop: 
##             Keep simulating the motion until we have all the data we want
while ( j < Number_samples_required ) 

  ## run through all (left-hand) particles finding the time until collision with their 
  ## right neighbours.
  ## Initialize the collisiontime vector to large positive numbers. 
  collisiontime=ones(1,N) * 2 * time_to_next_clock ;
  ## The Nth particle does not have a right neighbour, so we can use its "collision time" 
  ## to store the time of the next clock event, which is when we will gather data for histograms
  collisiontime(1,N) =  time_to_next_clock ; ## this is used to indicate a clock event
  for n=1:N-1
    dv = v(1,n) - v(1,n+1) ;
    ## if the right neighbour is coming towards us (in our frame of reference)
    ## then work out when the collision would be
    if ( dv > 0.0 ) 
      collisiontime(n) = (x(1,n+1)-x(1,n))/dv ; 
    endif
  endfor

  ## decide which collision needs to be rendered next
  [a,b] = min( collisiontime ) ;

  ## simulate the motion for a time given by the min
  dt = a ; 
  x = x + v * dt ;
  time_to_next_clock = time_to_next_clock - dt ;
  
  left = b ;  ## b holds the index of the particle that has its collision next
  if( left < N ) 
    right = b+1;
    ## change the velocities of the dudes involved in a collision
    [ v(1,left), v(1,right) ] = collide(  v(1,left), v(1,right) , im(1,left), im(1,right) ) ;
##    plot( x , y , "x6" , x , v , "o3");
  else
    ## no collision, it's clock time, so time to take a sample!
    j ++ ;
    if(doplot) 
      plot( x , y , "x5" , x , v , "o3");
      pause(0.075);
    endif
    time_to_next_clock = Tc ;
    logbook(j,:) = [x,v ] ; 
  endif
endwhile

bounce1.m is invoked by RUNbounce0.m

RUNbounce0.m

x = [ 1:10 ] ; N = 10 ; 
v = zeros(size(x) ) ;
im = zeros(size(x) ) ;
im(2:5) = 1.0/9.0 * [1.1 1.0 0.9 1.0 ] ;
im(6:9) = 1.0     * [1.1 1.0 0.9 1.0 ] ;

v(2:5) = 2.0  ; 
v(3) = 2.21 ; 
v(4) = 2.621 ; 
v(6:9) = -0.2 ; 
v(7) = -0.27 ;
v(8) = -0.37 ;
## to make a nice brief movie use this
Tc = 0.030123091;  Number_samples_required = 1000 ;  doplot = 1 ; 

bounce1;

## see RUNbounce1.m for a run that gives nice histograms

Here's a slightly fancier solution in which the particles have a finite radius, to make better looking movies.
bounce2.m is invoked by cradle2.m

more off
## Given a state (x,v) of length N and inverse-masses im(1..N)
## where x(1) = left-hand wall  and im(1) = 0 
## and   x(N) = right hand wall and im(N) = 0 

## determine which of the inter-particle collisions will take place next
## -- or is the next event one of the data-collection times {Tc, 2 Tc, 3Tc...}? --
## and simulate the motion until that event, and implement the collision. 

time_to_next_clock = Tc ;

j = 0 ; ## counter: will run up to Number_samples_required 
logbook = zeros(Number_samples_required , 2*N ) ; 

########################################################################
## set up graph
figure(1);
y = zeros(size(x)) + ytweak ;
## look at the current velocities to guess a good range for the velocity graph
vmax=max(v);
vmin=max(-v) ;
vmax=max(vmax,vmin);
axis([x(1) (x(N)+xtweak(N)) -vmax*5.8 vmax*5.8]);
gset border 15 
gset pointsize 8

gset nokey

########################################################################

## initialize lastplottime so that we think it is time to plot right away
lastplottime = time - time_between_plots ; 

################################################################################
## Main loop: 
##             Keep simulating the motion until we have all the data we want
while ( j < Number_samples_required ) 

  ## run through all (left-hand) particles finding the time until collision with their 
  ## right neighbours.
  ## Initialize the collisiontime vector to large positive numbers. 
  collisiontime=ones(1,N) * 2 * time_to_next_clock ;
  ## The Nth particle does not have a right neighbour, so we can use its "collision time" 
  ## to store the time of the next clock event, which is when we will gather data for histograms
  collisiontime(1,N) =  time_to_next_clock ; ## this is used to indicate a clock event
  for n=1:N-1
    dv = v(1,n) - v(1,n+1) ;
    ## if the right neighbour is coming towards us (in our frame of reference)
    ## then work out when the collision would be
    if ( dv > 0.0 ) 
      collisiontime(n) = (x(1,n+1)-x(1,n))/dv ; 
    endif
  endfor

  ## decide which collision needs to be rendered next
  [a,b] = min( collisiontime ) ;

  ## simulate the motion for a time given by the min
  dt = a ; 
  x = x + v * dt ;
  time_to_next_clock = time_to_next_clock - dt ;
  
  left = b ;  ## b holds the index of the particle that has its collision next
  if( left < N ) 
    right = b+1;
    ## change the velocities of the dudes involved in a collision
    [ v(1,left), v(1,right) ] = collide(  v(1,left), v(1,right) , im(1,left), im(1,right) ) ;
##    plot( x , y , "x6" , x , v , "o3");
  else
    ## no collision, it's clock time, so time to take a sample!
    j ++ ;
    if(doplot) 
      xyv = [ (x+xtweak)' ,y' ,v' ] ;
      ## elegant pausing rule that makes frame rate steady
      timetillnextplot = lastplottime + time_between_plots - time ;
      if ( timetillnextplot > 0.0 ) 
	pause( timetillnextplot ) ;
	timetillnextplot = 0.0 ; 
      else
	## we are behind by the negative quantity timetillnextplot!
	## to try to catch up, we will add this to "lastplottime" later
      endif
      gplot xyv u 1:2 w p 5 6, xyv u 1:3 w p 6 2 
      lastplottime = time + timetillnextplot ; 
    else
      if( rem(j,1000) == 0 )
	j
      endif
    endif
    time_to_next_clock = Tc ;
    logbook(j,:) = [x,v ] ; 
  endif
endwhile

bounce2.m is invoked by cradle2.m

cradle2.m

x = [ 0   14 15 16  17 18 19  20 21     35 ] ; N = 10 ; 
## extra displacement to add to the particles so they appear to bounce
## off each other
xtweak = [ 0 : (N-1) ] * 5 ; 
ytweak = 6;

## in cradle2 the masses are 4,9,1,1,1,1,1,1,1
## and the initial situation has all the energy in the left two masses
v = zeros(size(x) ) ;
im = zeros(size(x) ) ;
im(2:9) = 1.0 ;
im(2) = 1.0/4.123 ;
im(3) = 1.0/9.0123 ;

v(2:3) = - 2.0  ; 

## this specifies how often in simulation time we make a plot
Tc = 0.2030123091;  Number_samples_required = 30000 ;  doplot = 1 ; 
## this specifies how often in real time the plots will appear
time_between_plots = 0.025;

bounce2;