プログラムを作って熱力学を理解する-温度、エネルギー、エントロピー-

初めに

今回はこの記事で使いましたエネルギーやエントロピーによるGの変化のアニメーションのコードの解説になります。

コード

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import os
 
T_list = np.arange(0, 1010, 10)
 
 
def H(x):
    y = 30 * np.sin(x)
    return y
 
 
def S(x):
    y = np.cos(x)
    return y
 
 
def G(x, y, T):
    z = x - T * y
    return z
 
 
H_list = np.arange(-10, 10, 1)
S_list = np.arange(0, 10, 1)
X_, Y_ = np.meshgrid(H_list, S_list)
H = H(H_list)
S = S(S_list)
X, Y = np.meshgrid(H, S)
Z = G(X, Y, T_list[0])
 
fig, ax = plt.subplots()
im = ax.contourf(X_, Y_, Z)
colorbar = plt.colorbar(im)
colorbar.set_label("G")
plt.xlabel("H")
plt.ylabel("S")
 
text = ax.text(
    0.05,
    0.95,
    "",
    transform=ax.transAxes,
    va="top",
    color="white",
    fontsize=15,
    fontweight="bold",
)
 
X0, Y0 = np.argmin(Z, axis=0)[0], np.argmin(Z, axis=1)[0]
X0, Y0 = np.random.choice(Z.shape[0], size=1), np.random.choice(Z.shape[1], size=1)
 
 
def Data(X0, Y0, T):
    Z = G(X[X0, Y0], Y[X0, Y0], T)
    return Z
 
 
if os.path.isfile("*.csv"):
    os.remove("*.csv")
Z0_min = Data(X0, Y0, T_list)
df = pd.DataFrame([T_list, Z0_min]).T
df.columns = [
    "T",
    "G_" + str(round(X_[X0, Y0][0], 3)) + "_" + str(round(Y_[X0, Y0][0], 3)),
]
df.to_csv("PointData.csv", index=None)
 
for T in T_list:
    df2 = pd.DataFrame(G(X, Y, T), index=S_list, columns=H_list)
    df2.to_csv("Data_" + str(T) + ".csv")
 
 
def animate(i):
    Z = G(X, Y, T_list[i])
    im = ax.contourf(X_, Y_, Z)
    text.set_text("T = {}".format(T_list[i]) + " K")
    return im.collections + [text]  
 
 
ani = animation.FuncAnimation(fig, animate, frames=len(T_list), interval=200, blit=True)
ani.save("G_animetion.gif", writer="pillow")
plt.show()

解説

まず、温度T、エネルギーH、エントロピーS、自由エネルギーGを定義します。

温度は0~1000Kまで10Kずつ定義しています。

T_list = np.arange(0, 1010, 10)

エネルギーは-10~10まで1ずつ変化させていますが、変化の具合をsin関数で変化させます。特に理由がないですが、説明がしやすい関数にしました(笑)。ここは適当な関数を当てはめていろいろ試してみてください。

def H(x):
    y = 30 * np.sin(x)
    return y

エントロピーは0~10まで1ずつ変化させていますが、変化の具合をcos関数で変化させます。こちらも特に理由がないです。ここは適当な関数を当てはめていろいろ試してみてください。

def S(x):
    y = np.cos(x)
    return y

Gは定義のままの関数となります。ここは変えないでHやSを変えて実行してみましょう。ただ、物質によって温度のべき乗を考えたりする場合は変えてみてもいいかもしれません。

def G(x, y, T):
    z = x - T * y
    return z

これらの定義ができれば、後はpythonの3Dプロットを作り処理を行います。ここの解説は別の記事でします。この処理が終われば、温度ごとにアニメーションを作ります。それが以下のコードです。

def animate(i):
    Z = G(X, Y, T_list[i])
    im = ax.contourf(X_, Y_, Z)
    text.set_text("T = {}".format(T_list[i]) + " K")
    return im.collections + [text] 
 
 
ani = animation.FuncAnimation(fig, animate, frames=len(T_list), interval=200, blit=True)

defは関数で各温度毎に3Dグラフを作成しています。後はパラパラ漫画のようにすべてのグラフをしたのがこの記事のアニメーションになります。

HやSの関数を色々変えてみてGが温度に対してどう変化するか試してみてください。Gに対して理解が深まるかと思います。もし、おもしろい関数があったら教えてください。

今回はこの記事のアニメーションにしたプログラムの解説をしました。物理を勉強して、身につけたプログラミングを活用すると現象をより理解することができます。どちらも平行して勉強できることが理想ですが、大変だと思うので、本ブログではそれをサポートできればと思います。

2023/9/5 J.A.

コメント

タイトルとURLをコピーしました