Maxwell’s distribution using Box-Muller transform

normal distributionprobabilityprobability distributions

I am given the following problem to solve (This text is translated from Russian. So, there may be some translation issues):

… Another method to draw from the normal distribution is to draw two independent random numbers from the uniform distribution x1, x2 ∈ [0:0, 1:0). Then apply the following transformation:
$n_1 = \sqrt{-2\log{x_1}}\cos{2\pi x_2}$
$n_2 = \sqrt{-2\log{x_1}}\sin{2\pi x_2}$
resulting in two randomly independent numbers n1, n2 from a normal distribution with zero expected value and unit variance.

To change the distribution parameters to other parameters, e.g. the expected value for and the variance to, you should multiply the result of the draw by and add, i.e.
$N(\mu, \sigma) = \sigma N(0, 1) + \mu $
In the equation above, N(μ, σ) is a random variable with normal distribution with expected value μ and variance σ.

According to the Maxwell distribution, each component (x, y or z) of the velocity vector v is a random variable from a normal distribution with zero expected value, and variance
$\sqrt{\dfrac{k_B T}{m}}$
where m is the mass of the molecule, T is the temperature in Kelvin, kB is Boltzmann constant.

Your task: Draw 10,000 velocity vectors for the N2 nitrogen molecule at 300K. Calculate the average length of these vectors, and therefore the average value of the speed of the nitrogen molecule, using the formula:
$\bar{v} = \dfrac{1}{N}\sum_{i=1}^N\sqrt{v_{x_i}^2 + v_{y_i}^2 + v_{z_i}^2}$

I am not understanding the relationship of N(μ, σ) with n1 and n2, and how to get to vector v from there.

Can anyone explain?

Best Answer

Firstly, $n_1,\,n_2$ are $N(0,\,1)$ iids. Suppose $v_x,\,v_y,\,v_z$ are $N(0,\,a^2)$ iids, e.g. obtained by multiplying samples of $N(0,\,1)$ by a. By converting from Cartesian to spherical polar coordinates, you can easily show $v:=\sqrt{v_x^2+v_y^2+v_z^2}$ is Maxwell-Boltzmann distributed (I've borrowed the parameter $a$ from here). In particular, $v$ has a PDF proportional to $v^2\exp\frac{-v^2}{2a^2}$ on $\Bbb R^+$. Of course we want proportionality to $v^2\exp\frac{-mv^2}{2k_BT}$, so $a=\sqrt{\frac{k_BT}{m}}$. So:

  • Get $30000$ independent samples of $U(0,\,1)$ into $15000$ pairs
  • Use each such pair to make two $N(0,\,1)$ samples, so we have $30000$ of them
  • Multiply them by $\sqrt{\frac{k_BT}{m}}$
  • Group them into $10000$ triples
  • Calculate $v$ from the triples as above, giving $10000$ of them, then average

We can be more efficient by multiplying by $a$ at the end. A not particularly efficient implementation:

from numpy import cos, log, mean, pi, sin, sqrt
from numpy.random import random

twopi = 2*pi
m, k, T = 14*1.67e-27, 1.38e-23, 300
a = sqrt(k*T/m)

def box_muller():
    r, theta = sqrt(-2*log(random())), twopi*random()
    return [r*cos(theta), r*sin(theta)]

normals = [z for _ in range(15000) for z in box_muller()]
vs = [sqrt(sum(v*v for v in normals[3*i:3*i+3])) for i in range(10000)]
print(f'In SI units, a={a} while vbar={a*mean(vs)}')

You can actually show with some calculus that the exact result is $\sqrt{\frac{8k_BT}{\pi m}}$, which is very close to what we get.

Related Question