「エントロピーと秩序」のプログラムをpythonで書く～エントロピー

「エントロピーと秩序」のプログラムをpythonで書く第3弾。
「エントロピー」というタイトルがついているbasicのプログラムをpythonで書き直します。

前回100個の原子から成り立っている物質1(n1)と1500個の原子からなりたっている物質2(n2)にランダムに元気な原子を配置した時の、n1とn2の温度の揺らぎ可視化するプログラムをpython化しました。

今回のプログラムは、n1に配置されている100個の元気な原子がn2のほうへ一個ずつ移っていった時のエントロピーの変化を可視化するプログラムになっています。(モデル世界ということで、k=1で計算されています)
本家のプログラムにはグラフを表示する/しないを選択できる処理が入っていたのですが、グラフがないと面白みがないので省略しました。あと、各値をprintする処理も入っていたのですが、全体のエントロピーが最大となったポイントの値だけをグラフに表示するよう変更しています。

import math
import matplotlib.pyplot as plt

N = 100  # = int(input("ON状態にする原子数(1から100まで)"))

n2_values = []
s1_values = []
s2_values = []
s_values = []

max_s = 0
max_n2 = 0

for n2 in range(N):
    n1 = N - n2
    s1 = 0
    if n1 > 0:
        for i in range(0, n1):
            s1 += math.log((100 - i) / (n1 - i))
    s2 = 0
    if n2 > 0:
        for i in range(0, n2):
            s2 += math.log((1500 - i) / (n2 - i))
    s = s1 + s2
    if s > max_s:
        max_s = s
        max_n2 = n2
    n2_values.append(n2)
    s1_values.append(s1)
    s2_values.append(s2)
    s_values.append(s)

plt.plot(n2_values, s1_values, 'b', label="n1 entropy")
plt.plot(n2_values, s2_values, 'r', label="n2 entropy")
plt.plot(n2_values, s_values, 'y', label="total entropy")
plt.axvline(max_n2, label=f"max: {max_s:.1f} (n2={max_n2})", color="g")
plt.ylabel("entropy")
plt.xlabel("active atoms in n2")
plt.legend()
plt.show()

本家プログラムに

basic

210 FOR I=0 TO N1-1:
220 S1=LOG((100-I)/(N1-I))+S1
230 ...

という部分があり、basicのfor文の場合N1 = 1のときにI = 0で一回ループが回るのですが、pythonで

python

for i in range(n1 - 1):
    s1 += math.log((100 - i) / (n1 - i))
    ...

と書いてしまうと、n1 = 1の時にループが回らないので、どう書いたらよいのか少し悩みました。
このプログラムを実行すると、以下のようなグラフが表示されます。

グラフを見ると、n1から94個の元気な原子がn2に移った時に全体のエントロピーが最大となっています。
前回出てきた温度を計算する式を使って、元気な原子がn2に92～96個存在するときのn1とn2の温度を計算してみると、エントロピーが最大(n2=94)のときn1とn2の温度が一番近くなっていることがわかります。

>>> import math
>>> for i in range(92, 97):
...     n1 = 100 - i
...     n2 = i
...     t1 = 1 / math.log((100 - n1) / n1)
...     t2 = 1 / math.log((1500 - n2) / n2)
...     print(f"n2: {n2}, t1: {t1:.4f}, t2: {t2:.4f}, diff: {abs(t1 - t2):.4f}")
n2: 92, t1: 0.4094, t2: 0.3666, diff: 0.0429
n2: 93, t1: 0.3866, t2: 0.3681, diff: 0.0185
n2: 94, t1: 0.3634, t2: 0.3697, diff: 0.0062
n2: 95, t1: 0.3396, t2: 0.3712, diff: 0.0316
n2: 96, t1: 0.3147, t2: 0.3728, diff: 0.0581

「物事は自発的には全体のエントロピーが増大する方向にしか変化しない」
　 => 「部屋に熱いコーヒーを置いておくとやがて冷めて室温と同じになるが、それ以上に温度が下がったり、勝手に再び沸騰したりはしない」

ということがわかります。