### The Monte-Carlo Method

Learning calculus at school, we soon find out that while differentiation is relatively easy, at least for simple functions, integration is hard. So hard indeed that, in many cases, it is impossible to find a nice function that is the integral (or anti-derivative) of a given one. Thus, given ${f(x)}$ we can usually find ${d f /d x}$, whereas we may not be able to find ${\int f(x)\,d x}$.

Stan Ulam was a remarkable genius. Originally from Poland, he spent his working life in America. He played a key role in the development of the atomic and hydrogen bombs at Los Alamos. He developed a method of evaluating integrals that proved very powerful and useful. Although the integral of a function may be defined in a precise, deterministic way, the method of evaluating it is based on probability, and yields an answer that is “very likely” to be correct.

Integrating by Probability

In many problems we need to find the area under a closed curve. In Fig. 1 (left) we show the function ${f(x) = 1/x}$ with the grey area showing the integral from 1 to 7.5. If we are unaware of the log-function, we cannot integrate analytically (to prove that the area is ${\log 7.5}$). Suppose we consider a rectangle around the curve, as in Fig. 1 (right). The sum of the red and blue areas is obviously the breadth by height, or 6.5. But we want the size of the red area.

Ulam’s idea was this: suppose we pick a point ${p}$ at random within the rectangle. The probability that it falls within the red area is

$\displaystyle \mbox{Prob}(p\mbox{\ in Red Area}) = \frac{\mbox{Red Area}}{\mbox{Total Area}} \,.$

where the total area is the the sum of the red and blue areas. If we do the experiment once, it tells us little. But, instead of choosing just one point at random, we choose a large number of them. Then, the proportion ${P}$ of these points falling below the curve should be approximately equal to the probability of a single randomly chosen point falling below the curve:

$\displaystyle P = \mbox{Proportion below curve} = \frac{\mbox{Red Area}}{\mbox{Total Area}} \,.$

If we do the experiment many times, we will get a value for ${P}$. Then, since we know the area of the rectangle, we can solve the equation for the red area:

$\displaystyle \mbox{Red Area} = P \times \mbox{Total Area}$

How could we do this experiment in practice? Well, one way might be to make a kind of rectangular dart-board with the function ${f(x)}$ drawn on it. We might then throw darts repeatedly at the board, and count the fraction of them that land on points below the function.

In reality, we do not play darts; we have a mathematical description of the function${f(x)}$ and the rectangle. We generate a pair of random numbers ${(x,y)}$ and check if it falls below or above ${f(x)}$. We do this many times, and keep a note of how many times the point falls below ${f(x)}$. Then the proportion ${P}$ of these to the total number gives us an approximate value for the area integral ${\int f(x)\,d x}$.

Why Monte-Carlo?

The method was dubbed the Monte-Carlo Method by Nicholas Metropolis, a colleague of Ulam, because of its association with chance and gambling. There are many simpler ways to evaluate an integral like ${\int f(x)\,d x}$. However, in many problems in physics, multi-dimensional integrals must be computed, sometimes in more than 100 dimensions, and the Monte-Carlo method is very valuable in producing an accurate result. The random points are generated rapidly using a PNRG, or Pseudo-Random Number Generator.

The Monte-Carlo Method is one of the most popular and widely-used methods of obtaining numerical answers to problems in physics, chemistry, engineering and other disciplines.

As easy as Pi

The use of a probabilistic method to compute a deterministic quantity goes back to Georges-Louis Leclerc, Comte de Buffon in 1733. He devised a method of calculating ${\pi}$ by dropping a stick or needle on a lined floor (see earlier post here). We will now look at another way to use the Monte-Carlo approach to estimate ${\pi}$.

Suppose a unit circle is drawn within a square and points are chosen randomly within the square; think of a dart-board in a square frame, with darts thrown from a distance. The area of the square is 4 and the area of the circle is ${\pi}$. So, a proportion ${P=\pi/4}$ of the darts will fall within the circle. Fig. 2 (top left)shows the set-up, Fig. 2 (top right) shows 250 random points within the square. Fig. 2 (bottom left) shows 2500 random points and Fig. 2 (bottom right) shows 25000 random points within the square.

We used Mathematica to generate random number pairs ${(x,y)}$ in the range ${[-1,+1]\times[-1,+1]}$ and count the proportion with ${r^2=x^2+y^2<1}$. This proportion should approach ${\pi/4}$ when the number of points is large.

For N=250 it was 3.1554/4, for N=2500 it was 3.1901/4 and for N=25000 it was 3.1422/4. We see that it is approching the value 3.1416/4 or ${\pi/4}$.