Random numbers with Maxwell distribution in C

Question

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

 normal_distribution

Is defined as a variable (or anything 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.

Answer 1

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)d³v = C exp(- m v²/2 k_BT) d³v

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 Chi² distribution with 3 degrees of freedom and sample v² via gsl_ran_chisq(r, 3).

2. Sample v² 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