COS 226 Programming Assignment

Molecular Dynamics Simulation of Hard Spheres

Simulate the motion of N colliding particles according to the laws of elastic collision using event-driven simulation. Such simulations are widely used in molecular dynamics (MD) to understand and predict properties of physical systems at the particle level. This includes the motion of molecules in a gas, the dynamics of chemical reactions, atomic diffusion, sphere packing, the stability of the rings around Saturn, the phase transitions of cerium and cesium, one-dimensional self-gravitating systems, and front propagation. The same techniques apply to other domains that involve physical modeling of particle systems including computer graphics, computer games, and robotics.

Hard sphere model. The hard sphere model (billiard ball model) is an idealized model of the motion of atoms or molecules in a container. In this assignment, we will focus on the two-dimensional version called the hard disc model. The salient properties of this model are listed below.

• N particles in motion, confined in the unit box.
• Particle i has known position (rx_i, ry_i), velocity (vx_i, vy_i), mass m_i, and radius σ_i.
• Particles interact via elastic collisions with each other and with the reflecting boundary.
• No other forces are exerted. Thus, particles travel in straight lines at constant speed between collisions.

This simple model plays a central role in statistical mechanics, a field which relates macroscopic observables (e.g., temperature, pressure, diffusion constant) to microscopic dynamics (motion of individual atoms and molecules). Maxwell and Boltzmann used the model to derive the distribution of speeds of interacting molecules as a function of temperature; Einstein used it to explain the Brownian motion of pollen grains immersed in water.

Simulation. There are two natural approaches to simulating the system of particles.

• Time-driven simulation. Discretize time into quanta of size dt. Update the position of each particle after every dt units of time and check for overlaps. If there is an overlap, roll back the clock to the time of the collision, update the velocities of the colliding particles, and continue the simulation. This approach is simple, but it suffers from two significant drawbacks. First, we must perform N² overlap checks per time quantum. Second, we may miss collisions if dt is too large and colliding particles fail to overlap when we are looking. To ensure a reasonably accurate simulation, we must choose dt to be very small, and this slows down the simulation.
• Event-driven simulation. With event-driven simulation we focus on those times at which interesting events occur. In the hard disc model, all particles travel in straight line trajectories at constant speeds between collisions. Thus, our main challenge is to determine the ordered sequence of particle collisions. We address this challenge by maintaining a priority queue of future events, ordered by time. At any given time, the priority queue contains all future collisions that would occur, assuming each particle moves in a straight line trajectory forever. As particles collide and change direction, some of the events scheduled on the priority queue become "stale" or "invalidated", and no longer correspond to physical collisions. We can adopt a lazy strategy, leaving such invalidated collision on the priority queue, waiting to identify and discard them as they are deleted. The main event-driven simulation loop works as follows:
  • Delete the impending event, i.e., the one with the minimum priority t.
  • If the event corresponds to an invalidated collision, discard it. The event is invalid if one of the particles has participated in a collision since the event was inserted onto the priority queue.
  • If the event corresponds to a physical collision between particles i and j:
    • Advance all particles to time t along a straight line trajectory.
    • Update the velocities of the two colliding particles i and j according to the laws of elastic collision.
    • Determine all future collisions that would occur involving either i or j, assuming all particles move in straight line trajectories from time t onwards. Insert these events onto the priority queue.
  • If the event corresponds to a physical collision between particles i and a wall, do the analogous thing for particle i.

  This event-driven approach results in a more robust, accurate, and efficient simulation than the time-driven one.

Collision prediction. In this section we present the formulas that specify if and when a particle will collide with the boundary or with another particle, assuming all particles travel in straight-line trajectories.

• Collision between particle and a wall. Given the position (rx, ry), velocity (vx, vy), and radius σ of a particle at time t, we wish to determine if and when it will collide with a vertical or horizontal wall.Since the coordinates are between 0 and 1, a particle comes into contact with a horizontal wall at time t + Δt if the quantity ry + Δt vy equals either σ or (1 - σ). Solving for Δt yields:
  

  

  An analogous equation predicts the time of collision with a vertical wall.

• Collision between two particles. Given the positions and velocities of two particles i and j at time t, we wish to determine if and when they will collide with each other.Let (rx_i' , ry_i' ) and (rx_j' , ry_j' ) denote the positions of particles i and j at the moment of contact, say t + Δt. When the particles collide, their centers are separated by a distance of σ = σ_i + σ_j. In other words:
  

  σ²   =   (rx_i' - rx_j' )² + (ry_i' - ry_j' )²

  During the time prior to the collision, the particles move in straight-line trajectories. Thus,

  rx_i' = rx_i + Δt vx_i,   ry_i' = ry_i + Δt vy_i
  rx_j' = rx_j + Δt vx_j,   ry_j' = ry_j + Δt vy_j

  Substituting these four equations into the previous one, solving the resulting quadratic equation for Δt, selecting the physically relevant root, and simplifying, we obtain an expression for Δt in terms of the known positions, velocities, and radii.

  

  

  If either Δv ⋅ Δr ≥ 0 or d < 0, then the quadratic equation has no solution for Δt > 0; otherwise we are guaranteed that Δt ≥ 0.

Collision resolution. In this section we present the physical formulas that specify the behavior of a particle after an elastic collision with a reflecting boundary or with another particle. For simplicity, we ignore multi-particle collisions. There are three equations governing the elastic collision between a pair of hard discs: (i) conservation of linear momentum, (ii) conservation of kinetic energy, and (iii) upon collision, the normal force acts perpendicular to the surface at the collision point. Physics-ly inclined students are encouraged to derive the equations from first principles; the rest of you may keep reading.

• Between particle and wall. If a particle with velocity (vx, vy) collides with a wall perpendicular to x-axis, then the new velocity is (-vx, vy); if it collides with a wall perpendicular to the y-axis, then the new velocity is (vx, -vy).
• Between two particles. When two hard discs collide, the normal force acts along the line connecting their centers (assuming no friction or spin). The impulse (Jx, Jy) due to the normal force in the x and y directions of a perfectly elastic collision at the moment of contact is:
  

  

  and where m_i and m_j are the masses of particles i and j, and σ, Δx, Δy and Δ v ⋅ Δr are defined as above. Once we know the impulse, we can apply Newton's second law (in momentum form) to compute the velocities immediately after the collision.

  vx_i' = vx_i + Jx / m_i,    vx_j' = vx_j - Jx / m_j
  vy_i' = vy_i + Jy / m_i,    vy_j' = vy_j - Jy / m_j

Data format. Use the following data format. The first line contains the number of particles N. Each of the remaining N lines consists of 6 real numbers (position, velocity, mass, and radius) followed by three integers (red, green, and blue values of color). You may assume that all of the position coordinates are between 0 and 1, and the color values are between 0 and 255. Also, you may assume that none of the particles intersect each other or the walls.

N                            
rx ry vx vy mass radius r g b
rx ry vx vy mass radius r g b
rx ry vx vy mass radius r g b
rx ry vx vy mass radius r g b

Your task. Write a program MDSimulation.java that reads in a set of particles from standard input and animates their motion according to the laws of elastic collisions and event-driven simulation.

---

Perhaps everything below is just for the checklist??

Particle collision simulation in Java. There are a number of ways to design a particle collision simulation program using the physics formulas above. We will describe one such approach, but you are free to substitute your own ideas instead. Our approach involves the following data types: MinPQ, Particle, CollisionSystem, and Event.

• Priority queue. A standard priority queue where we can check if the priority queue is empty; insert an arbitrary element; and delete the minimum element.
• Particle data type. Create a data type to represent particles moving in the unit box. Include private instance variables for the position (in the x and y directions), velocity (in the x and y directions), mass, and radius. It should also have the following instance methods
  • public Particle(...): constructor.
  • public double collidesX(): return the duration of time until the invoking particle collides with a vertical wall, assuming it follows a straight-line trajectory. If the particle never collides with a vertical wall, return a negative number (or +infinity = DOUBLE_MAX or NaN???)
  • public double collidesY(): return the duration of time until the invoking particle collides with a horizontal wall, assuming it follows a straight-line trajectory. If the particle never collides with a horizontal wall, return a negative number.
  • public double collides(Particle b): return the duration of time until the invoking particle collides with particle b, assuming both follow straight-line trajectories. If the two particles never collide, return a negative value.
  • public void bounceX(): update the invoking particle to simulate it bouncing off a vertical wall.
  • public void bounceY(): update the invoking particle to simulate it bouncing off a horizontal wall.
  • public void bounce(Particle b): update both particles to simulate them bouncing off each other.
  • public int getCollisionCount(): return the total number of collisions involving this particle.
• Event data type. Create a data type to represent collision events. There are four different types of events: a collision with a vertical wall, a collision with a horizontal wall, a collision between two particles, and a redraw event. This would be a fine opportunity to experiment with OOP and polymorphism. We propose the following simpler (but somewhat less elegant approach).
  • public Event(double t, Particle a, Particle b):   Create a new event representing a collision between particles a and b at time t. If neither a nor b is null, then it represents a pairwise collision between a and b; if both a and b are null, it represents a redraw event; if only b is null, it represents a collision between a and a vertical wall; if only a is null, it represents a collision between b and a horizontal wall.
  • public double getTime():  return the time associated with the event.
  • public Particle getParticle1():  return the first particle, possibly null.
  • public Particle getParticle2():  return the second particle, possibly null.
  • public int compareTo(Object x):  compare the time associated with this event and x. Return a positive number (greater), negative number (less), or zero (equal) accordingly.
  • public boolean wasSuperveningEvent(): return true if the event has been invalidated since creation, and false if the event has been invalidated.

  In order to implement wasSuperveningEvent, the event data type should store the collision counts of two particles at the time the event was created. The event corresponds to a physical collision if the current collision counts of the particles are the same as when the event was created.

• Particle collision system. The main program containing the event-driven simulation. Follow the event-driven simulation loop described above, but also consider collisions with the four walls and redraw events.

Data sets. Some possibilities that we'll supply.

• One particle in motion.
• Two particles in head on collision.
• Two particles, one at rest, colliding at an angle.
• Some good examples for testing and debugging.
• One red particle in motion, N blue particles at rest.
• N particles on a lattice with random initial directions (but same speed) so that the total kinetic energy is consistent with some fixed temperature T, and total linear momentum = 0. Have a different data set for different values of T.
• Diffusion I: assign N very tiny particles of the same size near the center of the container with random velocities.
• Diffusion II: N blue particles on left, N green particles on right assigned velocities so that they are thermalized (e.g., leave partition between them and save positions and velocities after a while). Watch them mix. Calculate average diffusion rate?
• N big particles so there isn't much room to move without hitting something.
• Einstein's explanation of Brownian motion: 1 big red particle in the center, N smaller blue particles.

Things you could compute.

• Brownian motion. In 1827, the botanist Robert Brown observed the motion of wildflower pollen grains immersed in water using a microscope. He observed that the pollen grains were in a random motion, following what would become known as Brownian motion. This phenomenon was discussed, but no convincing explanation was provided until Einstein provided a mathematical one in 1905. Einstein's explanation: the motion of the pollen grain particles was caused by millions of tiny molecules colliding with the larger particles. He gave detailed formulas describing the behavior of tiny particles suspended in a liquid at a given temperature. Einstein's explanation was intended to help justify the existence of atoms and molecules and could be used to estimate the size of the molecules. Einstein's theory of Brownian motion enables engineers to compute the diffusion constants of alloys by observing the motion of individual particles. Here's a flash demo of Einstein's explanation of Brownian motion from here.
• Free path and free time. Free path = distance a particle travels between collisions. Plot histogram. Mean free path = average free path length over all particles. As temperature increases, mean free path increases (holding pressure constant). Compute free time length = time elapsed before a particle collides with another particle or wall.
• Collision frequency. Number of collisions per second.
• Root mean-square velocity. Root mean-square velocity / mean free path = collision frequency. Root mean square velocity = sqrt(3RT/M) where molar gas constant R = 8.314 J / mol K, T = temperature (e.g., 298.15 K), M = molecular mass (e.g., 0.0319998 kg for oxygen).
• Maxwell-Boltzmann distribution. Distribution of velocity of particles in hard sphere model obey Maxwell-Boltzmann distribution (assuming system has thermalized and particles are sufficiently heavy that we can discount quantum-mechanical effects). Distribution shape depends on temperature. Velocity of each component has distribution proportional to exp(-mv_x^2 / 2kT). Magnitude of velocity in d dimensions has distribution proportional to v^(d-1) exp(-mv^2 / 2kT). Used in statistical mechanics because it is unwieldy to simulate on the order of 10^23 particles. Reason: velocity in x, y, and z directions are normal (if all particles have same mass and radius). In 2d, Rayleigh instead of Maxwell-Boltzmann.
• Pressure. Main thermodynamic property of interest = mean pressure. Pressure = force per unit area exerted against container by colliding molecules. In 2d, pressure = average force per unit length on the wall of the container.
• Temperature. Plot temperature over time (should be constant) = 1/N sum(mv^2) / (d k), where d = dimension = 2, k = Boltzmann's constant.
• Diffusion. Molecules travel very quickly (faster than a speeding jet) but diffuse slowly because they collide with other molecules, thereby changing their direction. Two vessels connected by a pipe containing two different types of particles. Measure fraction of particles of each type in each vessel as a function of time.
• Time reversibility. Change all velocities and run system backwards. Neglecting roundoff error, the system will return to its original state!
• Maxwell's demon. Maxwell's demon is a thought experiment conjured up by James Clerk Maxwell in 1871 to contradict the second law of thermodynamics. Vertical wall in middle with molecule size trap door, N particles on left half and N on right half, particle can only go through trap door if demon lets it through. Demon lets through faster than average particles from left to right, and slower than average particles from right to left. Can use redistribution of energy to run a heat engine by allowing heat to flow from left to right. (Doesn't violate law because demon must interact with the physical world to observe the molecules. Demon must store information about which side of the trap door the molecule is on. Demon eventually runs out of storage space and must begin erasing previous accumulated information to make room for new information. This erasing of information increases the entropy, requiring kT ln 2 units of work.)

Cell method. Useful optimization: divide region into rectangular cells. Ensure that particles can only collide with particles in one of 9 adjacent cells in any time quantum. Reduces number of binary collision events that must be calculated. Downside: must monitor particles as they move from cell to cell.

Extra credit. Handle multi-particle collisions. Such collisions are important when simulating the break in a game of billiards.

This assignment was developed by Ben Tsou and Kevin Wayne.
Copyright © 2004.

Last modified on January 31, 2009.

Copyright © 2000–2019 Robert Sedgewick and Kevin Wayne. All rights reserved.

https://introcs.cs.princeton.edu/java/98simulation/

This section under major construction.

In 1953 Enrico Fermi, John Pasta, and Stanslaw Ulam created the first "computer experiment" to study a vibrarting atomic lattice. Nonlinear system couldn't be analyzed by classical mathematics.

Simulation = analytic method that imitates a physical system. Monte Carlo simulation = use randomly generated values for uncertain variables. Named after famous casino in Monaco. At essentially each step in the evolution of the calculation, Repeat several times to generate range of possible scenarios, and average results. Widely applicable brute force solution. Computationally intensive, so use when other techniques fail. Typically, accuracy is proportional to square root of number of repetitions. Such techniques are widely applied in various domains including: designing nuclear reactors, predicting the evolution of stars, forecasting the stock market, etc.

Generating random numbers. The math library function Math.random generate a pseudo-random number greater than or equal to 0.0 and less than 1.0. If you want to generate random integers or booleans, the best way is to use the library Random. Program RandomDemo.java illustrates how to use it.

Random random = new Random(); boolean a = random.nextBoolean(); // true or false int b = random.nextInt(); // between -2^31 and 2^31 - 1 int c = random.nextInt(100); // between 0 and 99 double d = random.nextDouble(); // between 0.0 and 1.0 double e = random.nextGaussian(); // Gaussian with mean 0 and stddev = 1

Note that you should only create a new Random object once per program. You will not get more "random" results by creating more than one. For debugging, you may wish to produce the same sequence of pseudo-random number each time your program executes. To do this, invoke the constructor with a long argument.

Random random = new Random(1234567L);

The pseudo-random number generator will use 1234567 as the seed. Use SecureRandom for cryptographically secure pseudo-random numbers, e.g., for cryptography or slot machines.

Linear congruential random number generator. With integer types we must be cognizant of overflow. Consider a * b (mod m) as an example (either in context of a^b (mod m) or linear congruential random number generator: Given constants a, c, m, and a seed x[0], iterate: x = (a * x + c) mod m. Park and Miller suggest a = 16807, m = 2147483647, c = 0 for 32-bit signed integers. To avoid overflow, use Schrage's method.

Precompute: q = m / a, r = m % a Iterate: x = a * (x - x/ q) * q) - r * (x / q)

Exercise: compute cycle length.

Library of probability functions. OR-Objects contains many classic probability distributions and random number generators, including Normal, F, Chi Square, Gamma, Binomial, Poisson. You can download the jar file here. Program ProbDemo.java illustrates how to use it. It generate one random value from the gamma distribution and 5 from the binomial distribution. Note that the method is called getRandomScaler and not getRandomScalar.

GammaDistribution x = new GammaDistribution(2, 3); System.out.println(x.getRandomScaler()); BinomialDistribution y = new BinomialDistribution(0.1, 100); System.out.println(y.getRandomVector(5));

Queuing models. M/M/1, etc. A manufacturing facility has M identical machines. Each machine fails after a time that is exponentially distributed with mean 1 / μ. A single repair person is responsible for maintaining all the machines, and the time to fix a machine is exponentially distributed with mean 1 / λ. Simulate the fraction of time in which no machines are operational.

Diffusion-limited aggregation.

Diffuse = undergo random walk. The physical process diffusion-limited aggregation (DLA) models the formation of an aggregate on a surface, including lichen growth, the generation of polymers out of solutions, carbon deposits on the walls of a cylinder of a Diesel engine, path of electric discharge, and urban settlement.

The modeled aggregate forms when particles are released one at a time into a volume of space and, influenced by random thermal motion, they diffuse throughout the volume. There is a finite probability that the short-range attraction between particles will influence the motion. Two particles which come into contact with each other will stick together and form a larger unit. The probability of sticking increases as clusters of occupied sites form in the aggregate, stimulating further growth. Simulate this process in 2D using Monte Carlo methods: Create a 2D grid and introduce particles to the lattice through a launching zone one at a time. After a particle is launched, it wanders throughout with a random walk until it either sticks to the aggregate or wanders off the lattice into the kill zone. If a wandering particle enters an empty site next to an occupied site, then the particle's current location automatically becomes part of the aggregate. Otherwise, the random walk continues. Repeat this process until the aggregate contains some pre-determined number of particles. Reference: Wong, Samuel, Computational Methods in Physics and Engineering, 1992.

Program DLA.java simulates the growth of a DLA with the following properties. It uses the helper data type Picture.java. Set the initial aggregate to be the bottom row of the N-by-N lattice. Launch the particles from a random cell in top row. Assume that the particle goes up with probability 0.15, down with probability 0.35, and left or right with probability 1/4 each. Continue until the particles stick to a neighboring cell (above, below, left, right, or one of the four diagonals) or leaves the N-by-N lattice. The preferred downward direction is analogous to the effect of a temperature gradient on Brownian motion, or like how when a crystal is formed, the bottom of the aggregate is cooled more than the top; or like the influence of a gravitational force. For effect, we color the particles in the order they are released according to the rainbow from red to violet. Below are three simulations with N = 176; here is an image with N = 600.

Brownian motion. Brownian motion is a random process used to model a wide variety of physical phenomenon including the dispersion of ink flowing in water, and the behavior of atomic particles predicted by quantum physics. (more applications). Fundamental random process in the universe. It is the limit of a discrete random walk and the stochastic analog of the Gaussian distribution. It is now widely used in computational finance, economics, queuing theory, engineering, robotics, medical imaging, biology, and flexible manufacturing systems. First studied by a Scottish botanist Robert Brown in 1828 and analyzed mathematically by Albert Einstein in 1905. Jean-Baptiste Perrin performed experiments to confirm Einstein's predictions and won a Nobel Prize for his work. An applet to illustrate physical process that may govern cause of Brownian motion.

Simulating a Brownian motion. Since Brownian motion is a continuous and stochastic process, we can only hope to plot one path on a finite interval, sampled at a finite number of points. We can interpolate linearly between these points (i.e., connect the dots). For simplicitly, we'll assume the interval is from 0 to 1 and the sample points t₀, t₁, ..., t_N are equally spaced in this interval. To simulate a standard Brownian motion, repeatedly generate independent Gaussian random variables with mean 0 and standard deviation sqrt(1/N). The value of the Brownian motion at time i is the sum of the first i increments.

Geometric Brownian motion. A variant of Brownian motion is widely used to model stock prices, and the Nobel-prize winning Black-Scholes model is centered on this stochastic process. A geometric Brownian motion with drift μ and volatility σ is a stochastic process that can model the price of a stock. The parameter μ models the percentage drift. If μ = 0.10, then we expect the stock to increase by 10% each year. The parameter σ models the percentage volatility. If σ = 0.20, then the standard deviation of the stock price over one year is roughly 20% of the current stock price. To simulate a geometric Brownian motion from time t = 0 to t = T, we follow the same procedure for standard Brownian motion, but multiply the increments, instead of adding them, and incorporate the drift and volatility parameters. Specifically, we multiply the current price by by (1 + μΔt + σsqrt(Δt)Z), where Z is a standard Gaussian and Δt = T/N Start with X(0) = 100, σ = 0.04.

construction of BM.

Black-Scholes formula. Move to here?

Ising model. The motions of electrons around a nucleus produce a magnetic field associated with the atom. These atomic magnets act much like conventional magnets. Typically, the magnets point in random directions, and all of the forces cancel out leaving no overall magnetic field in a macroscopic clump of matter. However, in some materials (e.g., iron), the magnets can line up producing a measurable magnetic field. A major achievement of 19th century physics was to describe and understand the equations governing atomic magnets. The probability that state S occurs is given by the Boltzmann probability density function P(S) = e^-E(S)/kT / Z, where Z is the normalizing constant (partition function) sum e^-E(A)/kT over all states A, k is Boltzmann's constant, T is the absolute temperature (in degrees Kelvin), and E(S) is the energy of the system in state S.

Ising model proposed to describe magnetism in crystalline materials. Also models other naturally occurring phenomena including: freezing and evaporation of liquids, protein folding, and behavior of glassy substances.

Ising model. The Boltzmann probability function is an elegant model of magnetism. However, it is not practical to apply it for calculating the magnetic properties of a real iron magnet because any macroscopic chunk of iron contains an enormous number atoms and they interact in complicated ways. The Ising model is a simplified model for magnets that captures many of their important properties, including phase transitions at a critical temperature. (Above this temperature, no macroscopic magnetism, below it, systems exhibits magnetism. For example, iron loses its magnetization around 770 degrees Celsius. Remarkable thing is that transition is sudden.) reference

First introduced by Lenz and Ising in the 1920s. In the Ising model, the iron magnet is divided into an N-by-N grid of cells. (Vertex = atom in crystal, edge = bond between adjacent atoms.) Each cell contains an abstract entity known as spin. The spin s_i of cell i is in one of two states: pointing up (+1) or pointing down (-1). The interactions between cells is limited to nearest neighbors. The total magnetism of the system M = sum of s_i. The total energy of the system E = sum of - J s_i s_j, where the sum is taken over all nearest neighbors i and j. The constant J measures the strength of the spin-spin interactions (in units of energy, say ergs). [The model can be extended to allow interaction with an external magnetic field, in which case we add the term -B sum of s_k over all sites k.] If J > 0, the energy is minimized when the spins are aligned (both +1 or both -1) - this models ferromagnetism. if J < 0, the energy is minimized when the spins are oppositely aligned - this models antiferromagnetism.

Given this model, a classic problem in statistical mechanics is to compute the expected magenetism. A state is the specification of the spin for each of the N^2 lattice cells. The expected magnetism of the system E[M] = sum of M(S) P(S) over all states S, where M(S) is the magnetism of state S, and P(S) is the probability of state S occurring according to the Boltzmann probability function. Unfortunately, this equation is not amenable to a direct computational because the number of states S is 2^N*N for an N-by-N lattice. Straightforward Monte Carlo integration won't work because random points will not contribute much to sum. Need selective sampling, ideally sample points proportional to e^-E/kT. (In 1925, Ising solved the problem in one dimension - no phase transition. In a 1944 tour de force, Onsager solved the 2D Ising problem exactly. His solution showed that it has a phase transition. Not likely to be solved in 3D - see intractability section.)

Metropolis algorithm. Widespread usage of Monte Carlo methods began with Metropolis algorithm for calculation of rigid-sphere system. Published in 1953 after dinner conversation between Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller. Widely used to study equilibrium properties of a system of atoms. Sample using Markov chain using Metropolis' rule: transition from A to B with probability 1 if Δ E <= 0, and with probability e^-ΔE/kT if Δ E > 0. When applied to the Ising model, this Markov chain is ergodic (similar to Google PageRank requirement) so the theory underlying the Metropolis algorithm applies. Converges to stationary distribution.

Program Cell.java, State.java, and Metropolis.java implements the Metropolis algorithm for a 2D lattice. Ising.java is a procedural programming version. "Doing physics by tossing dice." Simulate complicated physical system by a sequence of simple random steps.

Measuring physical quantities. Measure magnetism, energy, specific heat when system has thermalized (the system has reached a thermal equilibrium with its surrounding environment at a common temperature T). Compute the average energy and the average magenetization over time. Also interesting to compute the variance of the energy or specific heat <c> = <E²> - <E>², and the variance of the magnetization or susceptibility <Χ> = <M²> - <M>². Determining when system has thermalized is a challenging problem - in practice, many scientists use ad hoc methods.

Phase transition. Phase transition occurs when temperature T_c is 2 / ln(1 + sqrt(2)) = 2.26918). T_c is known as the Curie temperature. Plot magnetization M (average of all spins) vs. temperature (kT = 1 to 4). Discontinuity of slope is signature of second order phase transition. Slope approaches infinity. Plot energy (average of all spin-spin interactions) vs. temperature (kT = 1 to 4). Smooth curve through phase transition. Compare against exact solution. Critical temperature for which algorithm dramatically slows down. Below are the 5000th sample trajectory for J/kT = 0.4 (hot / disorder) and 0.47 (cold / order). The system becomes magnetic as temperature decreases; moreover, as temperature decreases the probability that neighboring sites have the same spin increasing (more clumping).

Experiments.

• Start will above critical temperature. State converges to nearly uniform regardless of initial state (all up, all down, random) and fluctuates rapidly. Zero magnetization.
• Start well below critical temperature. Start all spins with equal value (all up or all down). A few small clusters of opposite spin form.
• Start well below critical temperature. Start with random spins. Large clusters of each spin form; eventually simulation makes up its mind. Equally likely to have large clusters in up or down spin.
• Start close to critical temperature. Large clusters form, but fluctuate very slowly.

Exact solution for Ising model known for 1D and 2D; NP-hard for 3d and nonplanar graphs.

Models phase changes in binary alloys and spin glasses. Also models neural networks, flocking birds, and beating heart cells. Over 10,000+ papers published using the Ising model.

Q + A

Exercises

1. Print a random word. Read in a list of words (of unknown length) from standard input, and print out one of the N words uniformly at random. Do not store the word list. Instead, use Knuth's method: when reading in the ith word, select it with probability 1/i to be the new champion. Print out the word that survives after reading in all of the data.
2. Random subset of a linked list. Given an array of N elements and an integer k ≤ N, construct a new array containing a random subset of k elements. Hint: traverse the array, either accepting each element with probability a/b, where a is the number of elements left to select, and b is the number of elements remaining.

Creative Exercises

1. Random number generation. Can the following for computing a pseudo-random integer between 0 and N-1 fail? Math.random is guaranteed to return a floating point number greater than or equal to 0.0 and strictly less than 1.0.
   

   double x = Math.random(); int r = (int) (x * N);

   That is, can you find a real number x and an integer N for which r equals N?

   Solution: No, it can't happen in IEEE floating point arithmetic. The roundoff error will not cause the result to be N, even if Math.random returns 0.9999999999.... However, this method does not produce integers uniformly at random because floating point numbers are not evenly distributed. Also, it involves casting and multiplying, which are excessive.

2. Random number test. Write a program to plot the outcome of a boolean pseudo-random number generator. For simplicity, use (Math.random() < 0.5) and plot in a 128-by-128 grid like the following pseudorandom applet. Perhaps use LFSR or Random.nextLong() % 2.
3. Sampling from a discrete probability distribution. Suppose that there are N events and event i occurs with probability p_i, where p₀ + p₁ + ... + p_N-1 = 1. Write a program Sample.java that prints out 1,000 sample events according to the probability distribution. Hint: choose a random number r between 0 and 1 and iterate from i = 0 to N-1 until p₀ + p₁ + ... + p_i > r. (Be careful about floating point precision.)
4. Sampling from a discrete probability distribution. Improve the algorithm from the previous problem so that it takes time proportional to log N to generate a new sample. Hint: binary search on the cumulative sums. Note: see this paper for a very clever alternative that generates random samples in a constant amount of time. Discrete.java is a Java version of Warren D. Smith's WDSsampler.c program.
5. Sampling from a discrete probability distribution. Repeat the previous question but make it dynamic. That is, after each sample, the probabilities of some events might change, or there may be new events. Used in n-fold way algorithm, which is method of choice for kinetic Monte Carlo methods where one wants to simulate the kinetic evolution process. Typical application: simulating gas reacting with surface of a substrate where chemical reaction occur at different rates. Hint: use a binary search tree.
6. Zipf distribution. Use the result of the previous exercise(s) to sample from the Zipfian distribution with parameter s and N. The distribution can take on integer values from 1 to N, and takes on value k with probability 1/k^s / sum_(i = 1 to N) 1/i^s. Example: words in Shakespeare's play Hamlet with s approximately equal to 1.
7. Simulating a Markov chain. Write a program MarkovChain.java that simulates a Markov chain. Hint: you will need to sample from a discrete distribution.
8. DLA with non-unity sticking probability Modify DLA.java so that the initial aggregate consists of several randomly spaced cells along the bottom of the lattice. This simulates string-like bacterial growth.
9. DLA with non-unity sticking probability Modify DLA.java to allow a sticking probability less than one. That is, if a particle has a neighbor, then it sticks with probability p < 1.0; otherwise, it moves at random to a neighboring cell which is unoccupied. This results in a gives more clustered structure, simulating higher bond affinity between atoms.
10. Symmetric DLA. Initialize the aggregate to be a single particle in the center of the lattice. Launch particles uniformly from a circle centered at the initial particle. Increase the size of the launch circle as the size of the aggregate increases. Name your program SymmetricDLA.java. This simulates the growth of an aggregate where the particles wander in randomly from infinity. Here are some tricks for speeding up the process.
    

11. Variable sticking probability. A wandering particle which enters an empty site next to an occupied site is assigned a random number, indicating a potential direction in which the particle can move (up, down, left or right). If an occupied site exists on the new site indicated by the random number, then the particle sticks to the aggregate by occupying its current lattice site. If not, it moves to that site and the random walk continues. This simulates snowflake growth.
12. Random walk solution of Laplace's equation. Numerically solve Laplace's equation to determine the electric potential given the positions of the charges on the boundary. Laplace's equation says that the gradient of the potential is the sum of the second partial derivatives with respect to x and y. See Gould and Tobochnik, 10.2. Your goal is to find the function V(x, y) that satisfies Laplace's equation at specified boundary conditions. Assume the charge-free region is a square and that the potential is 10 along the vertical boundaries and 5 along the horizontal ones. To solve Laplace's equation, divide the square up into an N-by-N grid of points. The potential V(x, y) of cell (x, y) is the average of the potentials at the four neighboring cells. To estimate V(x, y), simulate 1 million random walkers starting at cell (x, y) and continuing until they reach the boundary. An estimate of V(x, y) is the average potential at the 1 million boundary cells reached. Write a program Laplace.java that takes three command line parameters N, x, and y and estimates V(x, y) over an N-by-N grid of cells where the potential at column 0 and N is 10 and the potential at row 0 and N is 5.Remark: although the boundary value problem above can be solved analytically, numerical simulations like the one above are useful when the region has a more complicated shape or needs to be repeated for different boundary conditions.

    simulated annealing

13. Simulating a geometric random variable. If some event occurs with probability p, a geometric random variable with parameter p models the number N of independent trials needed between occurrence of the event. To generate a variable with the geometric distribution, use the following formula
    

    N = ceil(ln U / ln (1 - p))

    where U is a variable with the uniform distribution. Use the Math library methods Math.ceil, Math.log, and Math.random.

14. Simulating an exponential random variable. The exponential distribution is widely used to model the the inter-arrival time between city buses, the time between failure of light bulbs, etc. The probability that an exponential random variable with parameter λ is less than x is F(x) = 1 - e^{λ x} for x >= 0. To generate a random deviate from the distribution, use the inverse function method: output -ln(U) / λ where U is a uniform random number between 0 and 1.
15. Poisson distribution. The Poisson distribution is useful in describing the fluctuations in the number of nuclei that decay in any particular small time interval.
    

    public static int poisson(double c) { double t = 0.0; for (int x = 0; true; x++) { t = t - Math.log(Math.random()) / c; // sum exponential deviates if (t > 1.0) return x; } }

16. Simulating a Pareto random variable. The Pareto distribution is often used to model insurance claims damages, financial option holding times, and Internet traffic activity. The probability that a Pareto random variable with parameter a is less than x is F(x) = 1 - (1 + x)^-a for x >= 0. To generate a random deviate from the distribution, use the inverse function method: output (1-U)^-1/a - 1, where U is a uniform random number between 0 and 1.
17. Simulating a Cauchy random variable. The density function of a Cauchy random variable is f(x) = 1/(Π(1 + x²)). The probability that a Cauchy random variable is less than x is F(x) = 1/Π (Π/2 + arctan(x)). To generate a random deviate from the distribution, use the inverse function method: output tan(Π(U - 1/2)), where U is a uniform random number between 0 and 1.
18. Generate random point inside unit disc. Incorrect to choose set r uniformly between 0 and 1, θ uniformly between 0.0 and 2π, and use (x, y) = (r cosθ, r sinθ). If you do this, more points close to center of disc. Instead, set (x, y) = (√r cos&theta, √r sinθ) Alternatively, generate x and y uniformly between -1 and 1 and accept if x² + y² ≤ 1. Plot a random sequence of points using both methods and see the bias.
19. Flipping bits. As part of a genetic algorithm, suppose you need to flip N bits independently, each with probability p, where p is some very small constant.
    • Method 1: loop through N bits, generate a Bernouilli(p) random variable for each one and flip accordingly. Takes time proportional to N.
    • Method 2: generate a Geometric(p) random variable X_0 and flip bit X_0; genereate another Geometric(p) random variable an flip bit X_0 + X_1, and so on. Takes time proportional to Np.
    • Method 3: the number of bits to flip in Binomial(N, p). Determine how many bits to flip by approximating with a Gaussian(Np, sigma) random variable. Then flip Z bits, taking care not to avoid duplicates. Takes time proportional to Np, but less calls to transcendental functions.
20. Random point inside N-dimensional sphere. Write a program InsideSphere.java that takes a command line parameter N and computes a random point inside an N-dimensional sphere with radius 1. Generate N uniform random variables deviates x₁, ..., x_N and use this point if
    

    (x1)2 + ... + (xN)2 ≤ 1

    Otherwise repeat.

21. Random point on surface of an N-dimensional sphere. Write a program Sphere.java that takes a command line parameter N and computes a random point on the surface of an N-dimensional sphere with radius 1 using Brown's method. Brown's method is to compute N independent standard normal deviates x₁, x_N. Then
    

    ( x1/r, x2/r, ..., xN/r ), where r = sqrt((x1)2 + ... + (xN)2)

    has the desired distribution. Use Exercise xyz from Section 3 to compute standard normal deviates.

22. Potts model. The Potts model is a variant of the Ising model where each site has q possible directions. (q = 2 corresponds to Ising) The total energy of the system E = sum of - J sigma(s_i, s_j) over all neighbors. The Kronecker delta function δ(x, y) = 1 if x = y and 0 otherwise.
23. 2D Brownian motion. Simulate diffusion of particles in a fluid. Write a data type BrownianParticle.java that represents a particle undergoing a Brownian motion in two dimensions. To do this, simulate two indepedent Brownian motions X(t) and Y(t), and plot (X(t), Y(t)). Create a client program that takes a command line integer N, creates N particles at the origin, and simulates a Brownian motion for the N particles.
24. Brownian bridge. A Browian bridge is a constratined Brownian motion, which is required to begin at the origin at time 0, and end at the origin at time T. If X(t) is a Brownian motion then Z(t) = X(t) - (t/T)X(T) is such a process. To plot, store the intermediate values X(t) and plot after you've computed X(T).
25. Rainbow. In 1637 Rene Descartes discovered the first scientific explanation for the formation of rainbows. His method involved tracing the internal reflections when a light ray is sent through a a spherical raindrop. Simulate the generation of a rainbow according to model of large number of parallel rays hitting a spherical raindrop. When a light ray hit a raindrop, the ray is reflected and refracted. We use the HSB color format, and choose the hue h at random between 0 (red) and 1 (violet). We use 1.33 + 0.06 * h for the refraction index of hue h. For each ray, we plot a single point of light, according to physical laws of refraction and reflection. Each point of light is then plotted in a random color that the observer will see, either in the primary or secondary rainbow. To perform the simulation, we choose one of the 7 colors uniformly at random. Then, we choose a point (x, y) in the unit circle, centered at (0, 0) and set the impact parameter r = sqrt(x² + y²). The angle of incidence θ_i = arcsin(r) and, by Snell's law, the angle of refraction θ_r = arcsin (r / n), where n is the refraction index. If the light ray is totally reflected only once, it emerges at an angle of θ_p = 4θ_r - 2θ_i, contributing to the primary rainbow. If the light ray is totally reflected a second time, it emerges at an angle of θ_p = 6θ_r - 2θ_i - π, contributing to the secondary rainbow. The intensities I_p and I_s of the primary and secondary rays are calculated according to the following transmission and reflection formulas for electromagnetic waves across the boundary of two media.
    

    Ip = 1/2 (s(1-s)2 + p(1-p)2) Is = 1/2 (s2(1-s)2 + p2(1-p)2) p = (sin(θi-θr)/sin(θi+θr))2 r = (tan(θi-θr)/tan(θi+θr))2

    The color intensities I_p and I_s are used to determine the saturation in the HSB color format. Program Rainbow.java simulates this process.

    Rainbow site.

    Last modified on January 31, 2009.

    Copyright © 2000–2019 Robert Sedgewick and Kevin Wayne. All rights reserved.