Box-Muller法でマクスウェル分布に従う乱数生成

---

Box-Muller法の原理

互いに独立な $[0:1]$ の乱数 $a, b$ を用いて、平均 $\mu$、標準偏差 $\sigma$ の正規分布に従う乱数を次のように表すことができる (wiki様)。 $$Z_1 = \mu + \sigma \sqrt{-2 \ln a} \cos 2\pi b$$ $$Z_2 = \mu + \sigma \sqrt{-2 \ln a} \sin 2\pi b$$

$\ln a$ は $a=0$ で発散するので $1-a$ にしてもいい (0≦a＜1 に設定されている場合が多い)。

---

マクスウェル分布に従う乱数生成

3D空間の各次元 ($v_x, v_y, v_z$) のマクスウェル分布 $$f(v_x) = \sqrt{\frac{m}{2 \pi k T}} \exp{\left( -\frac{mv_x^2}{2kT} \right)}$$

は $\mu = 0$、$\sigma = \sqrt{kT/m}$ の正規分布なので、

$$v_x = \sqrt{- \frac{2kT}{m} \ln a} \cos 2\pi b$$ $$v_y = \sqrt{- \frac{2kT}{m} \ln a} \sin 2\pi b$$

で乱数を生成することができる。$v_z$ は他の一様乱数 $c,d$を用意して

$$v_z = \sqrt{- \frac{2kT}{m} \ln c} \cos 2\pi d$$

で用意する。

次に、 $$v = \sqrt{v_x^2 + v_y^2 + v_z^2}$$ として求めた速さは $$f(v) = \sqrt{\frac{m}{2 \pi kT}}^3 v^2\exp{\left( -\frac{mv^2}{2kT} \right)} $$ に従う。

---

数値計算

• $v_x, v_y, v_z$ を生成
• $v = \sqrt{v_x^2 + v_y^2 + v_z^2}$ を計算する
• プロット

program rnd_maxerll
implicit none

integer, parameter :: N = 300000
real(8), parameter :: k = 1.38d-23
real(8), parameter :: m = 1.67d-27
real(8), parameter :: pi = 4.d0*atan(1.d0)
real(8), parameter :: T = 500.d0
real(8), parameter :: sigma = sqrt( k*T/m ), mu=0.d0
real(8) :: a, b, c, d, vx, vy, vz, v
integer :: i, j, l


do i = 1, N
    call random_number(a)
    call random_number(b)
    call random_number(c)
    call random_number(d)

    vx = sigma * sqrt( -2.d0 * log(a) ) * cos(2.d0*pi*b)
    vy = sigma * sqrt( -2.d0 * log(a) ) * sin(2.d0*pi*b)
    vz = sigma * sqrt( -2.d0 * log(c) ) * cos(2.d0*pi*d)
    v  = sqrt( vx**2 + vy**2 + vz**2 )
    write(55,*) vx, vy, vz, v
end do

import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('fort.55')
fig, axes = plt.subplots(2, 1, figsize=(6, 6))

axes[0].hist(data[:, 0], bins=30, edgecolor='black')
axes[0].set_title(r'1D velocity $v_x$')
axes[0].set_xlabel(r'$v_x (m/sec)$')
axes[0].set_ylabel('N of Data')

axes[1].hist(data[:, 3], bins=30, edgecolor='black')
axes[1].set_title(r'3D velocity $v = \sqrt{v_x^2 + v_y^2 + v_z^2}$')
axes[1].set_xlabel(r'$v (m/sec)$')
axes[1].set_ylabel('N of Data')

plt.tight_layout()
plt.savefig('hist.png', format='png', dpi=300)

1枚目(上)は$v_x$の分布、1枚目(下)は$v$の分布。
2枚目は規格化して分布関数を構成し、理論式と比べた図。縦線は熱速度 (ピーク速度) $v_{th} = \sqrt{\frac{2kT}{m}}$ を表している。
両者、よく一致していることがわかる。