I asked on Math.SE and was advised to try here instead.

I need to draw from a Maxwell-Boltzmann velocity distribution to initialise a molecular dynamics simulation. I have the PDF but I'm having difficulty finding the correct CDF so that I can make random draws from it.

The PDF I am using using is:

$$f(v)=\sqrt \frac{m}{2\pi kT} \times exp \left( \frac{-mv^2}{2kT}\right) $$

I am told that to find the CDF from the PDF we perform:

$$CDF(x)= \int_{-\infty}^x PDF(x) dx $$

After integrating $ f(v) $ I get:

$$ CDF(v)= \sqrt \frac{m}{2\pi kT} \times \left( \frac{\sqrt\pi\times erf \left( \frac{mv}{2\pi kT} \right) }{2\times \left( \frac{m}{2kT} \right)} \right) $$

$$CDF(v)= _{-\infty} ^{x} \left[ {\sqrt \frac{m}{2\pi kT} \times \left( \frac{\sqrt\pi\times erf \left( \frac{mv}{2\pi kT} \right) }{2\times \left( \frac{m}{2kT} \right)} \right)} \right] $$

  1. After I reach this point I am unable to proceed as I do not know how to evaluate something between $x$ and ${-\infty}$.

  2. I am also concerned that I have not done the integration correctly.

  3. I want to implement the CDF in C++ in the end so I can draw from it. Does anyone know if there will be a problem with doing this because of the erf, or will I be alright with this GSL implimentation ?

Thanks for your time.

@bryansis2010 on Math.SE says that I can evaluate in the range $x$ to $0$ instead of $-\infty$ as we do not drop below 0 Kelvin.

Would this then make the CDF:

$$ CDF(v)= \sqrt \frac{m}{2\pi kT} \times \left( \frac{\sqrt\pi\times erf \left( \frac{mv}{2\pi kT} \right) }{2\times \left( \frac{m}{2kT} \right)} \right) $$

as $erf(0)=0$

Is this correct?

  • 2
    $\begingroup$ Hi RRs_Ghost - I suspect that whoever told to bring this question here was incorrect. I don't think it's on topic for us. Even though the function you're integrating comes from physics, you're still just asking how to do an integration, which itself is a pure math problem. I won't close this immediately so that other people have a chance to object if they would like to, but we may just be sending you back to Mathematics. (Also, for future reference: you shouldn't cross-post a question to multiple SE sites. Ask for it to be migrated if it's off topic on the first place you put it.) $\endgroup$ – David Z Feb 17 '13 at 16:34

The solution is to realise that that function is merely a Gaussian. In fact Each component of the velocity vector has a normal distribution with mean =0 and st-dev $\sqrt {kT/m}$.

All that is left to do at that point is to get the Gaussian CDF (well known) and sample from it, making sure to plug in our mas and temperature.

$$CDF(x)=\frac{1}{2}\times \left[ 1+ erf\left( \frac{x-\mu}{\sqrt{2{\sigma}^2}} \right) \right]$$

GSL implements gsl_cdf_ugaussian_P (double x) here


Not the answer you're looking for? Browse other questions tagged or ask your own question.