MOTIVATION: In atomic and molecular systems there are very many degrees of freedom. For example the atoms in the solid, the electrons in the solid, the values of some quantum field at all points in space time. As a detailed example lets consider a gas of A atoms at some temperature T, with a potential v between the atoms at sites i and j which has a partition function Z, that is a multiple integral (one integral for each of the A atoms) written as
We cannot evaluate this for large A. Lets make an estimate for A=20. Lets give each coordinate only 10 values, that means 10 3A points. Say we can do 107 evaluations per second. Then we need 1053 seconds which is more than 1034 times the age of the universe according to some accepted estimates.
In Monte Carlo integration we do not exhaustively cover all these degrees of freedom, but rather sample points randomly; i.e. we convert the integrals into sums. We start with an explanation in one dimension, although the true power of the method is for high dimensions.
A routine to do this is given in chap8a.f, with input of N and an estimate of the required precision.
Run this program and prepare a table of values of N, result and error. How does the error depend on the seed? Think about how to reduce this error. Please mail your table or its webpage location to Dr Adler at phr76ja@tx.technion.ac.il.
Think about one problem with this type of integration. Imagine a function that has a narrow peak about some x-value when all the random x values selected for the sum fall outside the peak? Recall back to the lecture on Gaussian Quadrature, where we discussed non evenly spaced intervals for the sums and how they helped in many cases. There is a similar scheme for Monte Carlo Integration, known as importance sampling. This type of algorithm has revolutionized the study of much of high-energy physics.
The advantage now is that if we choose a weight function that behaves similarly to f(x) then the integrand will be very smooth. Of course this depends on whether we can find a weight function and make the change of variable.
If the distribution of points in y is uniform then x has a distribution dy/dx=w(x), meaning that the points are concentrated about the most ``important'' values of s where w and hopefully f are large.
Lets look at the integral
again with a weight function w(x)=(1/3)(4-2x). This normalises ok, is +ve definite and decreases monotonically over the range of integration. Since f/w=3/4 at both x=0 and x=1 it is very close to f(x). The new integration variable will be y=(1/3)x(4-x) which inverts to give x=2-(4-3y)1/2 and the program chap8b.f, can be used for this integration. Run it for the same values of N as before and add the new results to your table. Please email the table to phr76ja@tx.technion.ac.il.
Lets now look at general dimension. Obviously everything generalises; we have to choose the random coordinates of each component for the vector points xi independently. An example is the integration which compares the area of a quadrant of the unit circle to that of the unit square, by use of the units step function, theta.
A program to do this can be found on the http server as chap8c.f. Run this and see how close to pi you can get.
Sometimes it is easier to use random numbers with a specified distributions. For example if you had to integrate
This means a subroutine such as
subroutine logdst(x)
integer seed
data seed/98347927/
x=-log(1-urand(seed))
return
end
would do the trick.
The preceding analysis requires us to have x(y) analytically, which is not always the case. For example if we had wanted to use w(x)=(6/5)(1-0.5x2) to evaluate the integral
we would have had to solve a cubic equation to find y(x). However this change of variables and inversion can always be carried out numerically. Lets think about tabulating the values x(j) for which the incomplete integral of w wrt x takes the uniformly spaced values y(j) defined as j/M for j=0,1,....,M spanning [0,1]. Thus y(j) is equal to the integral from 0 to x(j) of dx' w(x') and values of x(j) with j an integer chosen from 0,1,...,M with equal probability will then approximate the required distribution. (Special treatment is still needed at the end points.) Generating the x(j) remains the problem, but by integrating the differential equation dy/dx=w(x) thru a simple discretization this can be done, giving a recurrence relation that is started with x(0)=0.
A better method is Van Neumann rejection which is good in multiple dimensions, although we will take a one-dimensional example to explain it. If we want x values between 0 and 1 with distribution w(x), we select w'(x) a positive function > than w(x) over this region. For example w'(x) could be a constant greater than the maximum value of w anywhere in the region. If we now generate points in the two dimensional space that uniformly fill the area under the curve w'(x), and ``accept'' only those under w(x) then the accepted points will be distributed according to w(x). In practice two random variables are chosen, one (xi) distributed proportionally to w' and the other, t, uniformly on [0,1]. The value xi is accepted if t is less than the ratio w(xi)/w'(xi). If a point is rejected we go on and select a new one. This explanation is a little confusing; see more here.
Back to the Lecture 7 page