 
   
\documentclass[11pt]{article}
\usepackage{html}
\pagestyle{empty}
   
% Fit to page specifications:
%
%       Defaults for article:   textwidth = 5 in
%                               textheight = 7.25 in
%                               hoffset = 0 in
%                               voffset = 0 in
%
%       These settings center the text on the page, with right and
%       left margins of 1.75 in and top and bottom margins of 2 in
%       for standard U.S. letter size paper.
%       Changing hoffset moves the text to the *right* by hoffset.
%       Changing voffset moves the text *down* by voffset.
   
\textwidth 6.5in
\textheight 9.5in
\voffset -0.3in
\hoffset -0.85in
   
\begin{document}
   
\begin{rawhtml}
        <center><img src=../../gif_files/GRAD_BLUE_LINE_f.gif></center>
\end{rawhtml}
  
\begin{center}
        {\LARGE\bf Random Numbers in Physics }
\end{center}
  
\begin{rawhtml}
        <center><img src=../../gif_files/GRAD_BLUE_LINE_f.gif></center>
\end{rawhtml}



This section describes some of the uses of random numbers in physics
applications.  These applications range over a wide variety from
describing phenomena that are intrinsically at random, like Brownian
Motion, to performing integrals in multi-dimensional phase space.

Ref:  {\em Computational Physics}, Koonin and Meredith, 1990



\subsection{ Generating Random Variables }

The first task in these applications is to generate random numbers.
These random numbers could be uniformly distributed or have
arbitrarily probability distributions.  It is also of importance to
understand how to test those random numbers.

It is not quite obvious how to generate random numbers on a
deterministic computer.  What is clear is that computers can repeat tasks
and are consistent in doing so.  Therefore the generation of uniformly
distributed random numbers has to follow specific algorithms which
can be repeated.  Which means that, starting from a seed, a random
number algorithm always generate the same random numbers sequence when
run repetitively.  As such random numbers generated on a deterministic
computer are not really random, they are said to be pseudo-random.

A very common algorithm used to generate a sequence of pseudo-random
numbers from a seed number is via a linear congruential method.  The
numbers are generated from the sequence described in the following
equation

$$
x_{i+1} = ( a x_i + c ) mod m
$$
where $x_{i+1}$ is generated from $x_i$ in a recursive manner.  The numbers 
$a$,$c$  and $m$ are chosen  carefully according to the computer architecture 
and
often times are very large.  The book Numerical Recipes devotes a
chapter to this method of generating random numbers.  It describes
these methods as an art rather than a science.

This C functions 
\begin{rawhtml}
<a href=../random_numbers/ran2.c>ran2.c</a>
\end{rawhtml}
and 
\begin{rawhtml}
<a href=../random_numbers/ran3.c>ran3.c</a>
\end{rawhtml}
 provided by Numerical Recipes
follow such methods based on linear congruential method.  In addition
to following that the generation of a sequence through the linear
congruential equation, the methods shuffle the orders of the numbers
generated by using its own generated random numbers.  An
explanation is provided in Numerical Recipes.

The {\it gcc} compiler provides a function, {\em rand()} which generates
integer pseudo-random numbers between $0$ and  {RAND\_MAX}. The man 
page ( {\em man rand()} ) says:

{\it The  versions  of rand() and srand() in the Linux C Library use the
       same random number generator as  random()  and  srandom(),  so  the
       lower-order  bits  should  be  as  random as the higher-order bits.
       However, on older rand() implementations, and on current  implementations
       on  different  systems, the lower-order bits are much less
       random than the higher-order bits.  Do not  use  this  function  in
       applications  intended  to  be  portable  when  good  randomness is
       needed }.

For most applications the use of pseudo-random numbers is perfectly
adequate.  What is important is that these numbers be statistically random.
These numbers should be uniformly distributed  and show no
correlations in any form whatsoever among themselves.  

 A very
elementary test that those numbers must past to be acceptable is
uniform distribution.  This is best illustrated by using the random
numbers generated to specify points in a plane.  The
following program 
\begin{rawhtml}
<a href=../random_numbers/random_points.c>random_points.c</a>
\end{rawhtml}
generates a sequence of x, y, values using the C
language {\em rand()}function to generate random numbers scaled 
in the interval [0,1].  

{\em Exercise: } Plot the
sequence of points generated by the program by using GNUPLOT
and watch for
uniform coverage.

{\em Exercise: } 
Do the same for rand2 and rand3 from Numerical Recipes.  Watch out for
patterns in the coverage. There should be none of course.



\subsection{ Exponential Random Numbers}

The basic rule relating two continuous probability
distributions $f(y)$ and $w(x)$ is
$$
   f(y) dy = w(x) dx
$$

Let us assume a domain [0,1] in $y$ and $x$ and that the
$y$ variable is uniformly distributed, $f(y) = 1$. Then
$$
   y(x) = \int_0^x dx' w(x')
$$
with the condition
$$
  y(1) = 1 = \int_0^1 dx' w(x') 
$$

To get an exponential distribution set $ w(x) = exp(-x) $ and extend
the domain to $\infty$ to get the transformation
$$
   y(x) = 1 - exp( -x    )
$$
or
$$
  x = - ln( 1 - y )
$$



\subsection{ Gaussian Random Numbers}

A particularly interesting distribution of random numbers is that of
Gaussian probability.  This is a distribution of random numbers which
are generated according to a Gaussian density of points. The latter, with
mean $0$ and variance $\sigma =1$, is
$$
    w(x) = ( 2 \pi )^{- \frac{1}{2} } exp( - x^2/2 )
$$

A very efficient method for generating normally distributed variables
is to consider a Gaussian distribution in two dimensions ($x_1$, $x_2$) for
which the number of points in a differential area is proportional to
$$
   exp( - \frac{1}{2} ( x_1^2 +x_2^2 ) ) dx_1 dx_2
$$


In terms of the usual polar coordinates
$$
    r = ( x_1^2 + x_2^2 )^\frac{1}{2} , \theta = tan^{-1}\frac{x_2}{x_1}
$$
the distribution is
$$
     exp( - \frac{1}{2} r^2 ) r dr d\theta
$$
or, with a further change of coordinate $u  = \frac{1}{2} r^2$,
$$
   exp( - u ) du d\theta
$$


Hence, if we generate $u$  between $0$ and $\infty$  with an exponential
distribution and $\theta$ uniformly between $0$ and $2 \pi$, then the
corresponding values of $x_1$ and $x_2$
$$
   x_1 = ( 2 u )^\frac{1}{2} \cos{\theta}
$$
$$
   x_2 = ( 2 u )^\frac{1}{2} \sin{\theta}
$$
will be  distributed normally.

To generate an arbitrary Gaussian distribution with mean $M$ and
variance $\sigma$
$$
w(x) = ( 2 \pi \sigma^2 )^{- \frac{1}{2} } exp( - \frac{(x-M)^2}{\sigma^2} )
$$
we only need to take the numbers generated above, multiply them by
$\sigma$, and add to them $M$



The program 
\begin{rawhtml}
<a href=../random_numbers/gaussian_points.c>gaussian_points.c</a>
\end{rawhtml}
 generates such random numbers.  

{\em Exercise: }
The thing to try is to generate ten
thousand numbers and plot a histogram of these numbers.  To do so
write a bin code which accepts the Gaussian numbers from standard
input and bin those numbers in arbitrarily chosen bin values.  The
results of such a filter will be the basis of a histogram distribution
of the numbers.  You should find that the code {\em gaussian.c}
  generates random
numbers centered on $x=0$ and  width $\sigma = 1.0$.



\subsection{ Random Walks }



Brownian Motion describes the movement of small particles in suspension
in liquids.  This motion is observed to be in a zig-zag
path due to the collisions of the molecules in the liquid 
with the small particles, each
collision impacting momentum on the small particles.

Random walks can be generated on the computer using pseudo-random numbers.
  An easy way to generate a random walk is to do so on a
square lattice.  The model is that of a walker who starts at an
arbitrary location on a two dimensional lattice with equal spacing
in both directions, 1 and 2.  At that intersection the random walker
decides to go north, east, south or west at random with an equal
probability for each outcome.  Once decided to go in a specific
direction the random walker takes one step in the given direction.
This one step takes him to the adjacent site.  Then the random walker
takes subsequent steps following the same algorithm.  The program
\begin{rawhtml}
<a href=../random_walk.c>random_walk.c</a>
\end{rawhtml}
implements this procedure.

The basic block of {\em random\_walk.c} is a function to generate the random
decision for the walker step.  The intrinsic C function, {\em rand() },
generates a random integer between 0 and a maximum value {\em MAX\_RND}.  A
simple scaling brings the random into a double C variable between 0
and 1.  Assuming uniformity in the distribution of the numbers a
simple structure of if statements allows to decide if the random
number falls in the interval [0.0, 0.25] [0.25, 0.50], [0.50, 0.75] or
[0.75, 1.00].  Depending on the value of the random number generated
by {\em rand()} a number 0, 1, 2 or 3 is thus returned from the function.
This allows the algorithm to decide which way the random walker will
take his step.

{\em Random\_walk.c} can generate various random walks
as specified by the parameter {\em -w \#} on the command line.  
These walks are based
on the same sequence of random numbers generated by {\em rand()} except for
the fact of using different groups      of numbers from the sequence.  
Therefore the
random walks are all very different.  These are best illustrated via the
graphic techniques described in the next section.

Note that the random walks generated by {\em random\_walk.c} appear to have
different lengths.  That is, the last step takes the random walker at
different distances from the origin of the walk at position [0,0] on
the lattice.  This is perfectly normal given that the walks are
totally at random.

A random walk can also be generated via steps of non-uniform lengths.
For instance, more randomness can be introduced in the walk if the
random walker is to decide the length of the step to take as well as
its direction.  For instance, the random walker at the given location
could take a step based on random numbers with Gaussian distribution to
provide the length of the step and a random choice between the
directions north, east, south or west. The latter can also be
generalized to a random angle between $0$ and $ 2 \pi$.
 This will generate random walks of
very different topology than the random walk on a lattice.



\subsection{ It's all in the Graphics! }


It would be interesting to display random walks one after another 
as they are generated.  In order to do so the {\em random\_walk.c} code must
interact with gnuplot.  The concept to use is therefore a Unix
code-to-code communication pipe.  This requires to look a bit more at
C syntax.

The idea is to establish a communication channel from one program
to another.  In other words, to do an equivalent of the standard
{\em stdout} piping into {\em stdin} between codes, but
 {\it while the codes are executing} and not involving
a {\em stdin} and {\em stdout} pipe!  
The C language and Unix offers two functions to establish and
close such  open pipes, {\em popen()} and {\em pclose()}.
To facilitate the tasks of the programmer the C language uses a syntax
to read/write through a Unix pipe
analogous to writing and reading from a disk file.
This has the great advantage that you can use {\em fprintf()} into a pseudo file
stream to communicate through the pipe.

To establish a pipe requires to establish a file
handler, a file pointer, which will corresponded to the channel for
the pipe, open the pipe, read from or write to the
pipe and close the pipe once done.  

1) {\em FILE *fp}  establishes a file pointer.

2) {\em fp = popen( ``gnuplot'', ``w'' )} starts gnuplot 
(or any other application),
opens a pipe to this application and allows this pipe to be
accessed via the file pointer fp for writing.

3) {\em fprintf( fp, ``command'' ) } sends (writes) information 
(command) to gnuplot.

4) {\em fflush( fp )} empties the buffer associated with the file
pointer fp.  This makes sure that gnuplot receives all the
information in the command that have been transmitted through the pipe.

5) {\em pclose( fp )} closes the communication pipe to gnuplot. Note
that you must make gnuplot quit through the pipe.



These provide the logical steps for starting gnuplot, opening a pipe to it,
transmitting commands to it, and closing the pipe.  Since the calculations
can happen extremely fast and the piping/graphing not so fast, you may want to
slow down the procedure by making the code sleep a few moments in
between the actions taken by gnuplot.  {\em sleep(time)} allows 
to insert a {\it time} 
seconds delay between the commands forwarded to the
pipe.


The code  
\begin{rawhtml}
<a href=../pipes/gnuplot_pipe_demo.c>gnuplot_pipe_demo.c</a>
\end{rawhtml}
illustrates the procedure for very
simple data set.  It picks up a file and plots it. The code
\begin{rawhtml}
<a href=../pipes/moving_circles.c>moving_circles.c</a>
\end{rawhtml}
generates the data for this demo.  Note how the 
frames are separated in the data.

To display random walks in sequence using gnuplot requires a bit
more thinking.  The following steps need to be taken:  
read the data through
a pipe via a redirection from {\em random\_walk.c}.  Store the data into a
temporary file under /tmp.  Open a pipe to gnuplot and instruct it to
plot the data points that were stored in the file.  Receive new data in the
redirected input and replace the file under /tmp with this new data.
Instruct gnuplot to graph this new random walk.  Continue until the
end of the externally piped stream is reached.

Storing the file on the /tmp instead of the user
directory is good practice.  
On all computers /tmp is public scratched disk space usable by
any user on the system.  After you create the data file on /tmp you
should make sure the code removes that data file.

A signal must be transmitted to the code to indicate the end of each
random walk.  This is easy to do with a \# followed by a comment.  \#
followed by comments are ignored by gnuplot in plotting
graphs.  Therefore the file can be plotted by normal means if the
data is stored in this format.  But this constitutes a easy 
marker to recognize in between the frames to produce a
sequence of frames.

The code 
\begin{rawhtml}
<a href=../pipes/gnuplot_pipe.c>gnuplot_pipe.c</a>
\end{rawhtml}
implements the procedure described above.  To use it, type

{\em random\_walk -w 5 | gnuplot\_pipe }

This will create 5 random walks that will be displayed  in
sequence. 

{\em Exercise: }
Explore the various options in {\em gnuplot\_pipe.c}.
In particular, use the options {\em -x xmin xmax} and
{\em -y ymin ymax} to fix the x- and y-ranges for
a movie like display.





\subsection{ Diffusion Limited Aggregates (DLA) }


A very interesting physical phenomena is that of electrolysis.  
In this process a
chemical compound is dissociated via establishing an electrical
potential between different points in the compound.  As an example,
imagine a dish of cylindrical shape with a cathode and an anode given
by a point in the center of the dish and a piece of metal surrounding
the dish on the perimeter. Upon establishing a potential
difference between the electrodes
 the molecules of the compound will disassociate and the two
ion species will gravitate towards the two  electrodes.

An actual experiment based on zinc sulfate dissociation is shown in
\begin{rawhtml}
<a href=http://polymer.bu.edu/doc/ECD.html>experimental archive</a>.
\end{rawhtml}

The shape of the accumulated material on the center electrode is a
fractal.  It results from the ions produced by the separation of the
molecules in the liquid that attach themselves to the
atoms already accumulated at the center of the dish.

A figure of this type can be modeled on a computer.  The algorithm to
do so is based on random walks.  This is quite realistic
since the ions are themselves subjected to random walks in the
physical system.

The algorithm  is quite simple.  It consists
in launching random walks  from a circle at a large distance,
and letting those random walks come in contact with the aggregate near
the center electrode.  Upon getting in contact, in the sense of the
nearest neighbor on the lattice, the ions attach themselves to the
aggregate.  Repeat these steps until the aggregate is large enough.

Here are some details about the algorithm:


\begin{enumerate}

\item Establish a large circle

\item Launch a random walk from the large circle at a random angle
(uniform distribution)

\item At each step of the random walk check whether the random walker is near the
aggregate.  If so, attach the last point of the walk the aggregate on the
lattice.

\item At each step of the random walk check that the random walk 
 has not gone beyond
the large circle.  If so, stop the walk.

\item Repeat the same procedure of launching a new random walk in the large
circle until the maximum number of walks has been reached or the
aggregate has reached the specified size.

\end{enumerate}

This algorithm is summarized in the following pictures:

\begin{rawhtml}
       <img src=../DLA/DLA_nohit.jpg  height=200 width=200 border=0>
\end{rawhtml}
\begin{rawhtml}
       <img src=../DLA/DLA_hit.jpg  height=200 width=200 border=0>
\end{rawhtml}




The program 
\begin{rawhtml}
<a href=../DLA/dla.c>dla.c</a>
( <a href=../DLA/DLA.cpp>DLA.cpp</a> )
\end{rawhtml}
implements this algorithm.  

The program 
contains  options to make  poor man movies of the random walks building-up 
the fractal and of
the evolution of the
aggregate.  Namely the data is printed out in a format compatible
with {\em gnuplot\_pipe.c}  in order to
visualize the random walks that build the DLA fractal, option {\em -p},
or the evolution
of the fractal aggregate, option {\em -m}.

The resulting object is a fractal.  It looks as if it could be a tree
or a plant (maybe of strange shape).  Mathematically it is a
fractal in the sense that the dimension of this object is between 1
and 2.  This could be estimated by the block algorithm which consists
in covering the object with small cells and estimating the number of
cells needed to cover the object as the size of the cells becomes smaller and
smaller.  This procedure 
for a plane will give dimension 2 while for a line it would
give dimension 1.  This procedure for the DLA fractal would yield a
dimension that is neither 1 or 2 but rather a number in between.

Run {\em dla -h} to see the code options.



\subsection{ Monte Carlo Methods }

Systems with a large number of degrees of freedom are often studied in
physics. Many-body systems ranging from many electron atoms to a
chunk of metal are examples of such systems. 
The solution of those systems often requires 
integrals over multi-dimensional spaces, typically {\em 3A} dimensions. 
An example is the partition function for a system at temperature
{\em T} ($k$ is Boltzman constant and $E$ the energy),
$$
   Z = \int d^3r_1 ... d^3r_A \exp( - \frac{1}{k T} E( r_1, ..., r_A )  )
$$
We will use this function later to solve the Ising Model.

The difficulty with calculating  numerically integrals in multi-dimensional
spaces is the necessity to compute the integrand a very large number
of times. For instance, calculating $Z$ using Trapezoidal or
Simpson Rules with only 10 points per degree of freedom implies 
$10^{ 3 A}$
 grid points, and therefore function evaluations. Even for a
tiny system with $A = 20$ this corresponds to $10^{ 60}$ grid points!

This section deals with a method of integration which is based on a
random of exploration of phase space. It is dubbed a {\em Monte Carlo
Integration Method}. The name, of course, stems from the random or
chance character of the method in connection to the casino in Monte
Carlo. The method is based on a representative random sampling of
phase space. As such it is analogous to the prediction of the results
of an election based on a poll of a small number of voters.



\subsection{ The Monte Carlo Integration Strategy }

We will first look at an integral in one dimension to learn about the
Monte Carlo Method. Of course the method is normally used to  evaluate 
multi-dimensional integrals. 

Consider an integral over $x$  in a domain [0,1].
$$
I = \int_{0}^{1} f(x) dx
$$
An alternative to the Trapezoidal or
Simpson Integration Rules is reached by thinking of $I$ as resulting from the
product of the average of $f(x)$ over the domain by the
volume of the domain. Therefore $I$ can be approximated as
$$
     I \equiv \frac{1}{N} \sum_{i=1}^{N} f( x_i)
$$
$I$ is therefore proportional the average of $f(x)$ using a set of 
points $x_i$. 
A possible choice of points $x_i$ could be random positions with equal
probability everywhere within the domain [0,1]. 

The uncertainty
associated with this evaluation of the integral is given by
the variance in calculating the average of $f(x)$, namely
$$
   \sigma^2 = \frac{1}{N}  < ( f(x)  - < f(x) > )^2 >
$$
 or
$$
   \sigma^2 = \frac{1}{N} ( < f^2(x) > - < f(x) >^2 )
$$
 or 
$$
   \sigma^2 = \frac{1}{N}  \left[ \frac{1}{N} \sum_{i=1}^{N} f_i^2 -
                 \left(  \frac{1}{N} \sum_{i=1}^{N} f_i  \right)^2   \right]
$$
$\sigma$ measures the extent to which $f(x)$ deviates from its average
value over the region of integration. 

The calculation of $\sigma$ reveals
two important aspects of the Monte Carlo Integration. First, the
uncertainty in the estimate of the integral decreases as $N^{-1/2}$. If more
points are used to calculate $I$, better the answer will be. But this
is by far a slower convergence than the traditional Trapezoidal or Simpson Rules
that go as $1/N$ or $1/N^2$.

The second aspect  that $\sigma$ reveals is that the
error in calculating the integral depends on the fluctuations in $f(x)$. For
$f(x)$ constant, a small number of $N$ (1 value) would be enough to
calculate the average accurately. On the other hand, a very peaked $f(x)$
would be poorly represented by sampling with equally distributed
random numbers and would therefore produce a very poor value of the 
integral.


\subsection{ A better sampling }


The Monte Carlo Integration Method can be improved greatly by
a nonuniform sampling  if variations in the integrand $f(x)$ 
are smooth.

Multiply and divide the integrand $f(x)$ by a
positive weight function $w(x)$
$$
    I = \int_{0}^{1} dx w(x) \frac{ f(x) }{ w(x) }
$$
Let's assume further that $w(x)$ is normalized to 1 in the same domain
as the integral.
$$
        \int_{0}^{1} dx w(x) = 1
$$
Change variable from $x$ to $y(x)$ by
$$
       dy = w(x) dx
$$
 or 
$$
       y(x) = \int_{0}^{x} dx' w(x')
$$
Note that $ y(0) = 0 $ and $ y(1) = 1 $.

The integral now becomes
$$
     I = \int_{0}^{1} dy \frac{ f(x(y)) }{ w(x(y)) }
$$

The Monte Carlo evaluation of $I$ proceeds as before, e.g., averaging
the values of f(x) divided by $w(x)$ via a random sampling of points
uniformally distributed in $y$ in the interval [0,1], namely
$$
    I \equiv \frac{1}{N} \sum_{i=1}^{N} \frac{ f( x(y_i) ) }{ w( x(y_i) ) }
$$

If $w(x)$ has the same smooth variation as $f(x)$, e.g.,  $w(x)$
is large or small where $f(x)$ is large or small, the ratio $f/w$ will be
roughly a constant over the interval [0, 1]. This ought to improve the
accuracy in calculating the integral. 



There is a fundamentaly different interpretation of the
the results that we have just obtained. 
The change of variable in differential form
$$
       dy = w(x) dx
$$
 shows that a uniform distribution in $y$ ( since $w(y) = 1$ )
results in a distribution in $x$ given by the weight function $w(x)$. The
integral can therefore be interpreted as an average of $f(x)$ with
random numbers sampling according to the
probability weight function $w(x)$. In other words, sampling $f(x)$ with
random numbers distributed according to $w(x)$, the latter roughly
following the functional dependence of $f(x)$, will be much more accurate
than sampling with equally distributed random numbers.
 
We will see in the next subsection a method to generate random
numbers of arbitrary distribution. This will become a fundamental tool
to evaluate multi-dimensional integrals.
 
\subsection{ The Metropolis {\it Et Al} Algorithm }
 
To perform accurate Monte Carlo Integration  requires
sampling the function $f(x)$ using random numbers distributed according to
a positive definite weight function $w(x)$. The Metropolis {\it Et Al}
 algorithm
was devised to generate random numbers with arbitrary positive 
weight function $w(x)$. This algorithm, due to Metropolis, Rosenbluth,
Rosenbluth, Taller, and Taller, in 1953, has become a very important
tool for any Monte Carlo technique.
 
 Assume that $X$ is a
variable in the multi-dimensional space. We are looking for a
set of points ${X_i}$ in this phase space which follow a probability
density $w(X)$. The sequence of points $X_0$, $X_1$, $...$  corresponds 
to the points
visited by a guided random walker moving through this multi-dimensional
space as his walk becomes longer and longer.
 
The algorithm can be described as follows: assume the walker is
at the point $X_n$ in the sequence. The walker takes a tentative trial step
to a new point $X_t$. This point can be chosen in any convenient way, for
example using a uniform random choice within a multi-dimensional cube
of the small size $\delta$ about $X_n$. This trial step is then accepted or
rejected according to the ratio.
$$
         r = \frac{ w(X_t) } { w(X_n) }
$$

If $r > 1$ then the step is accepted ( $X_{n + 1} = X_t$). If $r < 1$ the trial
step is accepted with probability $r$. This is accomplished by comparing
$r$ with the random number $n$ uniformly distributed in the interval [0,1]
and accepting the step if $ r > n$, 
( $X_{n + 1} = X_t$). If the trial step is not accepted
then it is rejected,( $X_{n + 1} = X_n$),
  which means that it does not participate in taking
averages. 
Any arbitrary points $X_0$ can be used as the starting point
for the random walk. All steps are generated by the algorithm as
described.
 
It is possible to prove that this algorithm yields a distribution of
random numbers, set of points in multi-dimensional space,
according to the weight function $w(x)$.
 The proof is in one of the following subsections.
 
An interesting question is what to choose for the step size, $\delta$. This
choice is obviously arbitrary. Taking $\delta$ too large will result
in poor sampling since when near a maximum of $w(x)$ very few steps will
take the walker away from the maximum. If on the other hand the steps
are too small, the walker will stand litttle chance to explore
the space by going far from the initial point. 
Neither of these two choices yields efficient sampling of
the space. In practice $\delta$ is chosen so as to have roughly half of
the trial steps accepted by the criteria of the algorithm.
 
A great danger in sampling a multi-dimensional
 space via a random walk is to have
steps taken by the walker which are not independent from one
another. This could result from poor random number generation
algorithms or from the way the walk is built. If this happens then the
calculation of integrals by averaging over points in a random walk will
be innacurate.
 
 \subsection{ A Simple Guided Random Walk }


The Metropolis algorithm allows to generate a series of
states (or points in phase space
 if given a function of space coordinates)
distributed according to a given probability weight  function. 
This is illustrated in 
this section by generating a guided random walk based on the 
following probability profile:

$$ P(r) = r^2 \exp{( - r^2 )} $$

where $$ r( x, y ) = \sqrt{ x^2 + y^2 } $$ is the distance from the origin.
This is a gaussian with a hole punched out in its center. Load 
the following in Maple:
\begin{rawhtml}
<a href=../Simple_Guided_Walk/Probability_profile.mw>Probability_profile.mw</a>
\end{rawhtml}
. 
This weight function
 is simple yet it is difficult to generate
   random points corresponding to it directly. 
Therefore, we use the Metropolis algorithm to find sampling points in 
the $x - y$ plane according to the probability profile $P(r)$.

The program 
\begin{rawhtml}
<a href=../Simple_Guided_Walk/guided_walk.c>guided_walk.c</a>
\end{rawhtml}
generates the random walk starting from an arbitrary location in the 
plane. It implements the Metropolis algorithm to generate
 the random steps based 
on the given probability weight function. Note that the  latter
does no need to be normalized for the Metropolis algorithm to work.
This turns out to be a huge advantage of the method.

Use GNUPLOT to plot the location of the steps. Note that the 
points seem to be more numerous within the annulus as 
specified by the probability weight function $P(r)$.

It is instructive  to check the probability distribution of 
those generated steps as compared to the requested probability 
weight function. 
This requires calculating the density of points in the $x - y$ plane.
 A convenient way to bin the points  is via a set of concentric 
rings of equal areas centered onto the origin. This is so 
because the probability profile is axially symmetric. 
The number of points in each ring should follow the probability 
weight function. The program 
\begin{rawhtml}
<a href=../Simple_Guided_Walk/check_radial.c>check_radial.c</a>
\end{rawhtml}
bins the random steps according to the ring they belong to
as given by their distances from the origin. 
It generates a table to be used in a histogram plot. 
A subtle point has to be taken into account in displaying 
the bin contents at distances between the boundaries of 
the rings. They should be posted at a radius such that it divides 
the rings in equal ( half ) areas on each sides.

It is also interesting to analyze the capability of the Metropolis 
algorithm to generate random points exhibiting the 
axial symmetry of the probability weight function. The program 
\begin{rawhtml}
<a href=../Simple_Guided_Walk/check_centroid.c>check_centroid.c</a>
\end{rawhtml}
calculates the center of mass of those points. It should be at (0,0).
The distance from the origin is an elementary
 measure of the lack of quality of the random steps.

We can also calculate the axial distribution of the random steps. 
This is accomplished by calculating the number of points in wedges 
of equal areas centered onto the origin. The program 
\begin{rawhtml}
<a href=../Simple_Guided_Walk/check_angular.c>check_angular.c</a>
\end{rawhtml}
calculates the data to plot this histogram.

Use 3,000 steps and 50,000 steps to develop an intuition about
the Metropolis algorithm by comparing the results.

The Metropolis algorithm guides the random walk to generate a radial 
profile that is quite accurate. The angular profile, on the other hand,
 is not so accurate. This is because the probability profile has a 
localized radial dependence which is very effective in guiding 
the random walk. This is not the case in the angular direction 
since the probability profile is a constant in this direction 
due to its axial  symmetry.

The
\begin{rawhtml}
<a href=../Simple_Guided_Walk/guided_walk.c>guided_walk.c</a>
\end{rawhtml}
code uses an arbitrary step size, yielding a somewhat large acceptance
ratio.  The effectiveness of the Metropolis algorithm certainly
depends on this step size.  This should be investigated.  The
step size should be added to the line arguments list in the
code to simplify this analysis.







\subsection{ The Ising Model in Two Dimensions }
 
The Ising model in two dimension is a model for a magnetic
material.  It consist of classical spins on a equally spaced
lattice in 2 dimensions. Basic statistical 
physics will be used 
to computer the properties of such a model at temperature $T$.

The model consists of a set of spin degree of freedom interacting with
each other and with an external magnetic field. For example, these spins could
represent the magnetic moments of the atoms in a solid. The spins
are located on a 2D lattice of size $N_s = N_x N_y$. 
Each spin is labeled by $S_\alpha = S_{ij}$,
where $i$ and $j$ refer to the position on the
lattice, and $\alpha$ is a generic site label. These spins are
assumed classical (i.e., not following quantum mechanics rules) and
assumed to be  spin $\frac{1}{2}$. Therefore each spin is pointing up
($S_\alpha=+1$) or down ($S_\alpha=-1$).

The energy of the system is written
$$
       E(S) \equiv H(S) = - J' \sum_{ < \alpha \beta > } S_\alpha S_\beta
                     - B' \sum_\alpha S_\alpha
$$
The notation $ < \alpha \beta >$ means that 
the sum is over the nearest neighbor pairs
of spin.  The spin at site $i,j$
interacts with spins at $i \pm 1, j$ and $i, j \pm 1$ with strength $J'$. 
Each spin could also
interact with an external magnetic field of strength $B'$.

\begin{rawhtml}
        <center><img src=../Ising_Model/interaction.jpg height=400 width=400 border=0></center>
\end{rawhtml}

Periodic boundary conditions are applied.
The neighbor of a spin on the left of the lattice is taken as the
spin on the right of the lattice and vice versa for the spin on the
right. The same applies for the spins at sites at the bottom and the top of the
lattice. These boundary conditions model a system
of infinite dimensions. The topology of the lattice is
therefore that of torus.

If $J$ is positive the energy favors the parallel spin between
neighbors, this is {\it ferromagnetism}. If $J$ is negative, the spin will
tend to be anti-aligned or {\it antiferromagnetism}.

Any system in statistical equilibrium at temperature $T$ is such that
all possible configurations of the systems may occur with  probabilities
$$
     w(S) = \frac{ \exp( - \frac{ H(S) }{ k T } ) }{ Z }
$$
where the partition function  $Z$
$$
        Z( J', B' ) = \sum_{ S }  \exp( - \frac{ H(S) }{ k T } )
$$
normalizes $w(S)$. $T$ is the temperature in absolute degree and
$k$ the Boltzman constant. For our calculation it is
convenient to scale the strength of the interaction $J=\frac{J'}{kT}$
 and of the magnetic field $B=\frac{B'}{kT}$ in terms of the 
factor $kT$. Therefore high temperature
will correspond to  small $J$ and $B$ while at low
temperature $J$ and $B$ will be large.

The energy of the system depends on the $N_s = N_x N_y$ spins on the
lattice. Each spin can be in two possible configurations, $\pm 1$.
The total number of configurations is then $2^{N_s}$.
 This is an enormous number of terms to perform sums over. A very modest
lattices of $20,20$ corresponds to a number of configurations
$2^{400}$. Therefore the sums have to be done via 
Monte Carlo sampling. Note that
these sums take the role of the integrals we described in previous
sections.

The needed sums are based on the simple statistics idea of averaging the
given variables over the probabilities. The energy of the system is
therefore
$$
     E = \sum_S w(S)  H(S)
$$
The magnetization is
$$
     M = \sum_S w(S) \left(  \sum_\alpha S_\alpha \right)
$$
The specific heat at constant field is
$$
     C_B = \sum_S w(S)  H^2(S)   -  E^2
$$

The program 
\begin{rawhtml}
<a href=../Ising_Model/ising.c>ising.c</a>
\end{rawhtml}
implements this solution. It uses Monte Carlo
sampling obtained via the Metropolis et al algorithm based on the
probability of the spin configurations to calculate the energy
and magnetization of the system.

The code defaults to a {\it 20 X 20} lattice. It samples phase space by
flipping the spins at random. The issue becomes that of the step size as we
discussed before. In this case the choice is to flip spins or not to
flip spins. Flipping all the spins on the lattice at once
would likely lead to a step that is too large; configurations would be
missed in the sampling. The code uses the idea of a sweep of the
lattice where spins are systematically flipped from site to site
on the entire lattice with the proviso that each 
trial flip be accepted or rejected
according to the Metropolis algorithm. This is found by experience to
be a more effective way of sampling.

It also leads to an efficient way of calculating the ratio of
probability between the trial and the previous spin configuration. The
acceptance of the trial step depends on the ratio of the weight
function.
$$
         r = \frac{ w(X_t) } { w(X_n) } \equiv
                e^{-H(S_t) + H(S) }
$$
When only one spin is flipped the ratio of the exponential of the
energy depends only on the neighbors spins. The ratio takes the form
$$
     r = e^{ -2 S_\alpha( J f + B ) }
$$
where
$$
    f = S_{i+1,j} + S_{i-1,j} + S_{i,j+1} + S_{i,j-1}
$$
The $f$ factor can only take 5 different values,
namely $0, \pm 1, \pm 2, \pm4$. Only 10 values of $r$ can occur. This
makes for a very efficient coding by storing the exponential factors
ahead of time.

Another interesting question is to know when to stop sampling?
 This question is answered in the code by a special
way of doing statistic. Sampling is done at every NFREQ sweeps. These
sweeps are binned into groups with NSIZE members. For each group the
mean and the standard deviations are calculated. Groups values for the
average are then considered to be independent variables and a grand
average is calculated. The standard deviation of this grand average can be
calculated via the standard variance formula for the different averages
of the groups or via the average of the variances of the groups. When
sufficient sampling is accomplished, these two results should be the
same. This provides a measure of goodness of the sampling.

\subsection{ Exact solution of the Ising Model }

Onsager solved the Ising model analytically through a fantastic
application of mathematical physics.

Energy
$$
   E =  - N_s J \coth(2J) \left[ 1 + \frac{2}{\pi} \kappa' K_1(\kappa) \right]
$$


Magnetization
$$
    M = \pm N_s \frac{  ( 1+z^2)^{1/4} (1-6 z^2 + z^4 )^{1/8} }
                     {  ( 1-z^2)^{1/2} }
$$


Specific heat
$$
   C_B = N_s \frac{2}{\pi} (  J \coth(2J) )^2 \left(
        2  K_1(\kappa) - 2  E_1(\kappa)  - ( 1 - \kappa' )
               \left[  \frac{\pi}{2}  + \kappa'  K_1(\kappa) \right]
                                              \right)
$$


$z$, $\kappa$ and $\kappa'$ are
$$
    z = exp(- 2 J)
$$
$$
    \kappa = 2 \frac{ \sinh(2 J ) } { \cosh^2( 2 J ) }
$$
$$
    \kappa' = 2 \tanh^2( 2 J ) - 1
$$

The elliptic functions are
$$
   K_1(\kappa) \equiv \int_{0}^{\pi/2} \frac{ d \phi }
               { ( 1 - \kappa^2 \sin^2( \phi) )^{1/2} }
$$
$$
   E_1(\kappa) \equiv \int_{0}^{\pi/2}  d \phi
                ( 1 - \kappa^2 \sin^2( \phi) )^{1/2}
$$

The analytical solution of the model
in 3D is not known, and it may be impossible to derive.

\subsection{ Proof of the Metropolis Algorithm }

To prove that the Metropolis algorithm generates a sequence of
random numbers distributed according to $w(x)$
 consider a large number of walkers
starting from different initial points and moving independently.
If $N_n(X)$ is the density of those walkers at $X$ after $n$ steps
then the net number of walkers moving from point $X$ to point $Y$ in the
next step is
$$
      \Delta N(X) = N_n(X) P( X \to Y ) - N_n(Y) P( Y \to X )
$$
$$
       \Delta N(X) = N_n(Y) P( X \to Y )
              \left[ \frac{ N_n(X) }{ N_n(Y) }  -
                             \frac{ P( Y \to X ) }{ P( X \to Y ) }  \right]
$$
Here $P(X \to Y)$ is the probability that the walker 
will make a transition to $Y$
if he is at $X$, and the inverse for $P(Y \to X)$. 
The equation shows equilibrium when
$$
    \frac{ N_n(X) }{ N_n(Y) } = \frac{ N_e(X) }{ N_e(Y) } =
             \frac{ P( Y \to X ) }{ P( X \to Y ) }
$$
Note that $N(X)$ changes when the system is not in equilibrium via the
weight that tends to drive it toward equilibrium



To show that the algorithm leads
to an equilibrium distribution of walkers proportional to $w(X)$ take
the probability of making a step to $Y$ as
$$
     P( X \to Y ) =  T( X \to Y )  A( X \to Y )
$$
where $T$ is the probability of making a trial step from $X$ to $Y$ and 
$A$ is
the probability of accepting that step. If $Y$ can be reached from $X$ in
a single step, i.e. it is within cube of side $\delta$ centered about $X$, then
$$
      T( X \to Y )  =  T( Y \to X )
$$
so that the equilibrium distribution of the metropolis random walkers
satisfy.
$$
          \frac{ N_e(X) }{ N_e(Y) } =
             \frac{ A( Y \to X ) }{ A( X \to Y ) }
$$

From the algorithm,, if $w(X) > w(Y)$, then $A( Y \to X ) = 1 $ and
$$
         A( X \to Y ) = \frac{ w(Y) }{ w(X) }
$$
If $w(X) < w(Y)$, then $A( X \to Y ) = 1 $ and
$$
         A( Y \to X ) = \frac{ w(X) }{ w(Y) }
$$
Hence,, in either case the equilibrium population of the metropolis
random walkers satisfy
$$
          \frac{ N_e(X) }{ N_e(Y) } = \frac{ w(X) }{ w(Y) }
$$

The walkers are distributed with the correct
distribution. This proves that the walkers will generate a
distribution of points ${X_i}$ in space which will be distributed
according to the weight function $w(X)$.

\subsection{ Monte Carlo Codes }

Demo codes - integral  of
$$
   I = \int ( x - \frac{1}{4} )^2 \exp( - \frac{1}{2} x^2 )
$$
\begin{enumerate}

\item by Trapezoidal Rule
\begin{rawhtml}
<a href=../Ising_Model/sample_codes/int_trap.c>int_trap.c</a>
\end{rawhtml}

\item Monte Carlo Method
\begin{rawhtml}
<a href=../Ising_Model/sample_codes/int_mc.c>int_mc.c</a>
\end{rawhtml}

\item and Guided Monte Carlo Method
\begin{rawhtml}
<a href=../Ising_Model/sample_codes/int_guided_mc.c>int_guided_mc.c</a>
\end{rawhtml}

\end{enumerate}


Ising Model solution: 
\begin{rawhtml}
<a href=../Ising_Model/ising.c>ising.c</a>
\end{rawhtml}




\end{document}
