In the following we will see how computers are used to solve complex problems in science. Along the way we will need to learn a little statistics and about the numerical technique called the "Monte Carlo algorithm". Our goal is to solve the "wildfire problem" (also known as percolation): imagine a plain filled with patches of dry grass. A fire has started at one end. It has a certain probability of jumping from one patch of grass to another. The question we seek to answer is: What is the probability that the fire crosses the entire plain? To answer this we first need to take a little detour into random numbers.

If you just want to see the percolation applet, click here.

Our first task is to understand what randomness means on a computer. To aid in this
click this
link now.
A window should pop up which contains an applet. When you click "go", the applet generates (=computes. There is a subtlety here)
the numbers 1-6 randomly (representing a die). The histogram to the right keeps track of the total
number of times each number is generated. This happens in real time so the graph will be modified as
you watch it. Click on **go** now. When the applet
is done you'll see that 13 of the hundred dice throws were 6's. We expect this number to be 100/6 = 16 2/3, but of course our process was random, and there is no reason to get exactly what we expect. If you
click **go** again you'll get a different number of 6's.

Can we get a better estimate of 16 2/3? Yes, the more times we throw the dice (generate random numbers) the more likely it is that we'll up with an estimate of 0.1666. Try entering 1000 into the "total number of throws" box (you may have to hit enter while in the box to get it to register). Then hit "go". The number of 6's/1000 should be much closer to 0.1666 than it was when you only did 100 throws. ( I get 0.152 instead of 0.130 when I do it.)

This raises an interesting question. What if I had run the applet ten times with 100 throws each? How does this compare with running it once with 1000 throws? In the second case we know we get a better estimate of the true probability of throwing a 6. What about the first case? We can get a better estimate there too. Here's how: write down your answer for each of the ten runs and make a histogram of the results. It should look something like this.

If I repeat the experiment 100 times (and I did) I get the following histogram.

**Exercise**
Look up the formula for a normal distribution. See if you can estimate
the average and the standard deviation by comparing the formula to the graph
above. You may want to plot them on the same graph and adjust the strength,
mean, and standard deviation of the normal distribution so that it looks
like the data.

Now you know how random numbers are calculated on a computer and a little
about how to analyse the results of statistical experiments. Let's move
on to computing something.
This example involves calculating the area of a circle. In this case we know
that the answer is pi R^{2} (I have chosen R=1 in this example), but
in general it may be an area that is too complicated to evaluate conveniently
and we want to develop a computational method. Imagine throwing darts at a
square of size 2 by 2 with a circle in it of diameter 2. If all of our darts
land in the square what is the probability that one lands in the circle? The
answer is given by the ratio of the area of the circle to the square.
In this case pi (D/2)^{2} divided by D^{2}. Thus if we keep track of
the number of times a dart has landed in the circle, we can get an estimate
of the area of a circle! To see this in action, close the dice window and
open this window.

The box on the right is the square which we are throwing darts at. When you click **go**, each
dart throw is shown as a red dot if it lands in the circle and blue if it
lands outside the circle. The area of the circle is then calculated as the
number of darts landing (the third box) divided by the total number of darts
times the area of the box (in these units it is 4). This should be pi=3.141592654... My first run with 100 darts gives 3.08, not too bad. With 1000 darts I
get 3.16. You can make a histogram of results just as with the dice example, and
analysing the histogram will allow you to make a better estimate of pi and of
your error in the estimate.

Throwing darts is the basis of the **Monte Carlo method**. In practice it is
wasteful to throw darts into the square so we always try to find a shape
which is close to the one we want to calculate the area of (things get more
complicated in many dimensions.) But remember, we must know the area of the
reference shape! Otherwise we don't know what to multiply by. Choosing a
good reference shape is called **importance sampling** and it is a vital
part of the Monte Carlo method.

**Exercise**
Evaluate the
area of the circle 20 times with 100 throws each. Make a histogram of your
results. Make a rough estimate of the error in your best estimate of the
area.

**Exercise**
Write a computer code which calculates the volume of a sphere.

**Exercise**
Write a Monte Carlo code which evaluates the following integral.

**Exercise**
Evaluate the probability for a dart landing in a D-dimensional
super-sphere of radius 1 if it is embedded in a D-dimensional super-cube
with sides of length 2. You should get (take D to be even).
Plot your answer vs. D. You should see that it
is exceedingly unlikely to hit the super-sphere once D gets large.

We are now ready to tackle the wildfire problem. Our approach will be to make
a simple **stochastic** (random) model of the system. We will then analyse
it using an applet. Along the way we will see many surprises.

Here's the model: divide a square into smaller squares and color each square red with a certain probability, p. For example, we may have a 4x4 grid of squares and color each square red with probability 0.6. We will let the computer do this. In practice this means that the computer generates a random number between zero and one for each square. If the random number is smaller than 'p', the square is marked, otherwise it is left unmarked.

Here is an example of doing this on a 6x6 grid with p=1/2.

Eighteen squares are colored, this is right since we expect about 1/2 of the 36 squares to be marked. It is luck that it is exactly 1/2 in this case! Each square is marked separately so anywhere between 0 and 36 marked squares are possible, with 18 being the most likely. The applet has an algorithm which indicates whether marked squares are next to each other -- in other words it figures which red squares are in connected clumps. For convenience it colors the largest such clump blue. In the figure this is the 9 square clump in the upper right. An example on a grid of 120x120 squares is shown below.

Notice that things get a lot more complicated when there are a lot of squares!

If a connected area goes from top to bottom we say that the system has "percolated". This will be indicated by the applet as a green clump. Here is an example from a 120x120 grid:

Here's what we want to find out: for a given probability of marking each square what is the probability of percolating? Also, how does this change with the grid size? The answer will amaze you!

Before we start, let's sharpen our understanding of the model by looking at some special cases. First, if p=0, then no squares will ever be marked and the probability of percolation will be zero. If p=1 then all the squares will be marked and the probability of percolation will be unity. In between, we expect some intermediate answer.

Now consider a 2x2 grid. It is a fairly simple matter to prove that the
probability of percolation is p^{2}(2-p^{2}).

**Exercise **Prove this

Let's use the computer to check this result. Close the pi window and open
the percolation window now. Select a 2x2 grid using the pull down menu on the
top right (it should already show 2x2). Make sure that the next pull down menu
reads 'percolation' and that the "prob." window reads 0.3. Now click **go**.

What you'll see is a series of 2x2 grids flashing by (two per second). Each one
of these represents the result of one experiment where the squares are marked
with probability,p, in this case 0.3. Each one of these is called a
**configuration**. The wide "total number of configurations" box shows
how many configurations that the computer will generate. The smaller boxes
to the left indicate what is happening. From top to bottom: the number of
configurations which have been generated at that time, the number of
configurations which have successfully percolated, and the ratio (labelled
"stat"). When I ran this I found 3 of 20 configurations had percolated, giving
a ratio of 0.15.

If you now hit **plot** you'll get a little graph. The x-axis is the probability you ran at (0.3 in this case)
and the y-axis is the probability of percolation (0.15 in this case). You
should see a blue dot at the point (0.3,0.15). This is the experiment you
just ran. You should also see a blue curve. This is the analytic result given
above. Notice that your point doesn't lie exactly on the curve -- its
stochastic! As discussed above, the more configurations you generate the
closer you will get to the theoretical curve.

Now let's fill out the whole curve. if you click on **auto** the applet will automatically cycle through various
probabilities and plot the results. If you find the pace too slow, put a smaller
number into the "pause" box (this is the pause between configurations measured
in milliseconds).

**Exercise **Increase the
number of configurations and generate a new curve. It should be more accurate.

**Exercise **Think of a
method to associate an error with each point you calculate.

What happens as the grid size increases? The limit of large size systems is
of interest to physicists because it is a good approximation to real systems
which contain huge numbers of degrees of freedom (atoms, molecules, etc). In
fact we are often interested in the **bulk limit** where the system size
goes to infinity.

Let's get an idea of what happens as the grid size (we'll call it **L**
for an LxL grid) increases by running the **auto** experiment for a 4x4 grid. To do this select 4x4 in the top
pull-down menu. Compare the 2x2 and 4x4 plots. The 4x4 plot should look steeper
in the middle section than the 2x2. Repeat this for larger systems.

**Exercise **What is
happening as L increases? Try to put all of the plots onto a single graph
with each curve labelled by L. Speculate on the infinite L shape of the curve.

What you will find is that the curve gets sharper and sharper as L is increased.
In fact the infinite L curve is a step function: it is zero up to a special
value of p called p_{c} -- the critical probability. Above p_{c}
the percolation probability is 1. This is an example of a phenomenon called a
**phase transition**. For infinite systems there is a magical probability
below which you *never* percolate and above which you *always*
percolate.

Our task in this section is to find an accurate estimate of p_{c}.
If you run **auto** for the 120x120
system, you can get a rough estimate of the critical probability. When I do it,
it looks like p_{c} is between 0.5 and 0.6. To do better than this we need
to run lots of cases around p=0.55 checking to see if the percolation
probability is (near) 0 or 1.

We'd like to look at the percolation probability on an infinite system, but of
course this is impossible. The best we can do is watch what happens to the
probability as we increase L. Imagine fixing p at 0.5 and calculation the percolation probability for L=2,4,6,8,etc. Then plot the probabilities vs. L and try
to see where is is going, **extrapolate** to L=infinity. Better yet, plot
the percolation probability vs. 1/L. Then the L=infinity corresponds to the
intercept at the origin.

The applet does this for you: choose a probability (say 0.3) and select
"extrapolate" from the second pull down menu. Then click **go**. You will get a graph of the perc prob
vs 4/L for L=4,10,20,40,60. For p=0.3 it is pretty clear that perc prob(L=infinity) = 0.
Try it for p=0.7. It is clear that perc prob(L=infinity) = 1.

**Exercise ** See how
accurately you can determine p_{c}. Compare to my answer.

You will find that it is pretty hard to do the extrapolation near p=0.6. This is because you are near the phase transition and the capabilities of stochastic methods break down there. To make progress you may need to use larger lattices (120x120) and run for a lot of configurations.

**Exercise **Make your own
1/L graph including larger lattices with many configurations.

Before we move on, there is another subtlety which we shouldn't ignore. Why
did I plot the perc probs vs. 4/L? Why not 1/L^{2} or exp(-L)? Sometimes
physicists have theoretical reasons to choose one form over another. But often
we don't know and have to guess. The error induced in our conclusions this way
is called **systematic error**. There is little we can do about it except
try to estimate it reliably. The errors which arise because we are doing
numerical experiments with random numbers are called **statistical errors**.

**Exercise **Make some
extrapolation graphs using 1/L^{2} as a plotting variable. Does it do
better or worse than 1/L?

We now wish to analyse the structure of percolating clusters more carefully. To start, let's take a look at one on a 120x120 grid.

Notice that the percolating cluster (in green) is not very smooth. It is filled
with holes and its border wanders all over the place. If you look in the red
"hole" at the center bottom, you will see that it is full of holes too. In fact
if you blow up any region of this plot (and we had run it for a very large
grid size) you'll see a pattern that looks very similar to what you were
looking at before. (You may have seen demonstrations of this with something
called the **Mandelbrot set**.) This is an example of the phenomenon of
**scale invariance** which occurs at phase transitions. Scale invariance
means that the thing you are looking at looks the same at all magnifications.

The jaggedness of the border of the percolating region may remind you of something else you've probably heard of -- fractals. A fractal is a surface which is so bumpy that it is not really a line and not really an area -- it is something in between. A typical example is a coast line which looks jagged and bumpy no matter what scale you draw the map at.

More to come on the definition of the fractal dimension and an applet demonstration of it

The fractal dimension is expected to be 91/48 = 1.90.

We continue our analysis of the character of percolating clusters by examining the fraction of squares which are in the percolating cluster. We call this fraction F. It turns out that F obeys a simple formula near the critical probability:

Here is a number known as a **critical exponent**. It tells us how the fraction of squares in the percolating cluster goes
to zero as we approach p_{c} from above. Its value is 5/36 = 0.139.

The applet provides a way to see this scaling law. Select "F calculator" from the
second pull down menu. Choose a lattice size and then click **auto**. The applet will cycle through various
probabilities and calculate the average number of spins in each percolating
cluster, F. In the bulk limit we expect F to be zero below p_{c}. If
you run on a 40x40 or larger lattice you should confirm this. Above p_{c} it is less clear what will happen, except that F=1 when p=1. What you will
find is that F approaches zero rapidly as p approaches p_{c} from
above. This means that even though a cluster has percolated, the fraction of
spins in it is essentially zero at the critical point! Percolating clusters
must be very wispy things near p_{c}. Of course, this is a reflection
of their fractal nature.

**Exercise **How would you
measure ?

It turns out that critical exponents have a very general validity. Imagine that
we solve a variety of two-dimensional percolation problems, say with a square
lattice (as we have been doing), a triangular lattice, and a honeycomb lattice.
Now calculate for each theory. You'll find
that it doesn't change! This is an example of a phenomenon called
**universality**. What it means in this case is that the
structure of the percolating cluster doesn't depend on the details of the
model that you use to study percolation. This is a very deep property of
physical systems near phase transitions. It underlies a concept called
the **renormalization group** in particle physics and the epsilon-expansion
in condensed matter physics. In 1982 Ken Wilson won the Nobel Prize by
exploiting these concepts to invent a
way to make calculations with phase transitions.

**Exercise **Think about
the percolation problem in two dimensions with a triangular lattice. Do
you expect p_{c} to be
larger or smaller than it is for the
square lattice? Write a code to simulate it. What fractal dimension does
its percolating clusters have?

**Exercise **Examine percolation for a three dimensional square lattice (a grid of cubes). What happens
to p_{c}
now?

**Exercise **Examine universality by calculating different quantities in several different two-dimensional
models. Why isn't p_{c} universal? Can you think of any other
quantities which are universal?

This site is now a part of the Merlot Mathematics Digital Library

*Introduction to Computational Physics*, Giardino.*Introduction to Percolation Theory*, D. Stauffer.*Numerical Recipes in C*, W.H. Press et al.