Как я могу реализовать дистрибутив Максвелла? - PullRequest
Новая
       23

Как я могу реализовать дистрибутив Максвелла?

1 голос
/ 15 апреля 2020

Мне дано решение следующей задачи (Этот текст переведен с русского. Итак, могут возникнуть некоторые проблемы с переводом):

... Еще один способ извлечь из нормальное распределение состоит в том, чтобы извлечь два независимых случайных числа из равномерного распределения x1, x2 ∈ [0: 0, 1: 0). Затем примените следующее преобразование:
enter image description here
, в результате чего получаются два случайно выбранных числа n 1 , n 2 из нормального распределения с нулевым ожидаемым значением и единицей дисперсии.

Чтобы изменить параметры распределения на другие параметры, например, на ожидаемое значение и дисперсию, вы должны умножить результат нарисуйте и добавьте, то есть
enter image description here
В приведенном выше уравнении N (μ, σ) является случайной величиной с нормальным распределением с ожидаемым значением μ и дисперсия σ.

Согласно распределению Максвелла, каждый компонент (x, y или z) вектора скорости v является случайной величиной из нормального распределения с ожидаемое нулевое значение и дисперсия
enter image description here
где m - масса молекулы, T - температура в градусах Кельвина, k B - постоянная Больцмана.

Ваша задача: нарисовать 10000 вело Векторы ty для молекулы азота N 2 при 300К. Рассчитайте среднюю длину этих векторов и, следовательно, среднее значение скорости молекулы азота, используя формулу:
enter image description here 

public class Maxwell
{
    public double N1 { get; private set; }
    public double N2 { get; private set; }
    public void Compute(Random random)
    {
        double x1 = random.NextDouble();
        double x2 = random.NextDouble();

        N1 = Math.Sqrt(-2 * Math.Log(x1)) * Math.Cos(2 * Math.PI * x2);
        N2 = Math.Sqrt(-2 * Math.Log(x1)) * Math.Sin(2 * Math.PI * x2);
    }
}

public class Program
{
    static void Main(string[] args)
    {
        Random r = new Random();
        Maxwell m = new Maxwell();
        m.Compute(r);

        double n1 = m.N1;
        double n2 = m.N2;

        //.....? 
    }
}

Я не понимаю, как реализовать N (μ, σ) из n1 и n2 и как оттуда добраться до вектора v .

Может кто-нибудь помочь?

Редактировать: Я реализовал это на основе Эри c Ответ Липперта : 

using System;

public class CommonDistributions
{
    public static double Uniform(Random random)
    {
        return random.NextDouble();
    }

    static double Gaussian(Random random)
    {
        return Math.Sqrt(-2 * Math.Log(Uniform(random))) * Math.Cos(2 * Math.PI * Uniform(random));
    }
    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 = {sum/N}");

        string str = string.Empty;
    }
}

Я не уверен в значениях температуры и масса азота 2.

Было бы неплохо, если кто-то может прокомментировать код.

1 Ответ

4 голосов
/ 15 апреля 2020

В этой ситуации нужно подумать: «Предположим, у меня есть волшебный ящик c, который ответил на конкретный вопрос c, который я ему задал; каковы будут входы и выходы этого ящика?» и написать метод, который реализует этот блок .

Начните с самого простого блока. Нет входных данных, выходные данные - это равномерно распределенное число от нуля до единицы: 

static Random random = new Random();
static double Uniform() => random.NextDouble();

ОК, теперь у нас есть новый инструмент в нашем наборе инструментов. Какая у нас следующая волшебная c коробка? Нет входных данных, выходной это нормально распределенное число со средним нулем и стандартным отклонением: 

static double StandardNormal() =>
  Sqrt(-2 * Log(Uniform())) * Cos(2 * PI * Uniform());

И у нас есть другой инструмент. Что мы можем построить с этим? Входы: среднее и стандартное отклонение, выход, нормально распределенное число с этим средним и стандартным отклонением: 

static double Normal(double mean, double sigma) =>
  sigma * StandardNormal() + mean;

Хорошо, теперь что нам нужно? Дисперсия как функция массы и температуры: 

static double KB = 1.38064852e-23;
static double MaxwellVariance(double mass, double temperature) => 
  Sqrt(KB * temperature / mass);

Супер, мы движемся прямо. Теперь что нам нужно? Входные данные - это масса и температура, выходные данные - это одна случайная максвелловская составляющая скорости: 

static double MaxwellComponent(double mass, double temperature) =>
  Normal(0.0, MaxwellVariance(mass, temperature));

Что нам теперь нужно? Тип для представления вектора: 

struct Vector
{
  public double X { get; }
  public double Y { get; }
  public double Z { get; }
  public Vector(double x, double y, double z)
  {
    this.X = x;
    this.Y = y;
    thiz.Z = z;
  }
}

Что нам нужно дальше? Случайный вектор: 

static Vector MaxwellVector(double mass, double temperature) =>
  ...

Вы можете взять его отсюда? Что вам понадобится дальше? Снова продолжайте разбивать его на однострочники . Не увлекайся Нет премии за написание длинного кода, который вы не понимаете.

Техника здесь разделяй и властвуй . С этими проблемами вы почти всегда можете написать метод длиной менее пяти строк кода, который вычисляет только одну вещь . Так сделай это; вычислять только одну вещь каждый раз, и тогда у вас есть новый инструмент в вашем наборе инструментов для вычисления вещи next . Более того, у вас есть набор методов, каждый из которых (1), очевидно, корректен, потому что это всего лишь одна строка кода, и (2) тестируемый! Напишите набор тестов!

ОБНОВЛЕНИЕ: Вопрос был обновлен для реализации некоторых из этих идей, и выглядит он довольно неплохо. Есть дополнительный вопрос о температуре и массе.

Температура выглядит хорошо; 300K. Но масса совершенно не права. В инструкциях говорится, что нужно использовать массу одну молекулу , но вы ввели массу один моль молекул .

Помните, "моль" - это как "пара" или "дюжина". Пара - это две вещи, дюжина - двенадцать, моль - около 600000000000000000000000 вещей. Очевидно, молекула N2 не весит 28 граммов. Скорее 600000000000000000000000 молекул N2 весит 28 граммов.

Помните также, что метрические c единицы массы и объема были выбраны совершенно произвольно. Если вы возьмете окружность Земли, разделите ее на 4 миллиарда, сделайте коробку с кубиками c со сторонами этой длины и заполните ее водой - это масса одного грамма.

Мы выбрали значение, связанное с «моль», потому что оно обладает свойством, что моль молекулы того же типа имеет массу, равную массе атома c молекулы в граммах. Итак, в восемнадцати из этих маленьких коробочек есть один моль молекул воды. Использование молярной массы - это просто удобство, потому что это делает числа более «разумными» для наших целей; обычно мы привыкли думать о некотором количестве грамм воды, а не о каком-то количестве молекул воды, но ваша проблема касается только десяти тысяч молекул, а не десяти тысяч граммов. Итак, что вы хотите сделать, это разделить массу одного моля на количество вещей в моле, и это даст вам массу одной молекулы в граммах .

Следующая вещь сделать это сделать анализ единицы, чтобы определить, должна ли масса быть в граммах или килограммах! Для КБ мы имеем значение 1.38E-23, которое, как отмечает Википедия, имеет единицы Джоулей на Кельвин. Как мы это используем? Мы берем квадрат root КБ * Т / М. Какими единицами должен быть квадрат root? Это стандартное отклонение скорости , которое имеет единицы метров в секунду, поэтому нам нужно KB*T/M, чтобы иметь квадраты метров в секунду в квадрате.

  • КБ - Джоуль на Кельвин ; Это Кельвин, поэтому KB * T имеет единицы в Джоулях.
  • Джоули имеют единицы: килограмм метров в квадрате в секунду в квадрате.
  • Поэтому, чтобы получить квадраты в метрах в секунду в квадрате нам нужно разделить на килограмм , а не грамм .

Так что вам нужно, граммы на моль, деленные на молекулы на моль, чтобы получить граммы на молекулу , а затем преобразовать это в килограммы на молекулу.

Имеет смысл? Привыкнуть проводить анализ единиц для каждой проблемы . Это уловило так много моих ошибок, когда я был студентом-физиком еще в темные века.

В сторону: Говоря об анализе единиц, кое-что, что следует быть осторожным: отрывок из вашего текста называет стандартное отклонение дисперсия , но стандартное отклонение фактически определяется как квадрат root дисперсии. Такое использование чрезвычайно распространено, и вы должны из контекста определить, означает ли «дисперсия» «действительно дисперсию» или, в данном случае, стандартное отклонение.

То есть текст должен сказать "N (μ, σ) - случайная величина с нормальным распределением с ожидаемым значением μ и дисперсией σ 2 ". Или следует сказать, что «N (μ, σ) является случайной величиной с нормальным распределением с ожидаемым значением μ и стандартным отклонением σ». Следите за этим и читайте с защитой.

Еще одно замечание: вы могли заметить, что способ представления дистрибутивов очень "неуклюжий". Такое чувство, что вам нужно проделать большую работу, чтобы представить что-то довольно простое. Мое текущее исследование в probabilisti c языках , которые делают этот вид работы очень простым. На вероятностном языке c мы бы представили ваш рабочий процесс примерно так: 

IDistribution<double> Speed(double mass, double temp)
{
  IDistribution<double> c = 
    Normal.Distribution(0.0, MaxwellVariance(mass, temp))
  double x = sample c;
  double y = sample c;
  double z = sample c;
  return Sqrt(x*x + y*y + z*z);
}
...
double mean = Speed(mass, temp).Mean(10000);

(Если это похоже на асинхронный метод c с заменой Task<T> на IDistribution<T> и заменой await на sample, это потому, что это так; и асинхронные, и вероятностные c рабочие процессы могут быть реализованы как сопрограммы.)

Если вас интересует тема вероятностных c языков, у меня есть небольшое, но продолжительное введение который начинается здесь: https://ericlippert.com/2019/01/31/fixing-random-part-1/

...