2
\$\begingroup\$

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}\$

using System;

public class CommonDistributions
{
    public static int Uniform(Random random, int n)
    {
        return (int)(random.NextDouble() * n);
    }

    public static double Uniform(Random random, double lo, double hi)
    {
        return lo + random.NextDouble() * (hi - lo);
    }

    public static double Gaussian(Random random)
    {
        double r, x, y;
        do
        {
            x = Uniform(random, -1.0, 1.0);
            y = Uniform(random, -1.0, 1.0);
            r = x * x + y * y;
        }
        while (r >= 1 || r == 0);

        return x * Math.Sqrt(-2 * Math.Log(r) / r);
    }

    public static double Gaussian(Random random, double mu, double sigma)
    {
        return sigma * Gaussian(random) + mu;
    }
}

public class MaxwellBolzman
{
    static double KB = 1.38064852e-23;

    static double MaxwellVariance(double mass, double temperature)
    {
        return Math.Sqrt(KB * temperature / mass);
    }

    static double MaxwellComponent(Random random, double mass, double temperature)
    {
        double mu = 0.0;
        double sigma = MaxwellVariance(mass, temperature);

        return CommonDistributions.Gaussian(random, mu, sigma);
    }
    public static double Maxwell(Random random, double mass, double temperature)
    {
        double one = MaxwellComponent(random, mass, temperature);
        double two = MaxwellComponent(random, mass, temperature);
        double thr = MaxwellComponent(random, mass, temperature);

        return Math.Sqrt(one * one + two * two + thr * thr);
    }
}

public static class MainClass
{
    public static void Main(String[] args)
    {
        Random random = new Random();

        const int N = 10000;
        const int T = 300;//300K
        const double mass = 28.02;//28.02 g/mol

        double sum = 0.0;

        for (int i = 1; i < N; i++)
        {
            sum = sum + MaxwellBolzman.Maxwell(random, mass, T);
        }

        Console.WriteLine($"Maxwell-Boltzman Vector = {sum/N}");

        string str = string.Empty;
    }
}

Kindly, review the implementation.

I am not sure about the values of temperature and the mass of Nitrogen 2.

\$\endgroup\$
4
  • 2
    \$\begingroup\$ You have some very overheated gas! 300K means 300 degrees on Kelvin scale (about 27 centigrades), not 300000. \$\endgroup\$
    – vnp
    Apr 15, 2020 at 2:59
  • \$\begingroup\$ Please do not update the code in your question to incorporate feedback from answers, doing so goes against the Question + Answer style of Code Review. This is not a forum where you should keep the most updated version in your question. Please see what you may and may not do after receiving answers. \$\endgroup\$
    – Mast
    Apr 15, 2020 at 11:32
  • \$\begingroup\$ @Mast, is it better now? \$\endgroup\$
    – user366312
    Apr 15, 2020 at 11:36
  • \$\begingroup\$ No. Rule of thumb: the moment answers start coming in, you stop touching the code. Comments? Fine. Answers? No. No more code edits. If you have new code that you'd like reviewed, ask a follow-up question instead (this should all be explained in the link, please comment if anything is left unclear). I'd encourage you to wait at least 24h between asking questions though. \$\endgroup\$
    – Mast
    Apr 15, 2020 at 11:39

1 Answer 1

3
\$\begingroup\$
  • The implementation of Gaussian(Random random) does not match the description. You should add the comment about the polar form. I know the equivalence. Some reviewers may not.

  • Polar form comes with the price. The loop

        do
        {
            x = Uniform(random, -1.0, 1.0);
            y = Uniform(random, -1.0, 1.0);
            r = x * x + y * y;
        }
        while (r >= 1 || r == 0);
    

    is potentially infinite. I'd be scared to see it in the production code.

  • The Wiki article warns against the price of the basic form coming from the computation of \$\sin\$ and \$\cos\$. I don't buy it. sincos is cheap. Thou shalt profile.

  • The condition r == 0 is pretty much guaranteed to fail, and I wouldn't trust the result of Math.Log(r) / r for really small r. Do not compare floating point values for equality. Chose an \$\epsilon\$, and compare for r < eps.

  • Maxwell calls MaxwellComponent three times, and each time MaxwellComponent computes the same MaxwellVariance. Seems excessive.

\$\endgroup\$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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