\documentstyle[12pt]{article}
\begin{document}
\title{Monte Carlo Methods}
\author{Anne Greenbaum \thanks{Math/AMath 381, Winter quarter, 1998.
Ref. {\em Monte Carlo Methods}, by Kalos and Whitlock, Wiley, 1986.} }
\maketitle
\section{Examples of Monte Carlo Computations}
\begin{itemize}
\item
Compute $\frac{\pi}{4}$.
\vspace{1in}
\item
Buffon needle problem: Suppose a sheet of paper has lines a distance
$d$ apart. Drop a needle of length $L < d$ onto the paper. What is
the probability that the needle intersects a line? The answer is
now known analytically; it is $\frac{2 L}{\pi d}$. It is reported
that Buffon conjectured this value and actually did the experiment
many times to confirm it before devising a proof.
\vspace{1in}
\pagebreak
\item
Multidimensional integration:
\[
\int_0^1\!\!\int_0^1 \ldots \int_0^1 f( x_1 , \ldots , x_m )
\,d x_1\,d x_2 \ldots\,d x_m
\]
A standard quadrature rule would approximate this integral using a formula
of the form
\[
( \Delta x )^m \sum_{i_1 =1}^{n} \sum_{i_2 =1}^{n} \ldots \sum_{i_m =1}^{n}
f( x_{1_{i_1}} , x_{2_{i_2}} , \ldots , x_{m_{i_m}} ) ,
\]
and the work required would be $O( n^m )$. This is infeasible, say, for
$n=10$ and $m=30$.
In a Monte Carlo computation, we choose many independent samples from
$[0,1] \times [0,1] \times \ldots \times [0,1]$ and approximate the
integral by
\[
\frac{1}{N} \sum_{k=1}^{N} f( \vec{x}_k ) .
\]
The error is proportional to $N^{-1/2}$. This is rather slow convergence;
we must increase the number of points by a factor of $4$ in order to
obtain a factor of $2$ reduction in the error. But for $m$ large, it
is the only feasible method. Most Monte Carlo calculations can be
interpreted as evaluating a multidimensional integral.
The Monte Carlo computation of $\pi$ can be thought of as a 2-dimensional
integral: Evaluate the integral over the square $[-1,1] \times [-1,1]$
of a function $f(x,y)$ that is equal to $1$ if $x^2 + y^2 \leq 1$ and $0$
otherwise.
Areas and volumes of more complicated regions are easily
computed provided only that we can enclose the region in a box or other
object of known volume and then determine whether a point lies inside
or outside the region. For example, suppose we wish to find the
weight and center of mass of a region that is the intersection of
a torus with the edge of a large box. [Ref. {\em Numerical Recipes}
by Press, et.al., Cambridge University Press, 1992.] Let the region
be defined by the three simultaneous equations
\[
z^2 + \left( \sqrt{x^2 + y^2} - 3 \right)^2 \leq 1 ,~~~x \geq 1 ,~~~
y \geq -3 ,
\]
and assume that it has constant density $\gamma$.
We want to estimate the following integrals over the the complicated
region:
\[
\int \gamma\,dx\,dy\,dz ,~~~
\int x \gamma\,dx\,dy\,dz ,~~~
\int y \gamma\,dx\,dy\,dz ,~~~
\int z \gamma\,dx\,dy\,dz .
\]
The coordinates of the center of mass will be the ratio of the latter
three integrals to the first one (the weight).
The following piece of MATLAB code will provide such estimates:
\begin{verbatim}
w = 0; %Initialize weight and other integrals to zero.
wx = 0; wy = 0; wz = 0;
for i=1:n,
x = 1 + 3*rand; %Generate random points uniformly distributed
y = -3 + 7*rand; %in the enclosing box [1,4] X [-3,4] X [-1,1].
z = -1 + 2*rand;
if z^2 + ( sqrt(x^2+y^2) - 3 )^2 <= 1, % Is point in the torus?
w = w + gamma; % If so, add to various integrals.
wx = wx + x*gamma;
wy = wy + y*gamma;
wz = wz + z*gamma;
end;
end;
boxvol = 3*7*2; %Volume of enclosing box.
weight = boxvol*w/n, %Integrals are volume of box times fraction
xmoment = boxvol*wx/n, %inside region.
ymoment = boxvol*wy/n,
zmoment = boxvol*wz/n,
\end{verbatim}
This is fine, but the computed estimates of \verb+weight+, etc.
are not so useful unless we have some indication of how accurate they are.
Here we use the central limit theorem: For $n$ sufficiently large, the
random variables \verb+weight+, etc. produced by this code are normally
distributed with means equal to the actual desired integral values
and standard deviations given by $\sigma / \sqrt{n}$, where $\sigma$
is the standard deviation of the individual distributions. Since the
variance of a random variable $X$ is $E( X^2 ) - ( E(X) )^2$, these
standard deviations can be approximated, along with the integrals.
For example, to compute \verb+weight+ and also provide an error estimate
the following code could be used:
\begin{verbatim}
w = 0; %Initialize weight to zero.
varw = 0; %Initialize variance estimate to zero.
for i=1:n,
x = 1 + 3*rand; %Generate random points uniformly distributed
y = -3 + 7*rand; %in the enclosing box [1,4] X [-3,4] X [-1,1].
z = -1 + 2*rand;
if z^2 + ( sqrt(x^2+y^2) - 3 )^2 <= 1, % Is point in the torus?
w = w + gamma; % If so, add to w.
varw = varw + gamma^2; % Add to estimates of E(w^2).
end;
end;
boxvol = 3*7*2; %Volume of enclosing box.
weight = boxvol*w/n, %Integral is volume of box times fraction
%inside region.
sigw = boxvol*sqrt((varw/n - (w/n)^2)/n), %Standard deviation for weight
\end{verbatim}
Can you see how to compute similar error estimates for the other integrals?
If the above code were run many, many times, how often would you expect
the computed value of \verb+weight+ to lie within \verb+sigw+ of the
actual weight?
\end{itemize}
\section{Application: Radiation Transport}
Assumptions:
\begin{enumerate}
\item
Particles move in straight lines between collisions (true for neutral
radiation --- light, X-rays, neutrons; not for electrons and protons).
\item
Collisions result in absorption or scattering, with some loss of energy.
\item
Neglect effect of radiation on the medium.
\item
Ignore effects of polarization on scattering.
\item
Neglect crystal effects.
\end{enumerate}
Example: Model X-ray photons in a dental X-ray. Determine dose of
radiation received at a particular position; determine fraction of
radiation that passes through a slab of matter of known thickness.
Independent variables: position $(x,y,z)$, direction $( \Omega_x ,
\Omega_y , \Omega_z )$, energy $(E)$, and time of production $(t)$
of photons. Each independent variable is random and may be described
by a pdf.
Once in matter, a photon may interact with an individual atom.
Whether this interaction takes place and with which atom is random.
The interaction usually involves the photon exchanging energy with
an electron, effectively being scattered with reduced energy and
changed direction. The photon may also give up all its energy to
the electron, in which case it is considered to be absorbed and
the history terminates. We must posit pdf's for each stochastic
event, based on theoretical and measured properties of X-rays.
Scattered photons are subject to the same processes as the original
ones with the same energy. The life history of a photon is followed
until it is absorbed or too low in energy or too far away to have an
appreciable effect on the required answer.
\subsection{Characterization of the Source}
Assume a point source $( x_0 , y_0 , z_0 )$ and assume that each photon
produced initially has energy $E_0$. (This is not quite right but
close enough.) X-rays are produced by turning on an electron beam
at time $t_1$ and turning it off at time $t_2$. Assume a uniform
distribution of production times between $t_1$ and $t_2$. Take the
intensity of the X-rays to be uniform in all directions. The initial
state of a photon can be computed as follows:
\begin{verbatim}
x = x0; y = y0; z = z0; E = E0;
t = t1 + (t2-t1)*rand;
phi = 2*pi*rand;
omegaz = 1 - 2*rand;
omegax = sqrt(1 - omegaz^2)*cos(phi); omegay = omegax*tan(phi);
\end{verbatim}
\subsection{Tracing a Path}
Since X-rays travel in straight lines between collisions:
\begin{verbatim}
xnew = x + omegax * d;
ynew = y + omegay * d;
znew = z + omegaz * d;
tnew = t + d/v;
\end{verbatim}
\noindent
where $v$ is the radiation velocity and $d$ is the distance to the
next interaction with an atom. The probability per unit pathlength
of having an interaction, $\Sigma_T$, is a property of the material
and does not change with the distance the photon has traveled.
Consequently, the distribution of distances $d$ is exponential in
homogeneous materials:
\[
F(d) = \exp ( - \Sigma_T d ) \equiv \xi \Rightarrow d = - \log ( \xi ) /
\Sigma_T .
\]
In MATLAB notation, insert the line:
\begin{verbatim}
d = -log(rand)/sigmat;
\end{verbatim}
\subsection{Modeling a Collision}
There are two possibilities --- absorption or scattering. If
$\xi < \mbox{Prob(absorption)}$, absorption occurs; otherwise scattering.
Both neutrons and photons tend to scatter preferentially in the forward
direction; the azimuth is uniformly distributed.
\vspace{1in}
\noindent
Example:
\[
\rho ( \cos \theta )\, d \cos \theta = \left\{ \begin{array}{ll}
2 \cos \theta\, d \cos \theta , & 1 > \cos \theta > 0 \\
0 , & \cos \theta < 0 \end{array} \right.
\]
\[
\phi = 2 \pi \xi .
\]
\end{document}