Results from modeling 1D collisions

We have a one-dimensional box with particles in it. They follow the rules for collisions we derived previously. If the masses are all the same, the particles simply exchange velocities, and that isn't very interesting. So our model has particles of two different masses (M/m = 2, 3 or 4 in the simulations I've run). The question is, what is the long-run distribution of velocities in the population of particles?

At the end of this post, I have listed the code I wrote. But I'm not going to go through it, for two reasons. First, I know you would be bored (me too). Second, it isn't very polished, and has a few issues. But the results give a picture of the distribution that I think must be correct. Maybe one of you can tell me if that's true or not.

We can run the simulation with the option to report lots of information as it runs. I always use the same names for common variables. I use v for "verbose", L for "list", D for "dictionary", and so on.

Here is some output from the code listing below with v=True. We show the particle positions in 40 bins, updating only when positions change, and print detailed info when something more interesting happens (reversal, collision). The particles are numbered and identified with the mass type 'm' or 'M'.

62 .................4...1....2.....3.....5. 64 .................4...1....2....3......5. 65 .................4....1....2...3......5. 67 ..................4...1....2...3......5. 68 ..................4...1....2..3.......5. 70 ..................4....1....2.3.......5. 72 ...................4...1....2.3.......5. 73 ...................4...1....23........5. 74 ...................4...1....23.........5 75 ...................4....1....*.........5 collide (i = 77) {2M x = 29.3763 v = 0.2} {3M x = 29.0261 v = -0.2444} {2M x = 29.3763 v = -0.2444} {3M x = 29.0261 v = 0.2} 77 ....................4...1....*.........5 78 ....................4...1...23.........5 80 ....................4....1..23.........5 81 ....................4....1..2.3........5 82 .....................4...1.2..3........5 85 .....................4....12..3........5 86 .....................4....*....3.......5 87 ......................4...*....3.......5 collide (i = 89) {1m x = 26.792 v = 0.2} {2M x = 26.4429 v = -0.2444} {1m x = 26.792 v = -0.3926} {2M x = 26.4429 v = 0.0519} 91 ......................4..12.....3......5 92 .......................4.12.....3......5 93 .......................41.2.....3......5 96 .......................*..2......3.....5 41235 14235 collide (i = 97) {1m x = 23.6513 v = -0.3926} {4M x = 23.8064 v = 0.2} {1m x = 23.6513 v = 0.3975} {4M x = 23.8064 v = -0.1951} 97 .......................41.2......3.....5 14235 41235 99 .......................41..2.....3.....5 100 .......................4.1.2.....3.....5 101 ......................4..1.2......3....5

How do we test the code? This is a hard problem. To make a start, we can run a simulation with X (box size) and the number of particles very small, and look carefully whenever things change. Two other tests: in a 1D box, the order of particles should never change, and also the total kinetic energy should not change with time. The latter requirement is met by this code. However, the former invariant is occasionally broken. The code reports this when it occurs. In the line that starts "96" above, we see "41235 14235", which indicates that the order switched as shown. On the following cycle it switched back, but this doesn't always happen. It gets harder when there are 3 or more particles within the same small interval. I have not solved this problem, but these "quantum"-like events are rare, and I believe they don't invalidate the statistic we are going to measure. (I can't prove that, of course, but the numbers are quite small).

I've run the simulation with more particles and a slightly larger space to make the velocities change more frequently. For these tests, I assigned the same velocity to every particle. I want to show you statistics for these runs, with various mass ratios as mentioned above. Two things are of interest: the standard deviation of the distribution of velocities, and the histograms. Here are the results (in R) for the stdevs:

#M = 2m > sd(m) [1] 0.2467407 > sd(M) [1] 0.1719307 #M = 3m > sd(m) [1] 0.2849638 > sd(M) [1] 0.1620676 #M = 4m > sd(m) [1] 0.3177323 > sd(M) [1] 0.1573599

I will have to look up the values that statistical mechanics theory predicts. It seems like a good guess that the sd of the small particles goes like the sd of the large particles * √M/n. But, I would have to do some other tests to be sure. For example, we have a number of m and M particules that is random, but we did random.seed(157), so it is actually set to one specific value. Finally, notice that the distributions are all normal. I hope this does not come as a surprise.

And here are the plots.
M/m = 2


M/m = 3


M/m = 4

import random,sys class Particle: n = 0 X = None def __init__(self,X): Particle.n += 1 self.n = Particle.n self.kind = random.choice('Mm') self.x = random.random() * X self.v = 0.2 self.next = self.x + self.v if not Particle.X: Particle.X = X def __repr__(self): s = '{' + str(self.n) + self.kind + ' x = ' s += str(round(self.x,4)).rjust(7) + ' v = ' s += str(round(self.v,4)).rjust(7) + '}' return s.ljust(25) def randomSpeed(self,Z=10): s = random.random() / Z if random.random() > 0.5: s *= -1 return s def move(self,i): self.x += self.v next = self.x + self.v if next <= 0: if v: print 'reverse (i =', print str(i) + ')' print self self.x = abs(next) self.v *= -1 if v: print self if next >= X: if v: print 'reverse (i = ', print str(i) + ')' print self self.x = X - (next - X) self.v *= -1 if v: print self #--------------------------------------- def initBox(N,X): def add(L): def unique(L,p): for q in L: if int(q.x) == int(p.x): return False return True p = Particle(X) if unique(L,p): L.append(p) # backtrack on the count of p's else: Particle.n -= 1 L = list() while len(L) < N: add(L) return L def boxRep(L): L2 = ['.'] * X for p in L: x = int(p.x) if L2[x] == '.': L2[x] = str(p.n) else: L2[x] = '*' return ''.join(L2) #--------------------------------------- def collide(p1,p2,i): global v if v: print 'collide (i = ', print str(i) + ')' print p1, p2 if p1.kind == p2.kind: # exchange velocities p1.v,p2.v = p2.v,p1.v p1.next = p1.x + p1.v p2.next = p2.x + p2.v if v: print p1, p2 return switch = p1.kind == 'm' if switch: p1,p2 = p2,p1 a = p1.v b = p2.v # useful fractions are cached f13 = 1.0/3 f23 = 2.0/3 f43 = 4.0/3 f25 = 2.0/5 f35 = 3.0/5 f85 = 8.0/5 # the mass ratio M/m is built into this: r = 2 if r == 2: p1.v = a*f13 + b*f23 p2.v = -b*f13 + a*f43 if r == 3: p1.v = a*0.5 + b*0.5 p2.v = -b*0.5 + a*1.5 if r == 4: p1.v = a*f35 + b*f25 p2.v = -b*f35 + a*f85 p1.next = p1.x + p1.v p2.next = p2.x + p2.v if v: if switch: print p2, p1 else: print p1, p2 def handleCollisions(L,roundNum): # only allow one collision per round alreadyMoved = [] for i,p2 in enumerate(L[:]): for j in range(i): if i in alreadyMoved: continue if j in alreadyMoved: continue p1 = L[j] fromLeft = p1.x < p2.x r = round(p1.next,5) if r == round(p2.next,5): collide(p1,p2,roundNum) alreadyMoved.extend([i,j]) if fromLeft: if p1.next > p2.next: collide(p1,p2,roundNum) alreadyMoved.extend([i,j]) else: if p1.next < p2.next: collide(p1,p2,roundNum) alreadyMoved.extend([i,j]) #--------------------------------------- def order(L): def f(p): return p.x rL = sorted(L,key=f) rL = [str(e.n) for e in rL] if len(rL) < 10: return ''.join(rL) return ','.join(rL) #======================================= random.seed(157) X = 40 v = True L = initBox(5,X) oL1 = order(L) oLp = oL1 prev = boxRep(L) N = 1000 T = N/10 vL = list() def getKE(t): mL,ML = t s = [0.5 * v**2 for v in mL] S = [v**2 for v in ML] return sum(s) + sum(S) for i in range(N): handleCollisions(L,i) for p in L: p.move(i) mL = list() ML = list() if i > T: for p in L: if p.kind == 'M': ML.append(p.v) else: mL.append(p.v) vL.append((mL,ML)) curr = boxRep(L) oL2 = order(L) if not v: continue if curr != prev: print str(i).rjust(5) + ' ', print curr, if oLp != oL2: print oLp, oL2 else: print oLp = oL2 prev = curr print oL1 print order(L) KEo = getKE(vL[0]) print 'KEo', KEo for mL,ML in vL[1:]: KE = getKE((mL,ML)) if KE != KEo: KEo = KE print KEo FHm = open('m.txt','w') # some issue with M and m?? FHM = open('MM.txt','w') for t in vL: mL = [str(round(f,4)) for f in t[0]] ML = [str(round(f,4)) for f in t[1]] FHm.write('\n'.join(mL)+'\n') FHM.write('\n'.join(ML)+'\n') FHm.close() FHM.close() ''' setwd('Desktop') m=read.table('m.txt',head=F)[,1] M=read.table('MM.txt',head=F)[,1] sd(m) sd(M) par(mfrow=c(1,2)) hist(m,freq=F,xlim=c(-1,1),ylim=c(0,2.5), breaks=seq(-1.5,1.5,by=0.1),main='') plot(function(x) dnorm(x,0,sd(m)), min(m),max(m),lwd=3, add=T,col='red') hist(M,freq=F,xlim=c(-1,1),ylim=c(0,2.5), breaks=seq(-1.5,1.5,by=0.1),main='') plot(function(x) dnorm(x,0,sd(M)), min(M),max(M),lwd=3, add=T,col='red') '''