Unit 3: Simulations

Lesson 21: Equilibrium part 2

Notebook Authors

(fill in your two names here)

Team Roles

Facilitator: (fill in name)
Spokesperson: (fill in name)
Process Analyst: (fill in name)
Quality Control: (fill in name)

If there are only three people in your team, have one person serve as both spokesperson and process analyst for the rest of this activity.

At the end of this Lesson, you will be asked to record how long each Model required for your team. The Facilitator should keep track of time for your team.

Computational Focus: Simulation Statistics

In this Lesson, we'll continue to work on the Gas Diffusion model that we started last Lesson and think more about equilibrium.

Model 1: Monte Carlo Algorithm for Gas Diffusion

In the previous lesson, we selected a particle at random and moved it to the other side. This procedure, however, is cumbersome, because we need to keep track of the position of all $N$ particles. We can simplify our algorithm by realizing our only interest is the number of particles on each side. Actually, we only need to know $n_{left}$, the number of particles on the left side, since the number on the right side is the difference, $N−n_{left}$. Because each particle has the same likelihood of going through the hole, the probability per unit time that a particle moves from left to right equals the number of particles on the left side divided by the total number of particles, so the probability of a move from left to right is: $$\frac{n_{left}}{N}$$

Check your code:

At this point, if your RBox class and methods are all functioning correctly, and you can easily get the values for the total number of particles and the number of particles on the left side, you can move along. IF your RBox class and methods are all messed up and you'll have a hard time getting the info you need out of them, then you should probably rewrite the RBox class using the "counter" method described in general above, and in more detail below:

The algorithm for simulating the evolution of the model is given by the following steps:

  1. Generate a random number $r$ from a uniformly distributed set of random numbers in the unit interval $0 \le r < 1$.
  2. If $r < \frac{n_{left}}{N}$, increase the number of particles in the right box by one, otherwise increase the number of particles in the left box by one.
  3. Increase the time by one time unit.

Make sure that you have viable code before moving on! Check with an instructor if you have questions.

Critical Thinking Questions

1. Calculate and Explain the probability that a particle in the left box will move to the right box for each of the following (show your work, $\LaTeX$ is good for math):

1a. all the particles are in the left box

1b. all the particles are in the right box

1c. the number of particles in each box is equal

2. If there are more particles on the left side of the box then the right, describe how the number of particles in each box will change. Justify your answer (and it better have probabilistic language in it).

3. If $r > \frac{n_{left}}{N}$ what is the new value of $n_{left}$?

4. Paste your RBox code below. Fix it and/or make sure it works.


In [ ]:

4a. Examine your RBox code. Describe each of your instance variables, what they are they for? Both answer here and comment your code!


In [ ]:

5. Make sure that your RBox code uses the probability of a random particle from the left being chosen as part of your algorithm (either directly or indirectly). Modify it so that it calculates the probability of a random particle from the left being chosen, it can print that probability until you know it works, then comment it out (you'll use it later for the history). Explain here and comment your code!


In [ ]:

Team Programming part 1

Make sure that as you work on these programming tasks, you are using your best pair programming technique (both driver and navigator) and switching drivers frequently.

Your code must have docstrings and comments for you to receive full credit.

We are still using an Object Oriented Programming approach for this Lesson.

A. Copy and paste your working RBox class here and rename it RandomBox double-check that all of the methods still work.


In [ ]:

B. Modify your run_simulation method so that it creates a record of the entire simulation, by saving the number of particles in the left half of the box at each time unit of the simulation to a list. In other words, the elements in the list will be the number of particles in the left half at the start of the simulation, after the first time unit, after the second time unit, etc. This list should be an instance variable in your RandomBox class. Document a few brief tests.

C. Add a get_history method that takes no other parameters besides self. It should return the history instance variable (the list recording the number of left particles at each time unit of the simulation). If the simulation has not been run, there should only be one number in the list. Document a few brief tests.

D. Modify your run_simulation method to detect when the simulation has reached equilibrium and store the value of the time step, but it should not stop. Again, you can have it print this information for now. Document a few brief tests.

optional. You can code a method to plot the history to get a visual sense of how your simulation has gone. It is also nice to have a horizonal line at equilibrium.

Model 2: Equilibrium Fluctuations

You should recall from other science classes that equilibrium is not static, rather there are still changes in the system over time. These changes can create fluctuations around the equilibrium point.

In the individual assignment for previous lesson, you have analyzed or will visually analyze line plots of the time evolution of $n_{left}(t)$ (where $t$ is time). The simulation could be described as having two different portions: the simulation before the system reached equilibrium, and after the system reached equilibrium.

A more quantitative measure of the equilibrium fluctuations is the mean square fluctuations, $\Delta n^2$, which is defined as: $$\Delta n^2 = \langle n^2 \rangle - \langle n \rangle ^2$$
The angle brackets, $\langle \rangle$, denote an average taken after the system has reached equilibrium. And the $n$'s are $n_{left}(t)$.

The relative magnitude of the fluctuations is a better measure for comparing simulations of different numbers of particles. The relative magnitude of the fluctuations is:
$$\frac{\Delta n}{ \langle n \rangle}$$
where $\Delta n$ is the square root of the mean square fluctuations and $n=n_{left}(t)$ -- the number of particles in the left half of the box at a given time, $t$, during the simulation. In this case, it would be the average, since it is in angle brackets, $\langle \rangle$.

Critical Thinking Questions

6. What is another name for a similar calculation to "mean square fluctuations" that we used to describe data in our earlier Statistics Lesson?

7. What is another name for a similar calculation to "square root of the mean square fluctuations" that we used to describe data in our earlier Statistics Lesson?

8. When calculating the equilibrium statsitics above, do you think you should include all values of $n_{left}(t)$ starting at $t=0$? Why or why not?

9. Propose an appropriate Boolean condition that might be used to signal when the program should start calculating the desired equilibrium statistics.

10. The run_simulation method should execute the simulation until the time when meaningful statistics can be calculated. There are multiple ways of modifying the RandomBox class so it calculates the averages, the mean square fluctuations, and the relative magnitude of the fluctuations after an initial equilibration period. To make sure everyone has similar code for the next model (since it will matter), you should plan to modify the code by introducing two sequential loops (not nested loops).
10a. Generally describe how long the first loop that does not calculate statistics will execute.

10b. Generally describe how long the second loop that calculates statistics will execute.

Team Programming part 2

Make sure that as you work on these programming tasks, you are using your best pair programming technique (both driver and navigator) and switching drivers frequently.

Your code must have docstrings and comments for you to receive full credit.

F. Modify the run_simulation method so it stops after equilibrium has been reached or the simulation time has ended, whichever occurs first. If equilibrium has been reached, print a message reporting when equilibrium was reached and if equilibrium has not been reached print a different, approproate message. Document a few brief tests.

G. Once you are satisfied that your simulation code correctly identifies when equilibrium has been reached, add the second loop that continues to run the simulation after equilibrium. Don’t worry about keeping track of statistics yet. You should make sure your method still runs the simulation for the expected total number of time units. Document a few brief tests.

Model 3: Helper Methods

After adding the second loop, you may notice that some duplicate code in your run_simulation method. Good programmers try to avoid duplicate code by introducing helper methods to their programs. The duplicate code is moved to a helper method, and the other code makes multiple calls to the helper method.

Using helper methods have two advantages:

  • Corrections and modifications are easier when they only have to be made once to the helper method, rather than multiple times in a larger program.
  • Long or complex code is easier to read and debug when some of the code is moved into helper methods.

To keep the program as readable as possible, helper methods should be given names that describe the functionality of the code.

Moving code that references or modifies instance variables to helper methods require no modifications. However, moving code that references local variables requires passing the local variables as parameters to the helper method. Moving code that modifies local variables requires both passing the local variables as parameters and returning the new values as return values. For this reason, it is sometimes not worthwhile to move duplicate code that modifies a local variable to a helper method.

Critical Thinking Questions

11. Examine your run_simulation method and identify the duplicate code. What instance variables are referenced and/or modified by the duplicate code?

12. What local variables are referenced and/or modified by the duplicate code?

13. As a team, agree on the lines of duplicate code that should be moved to a helper method. What would be a descriptive name for that method (based on the functionality of the code)?

14. Write high level pseudocode for the helper method for the duplicate code from the question above, and call the helper method from inside both loops

Pseudocode for Statistics

Write high level pseudocode for the simulation statistics in 2 steps as detailed below
(the individual assignment will be to write the Python code).

Review the equation for the mean square fluctuations and write pseudocode to keep track of the running totals.

Finally, write pseudocode to calculate the averages, the mean square fluctuations, and the relative magnitude of the fluctuations.
Since this is a lot of code, we’ll end up putting it in another helper method (called calculate_statistics). The run_simulation method should call calculate_statistics, passing it all the values that it has tracked during the simulation.
The calculate_statistics method should have multiple parameters, and it should use those values to calculate the averages, mean square fluctuations, and the relative magnitude of the fluctuations.

Temporal Analysis Report

How much time did it require for your team to complete each part?

Model 1:

Team Programming 1:

Model 2:

Team Programming 2:

Model 3:

Pseudocode for Statistics: