WEEK01: introduction, mathematical software, parallel programming, benchmarking[demon image borrowed from SkiFree, a silly Windows game from the early 90s] This week we got an introduction to the course. Neil went through a number of open-source software tools for symbolic and numerical math work, along with programming, visualization, and 3D graphics. We also touched on benchmarking using an iterative pi calculation program; this is a useful tool for quantitatively comparing the performance of various computing platforms across a range of scales. Our assignment was to model Maxwell's Demon using tools of our choice. Maxwell's Demon is the creation of James Clerk Maxwell; in the thought experiment, a sealed container is divided into two zones separated by a door that can be opened or closed without expending any energy. The container is filled with gas molecules that have different energy levels (speeds). A demon controls the door, only opening it when fast molecules come near and closing it to prevent slow molecules from passing through. By doing so, the container eventually sorts the gas into hot and cold sections, thereby decreasing the overall system entropy and violating the second law of thermodynamics.
- Maxwell's Demon Jupyter Notebook
modelI decided to use Python along with
numpyto run the simulation. I did the development within a Jupyter Notebook for ease of rapid iteration and because it can output static HTML documents for posting. The project was a bit of a learning curve -- I'm used to semicolons rather than tabs (i.e., new to Python)-- but eventually I got the simulation to work.
I created a structured array called
particles that keeps track of each particle's x/y position, previous x/y position, x/y velocity, and color. These were initialized and filled with random values using the
random generator, then color transparency values were set to the velocity vector magnitude. I added a line to represent the demon's door, then created an animation loop that updates particle position and tests for wall and door collisions. Here is that animation playing from initialization through sorting; it's a bit scratchy but you can see the particles sorting themselves out across the door (animation loops every ~25 seconds or so):
The code is fairly simple; figuring out how to test and modify the structured array's values took a bit of time but it worked out:
%matplotlib notebook import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation # fixed random seed np.random.seed(8675309) # create an empty box for particles box = plt.figure(figsize=(6, 6)) ax = box.add_axes([0,0,1,1], frameon=False) ax.set_xlim(0,1), ax.set_xticks() ax.set_ylim(0,1), ax.set_yticks() # create structured array of particles n_particles = 500 particles = np.zeros(n_particles, dtype=[('position', float, 2), ('position_prev', float, 2), ('speed', float, 2), ('color', float, 4)]) # initialize particles: random speed, random position; color = sqrt(xspd^2 + yspd^2) particles['position'] = np.random.uniform(0,1,(n_particles, 2)) particles['speed'] = np.random.uniform(-0.005,0.005, (n_particles, 2)) particles['color'][:, 3] = np.sqrt(np.add(np.square(particles['speed'][:, 0]), np.square(particles['speed'][:, 0]))) * 100 # make scatter plot which will be updated scat = ax.scatter(particles['position'][:,0], particles['position'][:,1], facecolors=particles['color']) # add demon gate gate = ax.plot([.5, .5],[0,1]) def update(frame_number): # update particle position based on speed and direction particles['position'][:, 0] += particles['speed'][:, 0] particles['position'][:, 1] += particles['speed'][:, 1] for p in particles: # check for wall collisions and bounce if p['position'] > 1 or p['position'] < 0: p['speed'] *= -1 if p['position'] > 1 or p['position'] < 0: p['speed'] *= -1 # check for speed at gate and bounce if p['position'] > 0.5 and p['position_prev'] < 0.5 and p['color'] < 0.35: p['speed'] *= -1 if p['position'] < 0.5 and p['position_prev'] > 0.5 and p['color'] >= 0.35: p['speed'] *= -1 # update previous particle positions for gate check p['position_prev'] = p['position'] # update scatter plot scat.set_offsets(particles['position']) # make the animation using the update function as the director animation = FuncAnimation(box, update, interval=20) plt.show