Physics 305 Lab 6:  Monte Carlo Integration

1. Introduction

In this exercise we investigate multi-dimensional Monte Carlo integration.  A classical example of the Monte Carlo integration approach is to determine the area of a circular "pond" by throwing stones into it.

If you know:
  1.  the total area Atot over which you are throwing the stones, and
  2.  that the locations (including the pond and its surrounding area) where they landed are uniformly and randomly distributed
then the area of the pond is determined as the ratio of the number of stones that hit the pond compared to the total number, multiplied by the total area
Apond = (Nhit/Ntot )* Atot

Multi-dimensional Monte Carlo integration implements this by throwing random numbers into a multi-dimensional space and determining when the values of these random numbers fall inside or outside of the boundary of the volume.

For example, to determine the area of a circle (we denote this by V(2) , the generalized "volume" of the 2-sphere), we can express it as the integral of an area element such as an annular ring:
mc1.jpg

or we can use a Monte Carlo approach where we use the following algorithm:
  1. generate two uniform random deviates x and y  over an interval [0,1]
  2. calculate the Cartesian distance r to the point they correspond to, r2 = x2 + y2
  3.  test if r falls within the circle determined by the maximum radius of integration R
  4. if r <= R, count this test as a "hit"
  5. if r>R count the test as a "miss"
  6. After N total trials, the area of the circle is approximated by    (Nhit / N) * 4R2
where 4R2 is the area of the box with sides D=2R in length which circumscribes the circle.

Simple example: area of the unit circle

Write your own version of the program (which you will use as a template to modify later), which demonstrates how the above algorithm works on the unit circle. 

Using your favorite plotting program (identifying those points inside versus outside by either their marker, color, or both), produce a display like that
shown in Figure 1.  Include this plot in your report, along with a meaningful description of what is being plotted.

pond1.png
    Figure 1: A plot of the data generated by the circle "volume" calculation program (plot via gnuplot).




Exercise Part I:


Multi-dimensional integration where you might really want Monte Carlo: The Volume of a Hypersphere

Of course the example above is kind of trivial, since we can do the integral analytically. So we will try something more interesting.

There are many applications  in physics where integrations must be performed over many more than even 3-dimensions. For example,
 So in many cases, Monte Carlo methods are the only alternative. If you stay in physics (or even if not--such techniques are now widely used in financial modeling) you will most certainly need to become familiar with Monte Carlo methods.

You are now going to make a copy of your program and generalize in in the following way.

The problem is to calculate the  volume of an N-dimensional sphere:  a hypersphere. We might do this analytically if we knew the surface area element, because we could then use an integral like the one for the circle above, which integrates an annular element (the boundary) of dimension one less than the one we want.

But what is the 3-D annular element for a 4-D sphere??

Since this is not something we know off the top of our head, we can just proceed to use the Monte Carlo method to compare the number of "hits" we get within the sphere, to the total number of trials, and then multiply this fraction by the total volume of a  4-D hypercube:

V(4)hypercube  = D4

which we can easily calculate even if it is hard for us to visualize more than a 3-dimensional cube.

It is of course also  possible to make an analytical estimate using integral calculus, but it is not at all trivial; see ( here ) for an example of how to do it.


Exercise Part II:
  1. Complete your program so that it works for calculating volumes for any hypersphere dimension up to 10, and compares its answer with posted analytical values.
  2. To get a feel for what is going on, calculate the volume of 2 and 3-dimensional spheres and check the accuracy vs. number of trials Nmax. Show a plot of fractional error vs. N for 2D and 3D.
  3. Using enough trials to get a few % accuracy, calculate the volume for the unit N-hypersphere from 2 to 10 dimensions using your program. Plot the value as a function of the number of dimensions and note the behavior. Is it surprising?
  4. Compare the value you get at each dimension (up to D=7) to the analytical values.
  5. Check the precision you get as a function of Nmax and the number of dimensions, again up to 7-D. How many trials would you need to get a precision of 10-5 in 10-dimensions?