Monte Carlo Methods

Comp110 Course Module

Introduction to Monte Carlo methods

Many methods in science and engineering involve using a stochastic or random element.  These techniques are powerful and can be applied to many problems;  we will look at two examples in this module.  The algorithms used here are sometimes called Monte Carlo methods after the casino in that European city-state, and we shall discuss these methods not only as a tool for studying many particle systems, but also as a way to compute integrals.  Central to any stochastic method is generating random numbers;  these numbers are provided by an algorithm that creates a sequence of random numbers over and over again.  In Mathematica, the functions RandomReal and RandomInteger provides this functionality.  Generating true random numbers is actually quite difficult.  The sequences you get from RandomReal[ ] are not in every sense truly random, although they should be random enough for practical purposes.  Thus, we will assume that the numbers that come out from these Random functions are indeed random enough.  As from any result computed based on statistics from a random function, the answers you get will vary a little each time.

As a first simple example, we will write a function that estimates π using a Monte Carlo method.  Our algorithm produces these estimates by throwing n darts at a square with corners at {±1,±1}.  Given the coordinates {x, y} of the dart's point of impact, we then test whether the dart hit inside the unit circle where (x^2)+(y^2)≤1.  If the darts hit randomly (with a uniform distribution) inside the square, the chance that a dart lands inside the circle is the ratio of the area of the circle to the area of the square.  Since these areas are π and 4 respectively, the chance of the dart landing inside the circle is (π/4).  In the following exercises, you will use this method to estimate π.  Although computing the area of a circle this way seems to be a very inefficient way of computing π, this kind of approach is useful for computing high dimensional integrals of complicated functions.

Simulating diffusion

In Monte Carlo simulations, a statistical approach, rather than an exact calculation, is used either because the process is inherently random or because it is too complicated.  Radioactive decay is an example of an inherently random process.  The repeated scattering of a photon as it bounces through a fog is an example where the rules are simple, but the overall process is too difficult to analyze.  In both cases there is an element of randomness (time to decay, distance of travel to the next fog droplet) and a willingness to give up details in favor of a simple average quantity (average time to decay, average attenuation of an initial beam of light).

In this section we will model diffusion as a class of systems that is random or stochastic.  For example, if you dump a bottle of perfume in the center of a room, then the spreading of perfume molecules throughout the room can be modeled as a diffusive process.  If the air is still, then the perfume molecules will slowly spread out to fill the whole room.  If we were to watch one particular perfume molecule as it moves through the air, the molecule would undergo a complicated trajectory, moving for a short period in a straight line between collisions with air or other perfume molecules.  Each collision would change the direction of motion of the perfume molecule.  We could, in principle, construct a theoretical description of the motion of the perfume molecules, before we realize that there are more than (10^23)molecules to keep track of.  It seems unlikely that any computer would be powerful enough to complete a calculation of this size any time soon.  However, what we are really interested in is not the individual motion of all the particles, but rather an understanding of the mixing process:  what we really require is a trajectory of an average perfume molecule.  Since the number of particles is large, statistical fluctuations from the average will be small.  Statistical techniques provide the type of information we are after;  each perfume molecule follows a complicated trajectory as it collides repeatedly with other molecules. This trajectory can be described by a 'random walk', a process in which the particle moves one step at a time according to some simple rules.  For the diffusion model, each step corresponds to the particle's trajectory between collisions.

The random walk is a classic case where Monte Carlo methods are useful.  In a fairly general form of the problem, the random walker is assumed to take steps of random length, in a random direction, at random times, all subject to specified probabilities.  If one thinks of starting many such walkers from the origin at time t==0, they will obviously have reached different points by different routes at some later time.  One might then seek to characterize the average motion of  all the walkers by an appropriate average measure, perhaps their absolute distance from the origin.  This problem is an example where fairly simple rules can lead to complicated behavior.  The random walk and some closely related models in fact lead to good descriptions of atomic diffusion, propagation in scattering media, fluid penetration in porous solids, and other topics of current interest.
Many parameters of simple random walks can be calculated analytically, but we will use numerical simulation because it is more general.  Consider a walker confined to a straight line.  At each tick of the clock the walker takes a step of length 1 to the right or to the left with equal probability.  We simulate this process by keeping track of the position after each tick, assuming the walker starts at the origin.  For each step we generate a random number: either  -1 or 1.  The process is repeated for the desired number of steps and the final position noted.  By repeating this procedure a large number of times, we can generate a good sample of final positions, from which we can compute whatever we want.   For example, the average displacement will be zero, but the average distance from the origin will be found to increase like the square root of the number of steps taken.  As with any statistical quantity, there will be an error which will decrease as we accumulate more trials, and which we can estimate by calculating the standard deviation of our result.

One-dimensional random walks

1D algorithm

To make these ideas more concrete, here are a couple of Mathematica functions that will evaluate the average distance for a one-dimensional random walker.  The following function produces a list of random numbers each of which is either -1 or 1.

rand1D[]:=With[{r=RandomInteger[{1,2}]},Which[r==1,-1,r==2,1]]

To compute the final destination of the walker, we make use of the FoldList command that cumulatively sums the elements of a list.  For example, the  following procedure generates a list of numbers that describes the path of a random walker

walk1D[n_]:=FoldList[Plus,0,Table[rand1D[],{n}]]

A typical output, for two runs after 100 steps, can be viewed using the ListPlot command, which plots the position of the particle versus time.  Note that rerunning the ListPlot command will produce a different plot.

ListPlot[walk1D[100],Joined->True]

1D analysis

One interesting question that is worth analyzing is:  On average, how far is a random walker away from the origin after n steps?  In terms of our perfume analogy, we would like to understand how fast perfume disperses on average as a function of time.

Given a large number of random walks, let ≺(Subscript[x, n])≻ denote the average position of these walkers after n steps.  We can estimate ≺(Subscript[x, n])≻ in Mathematica using the function meanDistance that computes m random walks of length n and takes the Mean of last element in these lists.

meanDistance[n_,m_]:=Mean[Table[Last[walk1D[n]],{m}]]

A plot of ≺(Subscript[x, n])≻ for various values of n (averaged over 50 runs) shows that ≺(Subscript[x, n])≻ seems to be zero for all n.

ListPlot[Table[{n,meanDistance[n,50]},{n,100}],PlotRange->{-100,100}]

This observation makes sense since a walker is as likely to step to the left as to the right.  (Note that increasing the number m of runs moves the plot closer to zero.)  So the quantity ≺(Subscript[x, n])≻ is not so informative, since the position of the walker is treated as a signed distance and on average this signed distance just cancels out.

A more interesting quantity is the value ≺((Subscript[x, n])^2)≻, the average of the square of the distance of the walker after n steps.  By squaring the distance of the walker, we will avoid cancellation and hopefully will gain some insight on how fast on average the walker moves away from the origin.

The Mathematica function meanSquDistance computes ≺((Subscript[x, n])^2)≻ in a way similar to meanDistance.

meanSquDistance[n_,m_]:=Mean[Table[((Last[walk1D[n]])^2),{m}]]

Below is a plot of ≺((Subscript[x, n])^2)≻ for m=50 trials.

ListPlot[Table[{n,meanSquDistance[n,50]},{n,100}],PlotRange->{-100,100}]

Observe that roughly these values grow linearly.  If we increase the number of trials to m=200, the linear behavior becomes sharper.

ListPlot[Table[{n,meanSquDistance[n,200]},{n,100}],PlotRange->{-100,100}]

If you increase the number of runs over which you average, the above result approaches a straight line.  The goal in this notebook is to calculate how a large number of atoms spatially move as a function of time. From the above results one can see that

<((Subscript[x, n])^2)>=D n

where n is the number of steps and the factor D is known as the diffusion constant.  From the above figure, we can see that D is approximately 1.

An obvious generalization is to see if this relationship applies in higher dimensions.  In 2D, a random walk on a grid is often taken as a simple model of particle diffusion in the regular crystal lattice.  The quantity of physical interest is the net change in position after a large number of steps.  The purpose of this exercise is to see if equation (1) holds in two dimensions.