Calculating pi is all well and good, but not particularly useful. Let's consider a more involved example.
Let's take a box which has some atoms inside. The atoms are in the vapour phase, and we can consider it an ideal gas. In this case the Maxwell-Boltzmann distribution tells us how fast the atoms are moving. The speed distribution obeys the relationship,
$$ f_v(v) = \left(\frac{m}{2\pi k_B T}\right)^{3/2} \; 4\pi v^2 \; \exp({-mv^2/2k_B T}) $$,
that is, the fraction of atoms $f_v$ with velocity $v$, $k_B$ is Boltzmann's constant, $T$ is the thermodynamic temperature (in Kelvin) and $m$ is the mass of the atoms.
Let's have a look at that distribution graphically (it should be familiar!):
In [54]:
def MB_speed(v,m,T):
""" Maxwell-Boltzmann speed distribution for speeds """
kB = 1.38e-23
return (m/(2*np.pi*kB*T))**1.5 * 4*np.pi * v**2 * np.exp(-m*v**2/(2*kB*T))
fig = plt.figure()
ax = fig.add_subplot(111)
v = np.arange(0,800,1)
amu = 1.66e-27
mass = 85*amu
for T in [100,200,300,400]:
fv = MB_speed(v,mass,T)
ax.plot(v,fv,label='T='+str(T)+' K',lw=2)
ax.legend(loc=0)
ax.set_xlabel('$v$ (m/s)')
ax.set_ylabel('PDF, $f_v(v)$')
plt.draw()
In atomic physics experiments with thermal vapours, we use lasers to probe the atoms. The resonant frequency of (an electronic transition in) the atoms is very well defined and has a Lorentzian distribution centred on the resonant frequency and with a width that is related to the radiative decay rate of the transition.
The laser we use to probe the atomic resonance can be considered to be monochromatic, i.e. it has a single frequency and travels along a single axis, which we normally take to be the z-axis. The Doppler effect due to the motion of the atoms relative to the laser beam introduces a shift of the resonance frequency, which broadens the spectral line.
We know the speed distribution of the atoms, but we want to know the distribution of velocities in the z-axis, since this is the only axis that will produce a Doppler shift.
We will set up a Monte-Carlo simulation of atomic speeds, with velocity vectors uniformly distributed over all angles. Let's take Rubidium (Rb) atoms, with $m=87$ a.m.u., at a temperature of 100C, or 373K.
How do we generate the random speeds? There are a few ways of doing this, but one of the most general is using the cumulative distribution function (CDF), $C_v(v)$, which can be found from the probability density function (PDF, which for this case is $f_v(v)$), as follows:
$$ C_v(v) = \int_{-\inf}^{v} f_v(v') \; {\rm d}v' $$This can be done numerically if required, but fortunately in this case the solution is analytic, and is given by
$$ C_v(v) = {\rm erf}\left(\frac{v}{\sqrt{2k_B T/m}}\right) - \sqrt{\frac{2}{\pi}} \frac{ v \; \exp(-mv^2/2k_B T) }{\sqrt{k_B T / m}} $$where ${\rm erf}$ is the Error function. Let's look at that now:
In [55]:
from scipy.special import erf
def MB_CDF(v,m,T):
""" Cumulative Distribution function of the Maxwell-Boltzmann speed distribution """
kB = 1.38e-23
a = np.sqrt(kB*T/m)
return erf(v/(np.sqrt(2)*a)) - np.sqrt(2/np.pi)* v* np.exp(-v**2/(2*a**2))/a
fig = plt.figure()
ax = fig.add_subplot(111)
v = np.arange(0,800,1)
amu = 1.66e-27
mass = 85*amu
for T in [100,200,300,400]:
fv = MB_CDF(v,mass,T)
ax.plot(v,fv,label='T='+str(T)+' K',lw=2)
ax.legend(loc=0)
ax.set_xlabel('$v$ (m/s)')
ax.set_ylabel('CDF, $C_v(v)$')
plt.draw()
To generate our random speeds, we take a set of random numbers uniformly distributed between 0 and 1 and use the inverse of the CDF to get out a velocity. Because of the shape of the CDF, the uniform [0,1] range is mapped onto the velocity distribution. For example, the range [0,0.05] maps to a 100m/s (0 - 100 m/s) spread in velocities, whilst the range [0.2,0.4] (at 400 K) roughly maps a similar 100m/s (200 - 300 m/s) spread in velocity space. With this coarse approximation, there is a 5% chance to have a velocity between 0 and 100 m/s, whereas there is a 20% chance of it being between 200 and 300 m/s.
The main issue is that there is no analytic inverse CDF function, so we are forced to do it numerically. We will do this using an interpolation function, but there are other options available - see, e.g. http://www.scientificarts.com/maxwell/maxwellpdf.html or http://astro.physics.ncsu.edu/urca/course_files/Lesson20/index.html.
Fortunately, the scipy module includes tools for creating interpolation functions and is straightforward to do:
In [59]:
from scipy.interpolate import interp1d as interp
amu = 1.66e-27
mass = 85*amu
T = 400
# create CDF
vs = np.arange(0,2500,0.1)
cdf = MB_CDF(vs,mass,T) # essentially y = f(x)
#create interpolation function to CDF
inv_cdf = interp(cdf,vs) # essentially what we have done is made x = g(y) from y = f(x)
# this can now be used as a function which is
# called in the same way as normal routines
Now we can generate our speeds, just by calling inv_cdf. We will write it as a function so it can be called later on. In this function, we also generate the random angles for our velocity vectors, then convert this into cartesian coordinates.
In [114]:
def generate_velocities(n):
""" generate a set of velocity vectors in 3D from the MB inverse CDF function """
rand_nums = np.random.random(n)
speeds = inv_cdf(rand_nums)
# spherical polar coords - generate random angle for velocity vector, uniformly distributed over the surface of a sphere
# see http://mathworld.wolfram.com/SpherePointPicking.html for more info (note theta and phi are the other way around!)
theta = np.arccos(np.random.uniform(-1,1,n))
phi = np.random.uniform(0,2*np.pi,n)
# convert to cartesian units
vx = speeds * np.sin(theta) * np.cos(phi)
vy = speeds * np.sin(theta) * np.sin(phi)
vz = speeds * np.cos(theta)
return speeds, vx, vy, vz
Now we just need to test the code. Let's first make sure that our generated speed distribution matches what we expect from $f_v(v)$:
In [115]:
spd, vx, vy, vz = generate_velocities(100000)
fig = plt.figure()
ax = fig.add_subplot(111)
#generate histogram of velocities
ax.hist(spd,bins=50,normed=1,fc='b',alpha=0.4,lw=0.2)
#compare this histogram to f(v) - this is MB_speed that we wrote earlier
vs = np.arange(0,1500)
kB = 1.38e-23
amu = 1.66e-27
mass = 85*amu
T = 400
fv = MB_speed(vs,mass,T)
ax.plot(vs,fv,'k',lw=2)
ax.set_xlabel('Speed (m/s)')
ax.set_ylabel('PDF')
Out[115]:
Great - that agrees pretty well! Now let's look at the component in the z-axis, which is what we were interested in:
In [122]:
fig = plt.figure()
ax = fig.add_subplot(111)
h,b,c = ax.hist(vz,bins=50,normed=1,fc='r',alpha=0.4,lw=0.2)
print h.sum()*(b[1]-b[0])
ax.set_xlabel('$v_z$ (m/s)')
ax.set_ylabel('Probability density')
plt.draw()
So the component of velocity in the z-axis is symmetric about 0, i.e. there are just as many atoms moving away from the laser as towards it, so our resonance frequency will be broadened symmetrically.
Of course, this has an analytic form which is easy enough to derive (see wikipedia, for example), but hopefully this excercise has shown the principles behind Monte Carlo. A non-trivial application of this kind of simulation is in calculating the fraction of atoms that will be slowed in a magnetic Zeeman slower, which are used extensively in ultra-cold atomic physics experiments (see here for an example).
Let's add the analytic form on now; it has a Gaussian (normal) form:
In [124]:
u = np.sqrt(2*kB*T/mass)
vzt = np.arange(-4*u,4*u,u/50)
mbz = np.exp(-vzt**2/u**2)/np.sqrt(u*np.pi)
mbz /= (mbz.sum() * (vzt[1]-vzt[0]))
ax.plot(vzt,mbz,'k--',lw=2)
ax.set_xlim(-1000,1000)
display(fig)
MC lecture slides: http://www.phys.ubbcluj.ro/~zneda/edu/mc/mcshort.pdf
An alternate MC method for calculating pi: http://tommyogden.com/pi/
Ch 10 of Computational Physics, Jos Thijssen, CUP (2007)