PHYCS 6730 Lab Exercise: Monte Carlo Integration

These exercises play with Monte Carlo integration. We take examples from Numerical Recipes (2nd edition) section 7.6. The example files are in ~p6730/exercises/monte. The code is in C, with the random number generator in Fortran, but you should be able to read it.

The problem we study is to compute the volume of a truncated bagel (torus). The torus is described by two radii: a smaller one, r and a larger one, R. We start by drawing a circle of radius R in the x-y plane centered at the origin. We make this circle be the center of a cylinder that has been bent around to join itself. The cylinder has circular cross section with radius r. In albebraic terms, the torus is the set of x,y, and z satisfying 

   z2 + (rho - R)2 < r2
   rho = sqrt(x2+y2)
The truncated bagel is illustrated in Numerical Recipes in C online , courtesy of Cambridge University Press.

Exercise 1. Volume of truncated bagel 2D example

If we only want the volume and the center of mass coordinates in x and y, we can think of the integration problem as a two-dimensional integral over x and y with the integrand equal to the length in z at a given x and y. torus0.c does the 2D integral.

Examine the code and answer these questions:

  1. In terms of the algebraic notation in Numerical Recipes, what is the value of the volume integrand at x,y?
  2. What does the "hit rate" measure?
Run the code for 10,000 trials to estimate the volume and c.m. coordinates. Then run it for 1,000,000 trials. How much has the error decreased? Does the error decrease in proportion to 1/sqrt(N)?

Exercise 2.

The code torus1 does the same calculation, but samples in three dimensions instead with a constant integrand. Run this code with 10,000 and 1,000,000 trials and compare the results with those from Exercise 1. Does the error decrease in proportion to 1/sqrt(N)? Which is more accurate? Explain why in terms of the "hit rate"?

Exercise 3.

In the code torus2 we introduce a z-dependent variable mass density rho = exp(5*z) as in Numerical Recipes. There are two ways to handle this: (1) the naive way, as in torus2, sampling uniformly in x,y, and z and summing the density and (2) with importance sampling, as in torus3. To do that we change variables of integration to x,y,s, where 

  s = exp(5*z)/5
  ds = exp(5z)*dz
  z = log(5s)/5
Note that the original integrand contains exp(5z)*dz, but the new integrand contains just ds with the exponential factor eliminated. So uniform sampling in x, y, and s might be more efficient. This is called importance sampling.

Try both methods and discuss what you find.