"""
配置の数 W(ni) の分布と熱力学的極限への収束を可視化するスクリプト。
概要:
総要素数 $N$ のうち、特定の状態に属する要素数が $ni$ である場合の組み合わせ（配置の数） $W(ni)$ を計算し、
システムサイズ $N$ による分布の変化をプロットします。
詳細説明:
- 配置の数 $W = {}_N C_{ni}$ を計算します。
- 大きな $N$ を扱うため、スターリングの近似 $\ln N! \approx N \ln N - N$ を用いて計算を行います。
- 各 $N$ において最大値（$ni = N/2$）で正規化し、相対的な確率分布として表示します。
- $N$ が増大するにつれて、$ni/N = 0.5$ 付近の分布がいかに鋭くなるか（相対的な揺らぎの減少）を確認できます。
- これは、多粒子系においてマクロな物理量が一点に定まる「熱力学的極限」を数学的に示唆します。
理論式:
$$\ln W(ni) \approx N \ln N - ni \ln ni - (N-ni) \ln (N-ni)$$
関連リンク: :doc:`plotW_usage`
"""
import matplotlib.pyplot as plt
import numpy as np
from numpy import log, exp, sqrt
[ドキュメント]
def calculate_normalized_W(N, ni_array):
"""
総要素数 N における特定の状態の要素数 ni に対する組み合わせの数 W を計算します。
最大値で正規化した値を返します。
:param N: float: 総要素数。
:param ni_array: np.ndarray: 評価する要素数 ni の配列。
:returns: list: 正規化された配置の数 W/Wmax。
"""
# スターリングの近似を用いた ln(W_max) の計算 (ni = N/2 のとき)
lnW_max = N * log(N) - 2.0 * (N / 2.0) * log(N / 2.0)
W_list = []
for ni in ni_array:
# ln(ni!) の項
ni_term = ni * log(ni) if ni > 0 else 0
# ln((N-ni)!) の項
Nni_term = (N - ni) * log(N - ni) if (N - ni) > 0 else 0
# Wmaxで規格化した ln(W/Wmax)
lnW_rel = N * log(N) - ni_term - Nni_term - lnW_max
W_list.append(exp(lnW_rel))
return W_list
[ドキュメント]
def main():
"""
異なる N に対して W(ni) を計算し、グラフを生成・保存します。
"""
# Nの値のリスト (5から始まり、倍々で増やす)
N_values = [5.0 * 2**i for i in range(0, 32, 4)]
fig, ax = plt.subplots(figsize=(10, 7))
for N in N_values:
# 分布の中心付近を重点的に計算するための範囲設定
N_avg = N * 0.5
Nrange = 6.0 * sqrt(N)
# 数値的に安定した範囲で ni を生成
ni_min = max([1.0, N_avg - Nrange])
ni_max = N_avg + Nrange
log_ni = np.linspace(log(ni_min), log(ni_max), 500)
ni_array = exp(log_ni)
# 配置の数の計算
W_vals = calculate_normalized_W(N, ni_array)
x_vals = ni_array / N
ax.plot(x_vals, W_vals, label=f'N = {N:.1e}', linewidth=1.5)
# グラフの装飾
ax.set_title('Normalized Number of Configurations ($W/W_{max}$)', fontsize=16)
ax.set_xlabel('$n_i / N$ (Occupancy Ratio)', fontsize=14)
ax.set_ylabel('$W(n_i) / W_{max}$', fontsize=14)
# 物理的に意味のある範囲（0.5付近の収束を見るため）
ax.set_xlim([0.1, 1.0])
ax.set_xscale('log')
ax.legend(title='System Size N', fontsize=10)
ax.grid(True, which='both', linestyle='--', alpha=0.5)
plt.tight_layout()
# 画像の保存
save_path = 'plotW.png'
plt.savefig(save_path, dpi=300)
print(f"Graph saved as {save_path}")
plt.show(block=False)
input("Press ENTER to terminate>>")
if __name__ == "__main__":
main()