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

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

温度とは、超簡単に言ってしまうと「元気な原子が存在する割合」を示す値です。
実際の世界には超元気な原子やまあまあ元気な原子、全然元気のない原子といろいろな元気度の原子がいるのですが、このプログラムは「元気な原子(ON状態の原子)」と「全然元気のない原子(OFF状態の原子)」で構成された世界をモデリングしています。

100個の原子から成り立っている物質1(左上の部分)と1500個の原子から成り立っている物質2(それ以外の部分)が接触して平衡状態にあるとき、元気な分子は物質1と物質2にランダムに配置された状態になります。その瞬間瞬間を可視化すると、原子数の少ない物質1(左上の物質)のほうが温度に揺らぎがありますが、同じくらいの温度になっていることがわかります。

なるべく元のプログラムに忠実に書き直したかったのですが、大量に出てくるgoto分をforとかifに書き直したり、可視化する部分をmatplotlibに置き換えたりしているうちに結構変わっちゃいました。

import time
import math
import random
import numpy as np
import matplotlib.pyplot as plt

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

fig = plt.figure()
ax1 = fig.add_subplot(1, 2, 1)
ax1.set_xlim(-1, NUMBER_OF_EXCUTIONS)
ax1.set_ylim(-1, 1)
ax2 = fig.add_subplot(1, 2, 2) 
ax2.plot((-1, 10), (30, 30), (10, 10), (30, 41), color = 'grey')
ax2.set_xlim(-1, 41)
ax2.set_ylim(-1, 41)
sc1 = ax1.scatter([], [], c="b")
sc2 = ax2.scatter([], [], c="r")

tx_coordinates = []
ty_coordinates = []
colors = []

for i in range(NUMBER_OF_EXCUTIONS):
    n1 = 0
    n2 = 0
    n = 0
    f = np.zeros([40, 40])
    x_coordinates = []
    y_coordinates = []
    c = []
    while(n < N):
        x = random.randint(0, 39)
        y = random.randint(0, 39)
        if f[x, y] == 1:
            continue
        f[x, y] = 1
        x_coordinates.append(x)
        y_coordinates.append(y)
        if x < 10 and y >= 30:
            c.append("b")
            n1 += 1
        else:
            c.append("r")
            n2 += 1
        n += 1
    if n1 == 0:
        t1 = 0
    else:
        t1 = 1 / math.log((100 - n1) / n1)
    if n2 == 0:
        t2 = 0
    else:
        t2 = 1 / math.log((1500 - n2) / n2)
    sc2.set_offsets(np.c_[x_coordinates, y_coordinates])
    sc2.set_color(c)
    tx_coordinates.append(i)
    ty_coordinates.append(t1)
    colors.append("b")
    tx_coordinates.append(i)
    ty_coordinates.append(t2)    
    colors.append("r")
    sc1.set_offsets(np.c_[tx_coordinates, ty_coordinates])
    sc1.set_color(colors)
    plt.pause(0.1)
plt.waitforbuttonpress()

左側のグラフが温度の分布を示し、右側のグラフが元気な分子の配置を示しています。