はじめに
「拡散」という現象は、熱力学・流体力学・生物物理など幅広い分野に登場する普遍的な物理プロセスである。本記事では、ブラウン運動（ランダムウォーク） に基づく粒子拡散シミュレーションを Python で実装し、初期に局在していた粒子集団が時間とともに空間的に均一化していく過程を定量的に可視化する。
特に、
- 粒子の空間分布の時間発展
- x 方向の確率密度関数（KDE推定）の収束過程
- 均一分布（理論値）への漸近的接近
を重点的に観察する。
(今回の記事は、統計力学についてあまり知らない私が、生成AIを使って拡散の様子を調べたかったために生成AIを使って作成した。)
理論的背景
拡散方程式
連続極限において、拡散現象は以下の拡散方程式（Fick の第二法則） で記述される。
$$
\frac{\partial c}{\partial t} = D \nabla^2 c
$$
ここで $c(\mathbf{r}, t)$ は粒子の濃度場、$D$ は拡散係数である。
一次元において初期条件が $c(x, 0) = \delta(x - x_0)$ のとき、解析解は：
$$
c(x, t) = \frac{1}{\sqrt{4\pi D t}} \exp\left(-\frac{(x - x_0)^2}{4Dt}\right)
$$
すなわち、初期にデルタ関数的に集中していた分布が時間とともにガウス分布に広がり、$t \to \infty$ で一様分布へ収束する。
ランダムウォークとの対応
粒子の微視的運動をランダムウォークでモデル化すると、各時間ステップ $\Delta t$ での変位は：
$$
\Delta x \sim \mathcal{N}(0,\ \sigma^2), \quad \sigma = \sqrt{2D\Delta t}
$$
に従う。このとき、アンサンブル平均としての分布関数は拡散方程式の解と一致する（アインシュタイン関係）。
シミュレーション設定
| パラメータ | 値 | 説明 |
|---|---|---|
| $N$ | 100 | 粒子数 |
| $D$ | 0.1 | 拡散係数 |
| $\Delta t$ | 0.01 | 時間ステップ |
| $\sigma$ | $\sqrt{2D\Delta t} \approx 0.045$ | 1 ステップあたりの標準偏差 |
| ステップ数 | 80 | シミュレーション総ステップ |
| 系のサイズ | $[-1, 1] \times [-1, 1]$ | 2次元矩形領域 |
| 境界条件 | 反射境界 | 壁での折り返し |
| 初期配置 | $-1 < x < 0$, $-1 < y < 1$ | 左半分に一様ランダム配置 |
実装
環境
Python 3.11
numpy
scipy
matplotlib
japanize_matplotlib
全ソースコード
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.stats import gaussian_kde
import japanize_matplotlib
# --- パラメータ ---
N = 100
DT = 0.01
D = 0.1
STEPS = 80
SNAPSHOT_STEPS = [1, 20, 40, 60, 80]
X_MIN, X_MAX = -1.0, 1.0
Y_MIN, Y_MAX = -1.0, 1.0
DPI = 100
FIG_W = 1280 / DPI
FIG_H = 1600 / DPI
rng = np.random.default_rng(42)
x = rng.uniform(X_MIN, 0.0, N)
y = rng.uniform(Y_MIN, Y_MAX, N)
x_grid = np.linspace(X_MIN, X_MAX, 300)
uniform_density_kde = 1 / (X_MAX - X_MIN)
def reflect(pos, lo, hi):
over_hi = pos > hi
over_lo = pos < lo
pos[over_hi] = 2 * hi - pos[over_hi]
pos[over_lo] = 2 * lo - pos[over_lo]
return pos
fig, (ax_scatter, ax_density) = plt.subplots(
2, 1, figsize=(FIG_W, FIG_H), dpi=DPI,
gridspec_kw={"height_ratios": [1, 1.2]}
)
fig.patch.set_facecolor("#1a1a2e")
for ax in (ax_scatter, ax_density):
ax.set_facecolor("#16213e")
ax.tick_params(colors="white", labelsize=18)
ax.xaxis.label.set_color("white")
ax.yaxis.label.set_color("white")
ax.title.set_color("white")
for spine in ax.spines.values():
spine.set_edgecolor("#444466")
ax_scatter.set_xlim(X_MIN, X_MAX)
ax_scatter.set_ylim(Y_MIN, Y_MAX)
ax_scatter.set_aspect("equal")
ax_scatter.set_xlabel("x", fontsize=20)
ax_scatter.set_ylabel("y", fontsize=20)
ax_scatter.set_title("粒子位置", fontsize=22, pad=14)
ax_scatter.axvline(0, color="#888899", lw=1.5, ls="--", alpha=0.5)
scatter = ax_scatter.scatter(x, y, s=40, c="#00d4ff", alpha=0.75)
step_text = ax_scatter.text(0.03, 0.97, "step = 0",
transform=ax_scatter.transAxes,
fontsize=18, va="top", color="white")
ax_density.set_xlim(X_MIN, X_MAX)
ax_density.set_ylim(0, 2.8)
ax_density.set_xlabel("x", fontsize=20)
ax_density.set_ylabel("確率密度", fontsize=20)
ax_density.set_title("x 方向の密度分布の推移", fontsize=22, pad=14)
ax_density.axhline(uniform_density_kde, color="white", lw=2,
ls="--", alpha=0.6, label="均一分布 (理論値)")
ax_density.axvline(0, color="#888899", lw=1.5, ls="--", alpha=0.4)
snapshot_lines = []
snapshot_colors = plt.cm.cool(np.linspace(0.2, 0.9, len(SNAPSHOT_STEPS)))
kde_init = gaussian_kde(x, bw_method=0.25)
(kde_line,) = ax_density.plot(x_grid, kde_init(x_grid),
color="#ff6b6b", lw=3, label="現在の密度", zorder=5)
ax_density.legend(loc="upper right", fontsize=14,
facecolor="#16213e", labelcolor="white",
edgecolor="#444466")
plt.tight_layout(pad=2.5)
def update(frame):
global x, y
sigma = np.sqrt(2 * D * DT)
x += rng.normal(0, sigma, N)
y += rng.normal(0, sigma, N)
x = reflect(x, X_MIN, X_MAX)
y = reflect(y, Y_MIN, Y_MAX)
scatter.set_offsets(np.c_[x, y])
step_text.set_text(f"step = {frame + 1}")
kde = gaussian_kde(x, bw_method=0.25)
kde_line.set_ydata(kde(x_grid))
actual_step = frame + 1
if actual_step in SNAPSHOT_STEPS:
idx = SNAPSHOT_STEPS.index(actual_step)
color = snapshot_colors[idx]
ax_density.plot(x_grid, kde(x_grid),
color=color, lw=1.8, alpha=0.55,
label=f"step={actual_step}", zorder=3)
snapshot_lines.append(True)
ax_density.legend(loc="upper right", fontsize=13,
facecolor="#16213e", labelcolor="white",
edgecolor="#444466")
return scatter, step_text, kde_line
ani = animation.FuncAnimation(
fig, update, frames=STEPS, interval=40, blit=False
)
ani.save("diffusion_animation.gif", writer="pillow", fps=25,
savefig_kwargs={"facecolor": "#1a1a2e"})
plt.show()
実装のポイント解説
反射境界条件
壁を透過させず折り返す境界条件は、系の粒子数を保存しながら閉じた領域内での拡散を正確に再現する。
def reflect(pos, lo, hi):
over_hi = pos > hi
over_lo = pos < lo
pos[over_hi] = 2 * hi - pos[over_hi]
pos[over_lo] = 2 * lo - pos[over_lo]
return pos
壁を超えた分を折り返すことで、壁面での「完全弾性衝突」を模倣している。
カーネル密度推定（KDE）
ヒストグラムの代わりに KDE（Kernel Density Estimation） を用いることで、滑らかな確率密度関数を推定できる。
kde = gaussian_kde(x, bw_method=0.25)
density = kde(x_grid)
バンド幅 bw_method=0.25 は Scott の法則より小さい値を手動設定しており、N=100 という少数粒子においても過度な平滑化を抑制している。
スナップショット描画
SNAPSHOT_STEPS で指定したタイミングの密度曲線をアクスに残すことで、分布の時間発展を一枚のパネルで比較できる。
結果と考察
密度分布の時間発展
シミュレーション結果から、以下の傾向が確認された。
| ステップ | 観察された密度プロファイル |
|---|---|
| step = 1 | $-1 < x < 0$ に強くピーク |
| step = 20 | ピークが右方向へ拡散し始める |
| step = 40 | $x=0$ をまたいで右半分にも広がる |
| step = 60 | 全域でほぼ均一に近づく |
| step = 80 | 理論値（均一密度 = 0.5）に漸近 |
この収束速度は拡散係数 $D$ と系のサイズ $L$ に依存し、特性時間スケールは：
$$
\tau \sim \frac{L^2}{D} = \frac{4}{0.1} = 40 \text{ steps} \times \Delta t
$$
と見積もられ、シミュレーション結果と定性的に一致する。
粒子数による統計的揺らぎ
N=100 という有限粒子数のため、理論値への収束には確率的揺らぎが伴う。この揺らぎの大きさは中心極限定理より $\sim 1/\sqrt{N}$ のオーダーであり、N を増やすほど理論曲線へのフィットが改善される。
まとめ
本記事では、2次元矩形領域内でのランダムウォークによる粒子拡散を Python で実装・可視化した。
- 拡散方程式とランダムウォークの等価性 をシミュレーションで直接確認できた
- KDE によるリアルタイム密度推定で、均一分布への収束過程を定量的に追跡できた
- GIF アニメーションとして出力することで、動的な拡散過程を直感的に伝えられた
今後の発展として、以下が考えられる。
- 障害物や濃度勾配の導入
- 粒子間相互作用（排除体積効果など）の追加
- 3次元への拡張
- KL ダイバージェンスによる均一分布への収束の定量評価
参考文献
- Einstein, A. (1905). "Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen." Annalen der Physik, 322(8), 549–560.
- Chandrasekhar, S. (1943). "Stochastic problems in physics and astronomy." Reviews of Modern Physics, 15(1), 1–89.
- Risken, H. (1989). The Fokker-Planck Equation. Springer.