############################################################################
# 箱の中のN個の自由粒子（古典力学）
############################################################################
import os
import numpy as np
import random
import matplotlib.pyplot as plt
import matplotlib.animation as animation
random.seed(0)
#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 16 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
#################################################################
## 物理系の設定
#################################################################
#質量
m_min = 1.0
m_max = 10.0
#粒子の半径
radius_min = 0.002
radius_max = 0.02
#粒子数
NM = 100
#箱の長さ
L = 1.0
#粒子に加わる力
fex = np.array([0.0, 0.0, 0.0])
#初期位置と速度
v_max = 2.0
#粒子に加わる力
def F (self, t, r_t, v_t ):
return self.fex + self.fc - self.m * self.beta * v_t
#時間区間
t_min = 0
t_max = 5.0
#時間区間数
NT = 500
skip = 1
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#
lw = 0.01
######################################
# ４次のルンゲ・クッタ クラス
######################################
class RK4:
#コンストラクタ
def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
self.m = m
#現時刻の位置と速度
self.r = r0.copy()
self.v = v0.copy()
#１ステップ前の時刻の位置と速度
self._r = r0.copy()
self._v = v0.copy()
#差分
self.dr = np.array([0.0, 0.0, 0.0])
self.dv = np.array([0.0, 0.0, 0.0])
#時間間隔
self.dt = dt
#外力
self.fex = np.array([0.0, 0.0, 0.0])
#衝突力
self.fc = np.array([0.0, 0.0, 0.0])
#粘性係数
self.beta = 0
#半径
self.radius = 0
#内部
self.__v1 = np.array([0.0, 0.0, 0.0])
self.__v2 = np.array([0.0, 0.0, 0.0])
self.__v3 = np.array([0.0, 0.0, 0.0])
self.__v4 = np.array([0.0, 0.0, 0.0])
self.__a1 = np.array([0.0, 0.0, 0.0])
self.__a2 = np.array([0.0, 0.0, 0.0])
self.__a3 = np.array([0.0, 0.0, 0.0])
self.__a4 = np.array([0.0, 0.0, 0.0])
#速度を与えるメソッド
def V(self, t, r, v):
return v.copy()
##########################################################
#加速度を与えるメソッド
def A(self, t, r_t, v_t):
return F(self, t, r_t, v_t ) / self.m
##########################################################
#時間発展を計算するメソッド
def timeEvolution(self, t):
#１段目
self.__v1 = self.V(
t,
self.r,
self.v
)
self.__a1 = self.A(
t,
self.r,
self.v
)
#２段目
self.__v2 = self.V(
t + self.dt / 2.0,
self.r + self.__v1 * self.dt / 2.0,
self.v + self.__a1 * self.dt / 2.0
)
self.__a2 = self.A(
t + self.dt / 2.0,
self.r + self.__v1 * self.dt / 2.0,
self.v + self.__a1 * self.dt / 2.0
)
#３段目
self.__v3 = self.V(
t + self.dt / 2.0,
self.r + self.__v2 * self.dt / 2.0,
self.v + self.__a2 * self.dt / 2.0
)
self.__a3 = self.A(
t + self.dt / 2.0,
self.r + self.__v2 * self.dt / 2.0,
self.v + self.__a2 * self.dt / 2.0
)
#４段目
self.__v4 = self.V(
t + self.dt,
self.r + self.__v3 * self.dt,
self.v + self.__a3 * self.dt
)
self.__a4 = self.A(
t + self.dt,
self.r + self.__v3 * self.dt,
self.v + self.__a3 * self.dt
)
#差分の計算
self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0
self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0
#################################################################
## 時間発展の計算
#################################################################
#壁との衝突判定（戻値：法線ベクトル）
def CollisionWallDetection( particle ):
#衝突リスト
collisionlist = []
if(particle.r[0] + particle.radius > L):
collisionlist.append({ "type": 1, "normal" : np.array([-1.0, 0.0, 0.0]) })
if(particle.r[0] - particle.radius < 0):
collisionlist.append({ "type": 1, "normal" : np.array([1.0, 0.0, 0.0]) })
if(particle.r[1] + particle.radius > L):
collisionlist.append({ "type": 1, "normal" : np.array([0.0, -1.0, 0.0]) })
if(particle.r[1] - particle.radius < 0):
collisionlist.append({ "type": 1, "normal" : np.array([0.0, 1.0, 0.0]) })
if(particle.r[2] + particle.radius > L):
collisionlist.append({ "type": 1, "normal" : np.array([0.0, 0.0, -1.0]) })
if(particle.r[2] - particle.radius < 0):
collisionlist.append({ "type": 1, "normal" : np.array([0.0, 0.0, 1.0]) })
return collisionlist
#粒子の運動エネルギー
def T():
_T = 0
for nm in range(NM):
_T += 1.0/2.0 * particles[nm].m * np.sum( particles[nm].v**2 )
return _T
#粒子の位置エネルギー
def U():
_R = 0
for nm in range(NM):
_R += - np.dot(particles[nm].fex, particles[nm].r)
return _R
#系の時間発展の計算
def timeEvolution( t, timeScale ):
#複数衝突
cylinderCollisionList = []
for nm in range(NM):
#時間スケールの指定
particles[nm].dt = dt / timeScale
#ルンゲ・クッタ法による時間発展
particles[nm].fc = np.array([0.0, 0.0, 0.0])
particles[nm].timeEvolution( t )
#位置と速度を更新
particles[nm].r += particles[nm].dr
particles[nm].v += particles[nm].dv
#壁との衝突リスト
CollisionWallList = CollisionWallDetection( particles[nm] )
#衝突あり
if( len(CollisionWallList) > 0 ):
#２つ以上の壁との衝突がある場合 → 高精度計算モードへ
if( len(CollisionWallList) > 1 ):
return False
#以下、衝突が１箇所
normal = CollisionWallList[ 0 ]["normal"]
#いったん巻き戻しておく
particles[nm].r -= particles[nm].dr
particles[nm].v -= particles[nm].dv
#衝突力
if( CollisionWallList[ 0 ]["type"] == 1 ): #壁との衝突
#print(nm)
#相対速度
barV = - particles[nm].v
particles[nm].fc = (2.0 * particles[nm].m * np.dot(barV, normal) * normal ) / particles[nm].dt - np.dot(particles[nm].fex, normal) * normal
#位置と速度を更新
particles[nm].timeEvolution( t )
particles[nm].r += particles[nm].dr
particles[nm].v += particles[nm].dv
#粒子同士の衝突判定
for nm1 in range(NM):
for nm2 in range(NM):
if( nm1 >= nm2 ): continue
#２点間ベクトル
r12 = particles[nm2].r - particles[nm1].r
r12_length = np.sqrt(np.sum( r12**2 ))
if( r12_length < particles[nm1].radius + particles[nm2].radius ):
#print( t, nm1, nm2 )
#すでに衝突力が与えられている場合 → 高精度計算モードへ
if( (particles[nm1].fc[0] != 0 or particles[nm1].fc[1] != 0 or particles[nm1].fc[2] != 0)
or ( particles[nm2].fc[0] != 0 or particles[nm2].fc[1] != 0 or particles[nm2].fc[2] != 0 )):
#時間発展：失敗
return False
n12 = r12 / r12_length
v12 = particles[nm2].v - particles[nm1].v
m1 = particles[nm1].m
m2 = particles[nm2].m
barf = particles[nm2].fex / m2 - particles[nm1].fex / m1
f21 = 2.0 * m1 * m2 / ( m1 + m2 ) * np.dot(v12, n12) * n12 / particles[nm1].dt
f21 += m1 * m2 / ( m1 + m2 ) * np.dot(barf, n12) * n12
f12 = -f21
#衝突力の設定
particles[nm1].fc = f21
particles[nm2].fc = f12
#いったん巻き戻しておく
particles[nm1].r -= particles[nm1].dr
particles[nm1].v -= particles[nm1].dv
particles[nm2].r -= particles[nm2].dr
particles[nm2].v -= particles[nm2].dv
#時間発展を計算
particles[nm1].timeEvolution( t )
particles[nm2].timeEvolution( t )
#位置と速度を更新
particles[nm1].r += particles[nm1].dr
particles[nm1].v += particles[nm1].dv
particles[nm2].r += particles[nm2].dr
particles[nm2].v += particles[nm2].dv
#２点間ベクトル
r12 = particles[nm2].r - particles[nm1].r
r12_length = np.sqrt(np.sum( r12**2 ))
#計算後の２粒子が重なっている場合は失敗
if( r12_length < particles[nm1].radius + particles[nm2].radius ):
#時間発展：失敗
return False
#計算後、壁とあるいは粒子同士が重なっている場合は失敗
for nm1 in range(NM):
#壁との衝突リスト
CollisionWallList = CollisionWallDetection( particles[nm1] )
#衝突あり
if( len(CollisionWallList) > 0 ):
#時間発展：失敗
return False
for nm2 in range(NM):
if( nm1 >= nm2 ): continue
#２点間ベクトル
r12 = particles[nm2].r - particles[nm1].r
r12_length = np.sqrt(np.sum( r12**2 ))
if( r12_length < particles[nm1].radius + particles[nm2].radius ):
#時間発展：失敗
return False
#時間発展：正常
return True
#高精度計算モード
def highPrecisionMode(t, level):
print ( "レベル" , level, "t=", round(t, level+3 ))
timeScale = 10**level
for nm in range(NM):
particles[nm].r = particles[nm]._r.copy()
particles[nm].v = particles[nm]._v.copy()
for _n in range(0, 10):
t_level = t + dt / timeScale * _n
result_level = timeEvolution( t_level, timeScale )
if( result_level == False ):
#もっと高精度が必要
level += 1
highPrecisionMode(t_level, level)
level -= 1
else:
#時間発展前の位置と速度を保持
for nm in range(NM):
particles[nm]._r = particles[nm].r.copy()
particles[nm]._v = particles[nm].v.copy()
#################################################################
## 計算開始
#################################################################
#粒子
particles = []
for nm in range(NM):
radius = radius_min + ( radius_max - radius_min) * random.random()
m = m_min + ( m_max - m_min) * random.random()
theta_v = 2.0 * np.pi * random.random()
phi_v = 0 * 2.0 * np.pi * random.random()
_v_max = v_max * random.random()
v0 = np.array([
_v_max * np.sin( theta_v ) * np.cos( phi_v ),
_v_max * np.sin( theta_v ) * np.sin( phi_v ),
_v_max * np.cos( theta_v )
])
for i in range(1000):
r0 = np.array([
radius + (L - 2.0 * radius) * random.random(),
L / 2.0,
radius + (L - 2.0 * radius) * random.random()
])
flag = True
for p in particles:
#２点間ベクトル
r12 = r0 - p.r
r12_length = np.sqrt(np.sum( r12**2 ))
if( r12_length < radius + p.radius ):
flag = False
break
if(flag == False): continue
particles.append( RK4( m, dt, r0, v0 ))
particles[ len( particles ) - 1 ].radius = radius
break
#粒子外力の設定
for nm in range(NM):
particles[nm].fex = fex
####################################
## 描画用関数
####################################
#アニメーション作成用
ims=[]
#ピストン位置の配列
zs = []
def plot( t ):
#各コマを描画
for nm in range(NM):
markersize = particles[nm].radius * 1000
colorRG = 1.0 - particles[nm].m / m_max
if( nm == 0 ): img = plt.plot([particles[nm].r[0]], [particles[nm].r[2]], color = (colorRG, colorRG, 1.0), marker = 'o', markersize = markersize, linestyle='solid', linewidth = 0 )
else: img += plt.plot([particles[nm].r[0]], [particles[nm].r[2]], color = (colorRG, colorRG, 1.0), marker = 'o', markersize = markersize, linestyle='solid', linewidth = 0 )
time = plt.text(0.9, 1.12, "t = " + str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
#テキストをグラフに追加
img.append( time )
ims.append( img )
E = T() + U()
print( "t = {:.4f}".format(t)," E =", E )
####################################
## 各時刻における運動方程式の計算
####################################
for tn in range(len(ts)):
t = ts[tn]
if( tn%skip == 0 ):
plot(t)
#時間発展前の位置と速度を保持
for nm in range(NM):
particles[nm]._r = particles[nm].r.copy()
particles[nm]._v = particles[nm].v.copy()
level = 0
#時間発展の計算
result = timeEvolution( t, 1 )
#複数衝突を検知した場合
if( result == False ):
level += 1
highPrecisionMode(t, level )
#################################################################
## 分子運動アニメーション
#################################################################
plt.title( u"箱の中の粒子の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)
#描画範囲を設定
plt.xlim([-0.1, 1.1])
plt.ylim([-0.1, 1.1])
#箱の概形
plt.vlines([0-lw], -2*lw, L + 2*lw, "black", linestyles='solid', linewidth = 10 )
plt.vlines([L+lw], -2*lw, L + 2*lw, "black", linestyles='solid', linewidth = 10 )
plt.hlines([-lw], 0, L, "black", linestyles='solid', linewidth = 10 )
plt.hlines([L+lw], 0, L, "black", linestyles='solid', linewidth = 10 )
#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output100.mp4", writer="ffmpeg", dpi=300)
#グラフの表示
plt.show()
物理シミュレーションが趣味なシミュ君のページ