Random numbers with Maxwell distribution in C

user_8740387 0
Viewed 2268

I need a 10000 random number that distributed with the Maxwell distribution. For the normal distribution, I know that I have to use Box-Muller transform but my question is that why 

 normal_distribution

is defined as a variable (or anything that I don't know) by default? Is "normal_distribution" a function that gives numbers with normal distribution? If yes, is this possible for Maxwell distribution? If not, what should I do? In fact, I want to learn how to create random numbers with Maxwell distribution in C. Thanks for any tips.

distribution numerical-methods random c
1 Answers

Your good bet is to use GNU Scientific Library. It is reasonable written and tested computational library in pure C. Maxwell distribution is basically three normally distributed components squared and summed together.

Some untested code:

#include <math.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>

double sample_maxwell(gsl_rng* r, double sigma) {
    double vx = gsl_ran_gaussian_ziggurat(r, sigma);
    double vy = gsl_ran_gaussian_ziggurat(r, sigma);
    double vz = gsl_ran_gaussian_ziggurat(r, sigma);

    return sqrt(vx*vx + vy*vy + vz*vz);
}

int main() {
    gsl_rng_env_setup();

    const gsl_rng_type* T = gsl_rng_default;
    gsl_rng*            r = gsl_rng_alloc(T);

    printf ("generator type: %s\n", gsl_rng_name (r));
    printf ("seed = %lu\n", gsl_rng_default_seed);
    printf ("first value = %lu\n", gsl_rng_get (r));        

    double Temperature = 300.0;   // K
    double kBoltzmann  = 8.62e−5; // eV/K
    double mass        = 1.0e+9;  // 1GeV/c^2, roughly atomic hydrogen

    double sigma = sqrt(kBoltzmann*Temperature/mass);

    for(int k = 0; k != 100; ++k) {
        double v = sample_maxwell(r, sigma);
        printf("%e", v);
    }

    gsl_rng_free(r);        

    return 0;
}

Maxwell distribution came from statistical physics, where molecule of gas has speed v and Boltzmann exponent comes to velocity distribution

f(v)d3v = C exp(- m v2/2 kBT) d3v

To get Maxwell distribution from this one, you could just express distribution in spherical over v coordinates, power term comes from Jacobian, basically. But one can easily see it is equivalent to product of three Gaussian distribution for each component. I'll refer you to the following page https://farside.ph.utexas.edu/teaching/sm1/Thermalhtml/node87.html and formula 7.217 in particular. Here is another link to friendly discussion https://scicomp.stackexchange.com/questions/19969/how-do-i-generate-maxwell-boltzmann-variates-using-a-uniform-distribution-random.

There are alternative ways to sample Maxwell:

  1. Recognize, that sum of three squared Gaussian could be expressed as Chi2 distribution with 3 degrees of freedom and sample v2 via gsl_ran_chisq(r, 3).

  2. Sample v2 as result of Gamma distribution via gsl_ran_gamma_knuth(r, 3./2., 1.0).

And there is always Wikipedia where all of this is also stated: https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution

user_40446960
Related