In this simulation, we consider a box divided in half with by a membrane with a small hole in the center. The density of the gas is small enough so that the gas molecules do not interact. The total number of particles in the box is N, and the number on one side starts at n (leaving N-n on the other side).
This simulation implements the ApproachToEquilibriumParticlesInBox() algorithm described on pages 9-10 of Statistical and Thermal Physics by Gould and Tobochnik.
from __future__ import division #crucial to do this first!
First I define two functions: (1) propagate(n,N) generates a random number between 0 and 1 and compares it to the ratio on n/N (number on left/total number) and either boosts the number on the left by one or decreases it by one. (2) createRun(n, N, timeSteps) generates a complete run with length = timeSteps. At it's conclusion, it returns two arrays: left() = number of particles on left side of box at each time step, and right() = number of particles on right side of box at each time step.
def propagate(n,N): r = random.random() if r < n/N and n > 0: n = n - 1 elif r >= n/N and n>=0: n = n + 1 elif r < n/N and n == 0: n = n else: print "problem: should not have arrived here!" pass return n def createRun(n, N, timeSteps): t = np.arange(0,timeSteps) left =  right = for i in range(len(t)): n = propagate(n, N) left.append(n) right.append(N-n) return np.array(left), np.array(right)
Now I define a function runAndPlot() which creates a simulation with the option to run a number of different trials (numTrials); the routine then averages the results at each time step so that one can get an "ensemble" average (not really, since the ensemble average requires an infinite number of trials.)
def runAndPlot(nleft, N, timeSteps, numTrials): for i in range(numTrials): left, right = createRun(nLeft, N, timeSteps) if i==0: numLeft = left numRight = right if i >= 0: numLeft = np.vstack((numLeft, left)) numRight = np.vstack((numRight, right)) figsize(10,8) xlabel("Number of Steps") ylabel("Number on Left Half of Box") ylim(0,N) plot(np.average(numLeft,axis=0))
Here are three differnet runs with 1, 10, and 1000 trials all displayed together:
runAndPlot(1,10,100,1) runAndPlot(1,10,100,10) runAndPlot(1,10,100,1000)