【Python】重力＋張力による滑車振子運動
今回は下図のようなレール上を移動する滑車につながれたおもりの運動をシミュレーションするよ。まずはそのために必要な計算アルゴリズムの導出を行うよ。
滑車とおもりに加わる力は次のように定義するよ。
$$
\boldsymbol{F}_{\rm box} =m_{\rm box}\, \boldsymbol{g} - \boldsymbol{S} +\boldsymbol{N}_{\rm box}+ \boldsymbol{f}^{(\rm ex)}_{\rm box}
$$
$$
\boldsymbol{F} =m \boldsymbol{g} + \boldsymbol{S}+ \boldsymbol{f}^{\rm (ex)}
$$
$${\boldsymbol{S}}$$は張力で、$${\boldsymbol{f}^{(\rm ex)}_{\rm box}}$$と$${\boldsymbol{f}^{(\rm ex)} }$$は後に倒立振子運動をシミュレーションすることも考慮して加えておく滑車とおもりに外部から加える力（外力）だよ。$${\boldsymbol{N}_{\rm box}}$$はレールが滑車に与える垂直抗力だよ。今回は直線を想定してベクトルを$${\hat{\boldsymbol{t}}}$$として、
$$
\boldsymbol{N}_{\rm box} = -m_{\rm box}\, \boldsymbol{g} + \boldsymbol{S} +m_{\rm box}( \boldsymbol{g} \cdot \hat{\boldsymbol{t}})\, \hat{\boldsymbol{t}} - ( \boldsymbol{S} \cdot \hat{\boldsymbol{t}})\,\hat{\boldsymbol{t}}
$$
と考えられるね。まだ張力$${\boldsymbol{S}}$$は未知だけれども、滑車とおもりの加速度は次のように得られるね。
$$
\boldsymbol{a}_{\rm box} =\frac{\boldsymbol{F}_{\rm box}}{m_{\rm box}} =\left[ \left( \boldsymbol{g}-\frac{ \boldsymbol{S} }{ m_{\rm box}} \right)\cdot \hat{\boldsymbol{t}}\right]\hat{\boldsymbol{t}} + \frac{\boldsymbol{f}^{(\rm ex)}_{\rm box}}{m_{\rm box}}
$$
$$
\boldsymbol{a} =\frac{\boldsymbol{F}}{m} = \boldsymbol{g} + \frac{\boldsymbol{S}}{m} + \frac{\boldsymbol{f}^{\rm (ex)}}{m}
$$
この未知な$${\boldsymbol{S}}$$は単振子のときと同様、紐の長さが一定という拘束条件から導くことができるよ。
拘束条件（位置・速度・加速度）
$$
|\boldsymbol{r}-\boldsymbol{r}_{\rm box}| = L
$$
$$
(\boldsymbol{r}-\boldsymbol{r}_{\rm box}) \cdot (\boldsymbol{v}-\boldsymbol{v}_{\rm box}) = 0
$$
$$
(\boldsymbol{r}-\boldsymbol{r}_{\rm box}) \cdot (\boldsymbol{a}-\boldsymbol{a}_{\rm box}) + |\boldsymbol{v}-\boldsymbol{v}_{\rm box}|^2 =0
$$
この拘束条件に先の加速度を代入して、$${\boldsymbol{S}}$$を導くよ。
張力$\boldsymbol{S}$の導出
代入して少し整理した結果は次のとおりだよ。
$$
\boldsymbol{g} \cdot \boldsymbol{L}-(\boldsymbol{g} \cdot \hat{\boldsymbol{t}}) (\boldsymbol{L} \cdot \hat{\boldsymbol{t}}) + \frac{\boldsymbol{L} \cdot \boldsymbol{S}}{m}+ \frac{ (\boldsymbol{S} \cdot \hat{\boldsymbol{t}})(\boldsymbol{L} \cdot \hat{\boldsymbol{t}}) }{m_{\rm box}} + \frac{\boldsymbol{L} \cdot \boldsymbol{f}^{\rm (ex)}}{m} - \frac{\boldsymbol{L} \cdot\boldsymbol{f}^{(\rm ex)}_{\rm box}}{m_{\rm box}} +|\boldsymbol{v}-\boldsymbol{v}_{\rm box}|^2=0
$$
ただし、
$$
\boldsymbol{L} := \boldsymbol{r}-\boldsymbol{r}_{\rm box}
$$
と定義するね。なお張力$\boldsymbol{S}$の方向は必ずこの$${\boldsymbol{L} }$$と一致するはずだから、
$$
\boldsymbol{S} =S\, \frac{\boldsymbol{L} }{|\boldsymbol{L} |} = \frac{S}{L}\, \boldsymbol{L}
$$
と表すことができるね。これを代入して張力の大きさ$$S$$について解くと次のようになるね。
$$
S =- \frac{Lm m_{\rm box}}{ m ( \boldsymbol{L}\cdot \hat{\boldsymbol{t}})^2+m_{\rm box}L^2 } \left[ |\boldsymbol{v}-\boldsymbol{v}_{\rm box}|^2+ (\boldsymbol{g} \cdot \boldsymbol{L}) - (\boldsymbol{g} \cdot \hat{\boldsymbol{t}}) (\boldsymbol{L} \cdot \hat{\boldsymbol{t}}) + \frac{\boldsymbol{L} \cdot \boldsymbol{f}^{\rm (ex)}}{m} - \frac{\boldsymbol{L} \cdot\boldsymbol{f}^{(\rm ex)}_{\rm box}}{m_{\rm box}} \right]
$$
以上の結果から張力ベクトルは次のようになるね。
【計算アルゴリズム】滑車とおもり間に働く張力ベクトル
$$
\boldsymbol{S} =- \frac{m m_{\rm box}}{ m ( \boldsymbol{L}\cdot \hat{\boldsymbol{t}})^2+m_{\rm box}L^2 } \left[ |\boldsymbol{v}-\boldsymbol{v}_{\rm box}|^2+ (\boldsymbol{g} \cdot \boldsymbol{L}) - (\boldsymbol{g} \cdot \hat{\boldsymbol{t}}) (\boldsymbol{L} \cdot \hat{\boldsymbol{t}}) + \frac{\boldsymbol{L} \cdot \boldsymbol{f}^{\rm (ex)}}{m} - \frac{\boldsymbol{L} \cdot\boldsymbol{f}^{(\rm ex)}_{\rm box}}{m_{\rm box}} \right] \boldsymbol{L}
$$
上記アルゴリズムをもとに水平なレール上を運動する滑車につながれたおもりのシミュレーション結果を示すよ。
滑車振子運動シミュレーション #python #古典力学 pic.twitter.com/5vKdq3XC0o
— シミュ君 (@simu_kun) January 23, 2025
Pythonプログラムソース
Pythonプログラムソースを公開するよ。もし良かったら試してみてくださーい！
#################################################################
TITLE = "張力＋重力による滑車振子運動"
#################################################################
import os
import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import colorsys
###############################
# 実験室グローバル変数
###############################
TIME_START = 0.0 #計算開始時刻
TIME_END = 5.0 #計算終了時刻
DT = 0.01 #時間間隔
PLOT_FLAG = False #グラフ描画
ANIMATION_FLAG = True #アニメーション作成
ANIMATION_INTERVAL = 5 #アニメーション作成時の描画インターバル
OUTPUT_FLAG = False
OUTPUT_PASS = "./data_20250122/"
if( OUTPUT_FLAG and os.path.isdir(OUTPUT_PASS) == False) : os.mkdir(OUTPUT_PASS)
###############################
# 計算グローバル変数
###############################
#時刻
times = np.arange(TIME_START, TIME_END, DT)
#重力加速度
vec_g = np.array([0, 0, -10])
#質量
m = 1.0
#m_box = 5.0
#ひもの長さ
L0 = 3
#初期条件
vec_r0 = np.array([0, 0, -L0])
#vec_v0 = np.array([12, 0, 0])
vec_r0_box = np.array([0, 0, 0])
vec_v0_box = np.array([0, 0, 0])
#外力
vec_F = np.array([0, 0, 0])
vec_F_box = np.array([0, 0, 0])
#レールの法線ベクトル
hat_t = np.array([1, 0, 0.1])
hat_t = hat_t / np.sqrt( np.dot( hat_t, hat_t ))
def rail( x ):
y = hat_t[2]/hat_t[0] * x
return y
#滑車質量の違い
m_boxs = [3] #np.arange(1, 10, 0.05)
hat_n = np.array([0, 0, 1])
hat_n = hat_n / np.sqrt( np.dot( hat_n, hat_n ))
###############################
# 関数定義
###############################
#張力＋重力の微分方程式
def equations(t, X, m_box):
vec_r_box = np.array([ X[0], X[1], X[2] ])
vec_v_box = np.array([ X[3], X[4], X[5] ])
vec_r = np.array([ X[6], X[7], X[8] ])
vec_v = np.array([ X[9], X[10], X[11] ])
#相対位置・速度ベクトル
bar_r = vec_r - vec_r_box
bar_v = vec_v - vec_v_box
#紐の長さ
L = np.sqrt( np.dot(bar_r, bar_r) )
#相対速度の大きさ
v = np.sqrt( np.dot(bar_v, bar_v) )
#単位ベクトル
#hat_z = np.array([ 0, 0, 1 ])
#張力
vec_S = -m * m_box / ( np.dot(bar_r, hat_t)**2 * m + L**2 * m_box ) * ( v**2 + np.dot(vec_g, bar_r) - np.dot(vec_g, hat_t) * np.dot(bar_r, hat_t) + np.dot(vec_F, bar_r )/m - np.dot(vec_F_box, bar_r )/m_box ) * bar_r
#質点に加わる力
F_box = np.dot( m_box * vec_g - vec_S, hat_t) * hat_t + vec_F_box
F = m * vec_g + vec_S + vec_F
#ニュートンの運動方程式
vec_a_box = F_box / m_box
vec_a = F / m
return [
vec_v_box[0], vec_v_box[1], vec_v_box[2], vec_a_box[0], vec_a_box[1], vec_a_box[2],
vec_v[0], vec_v[1], vec_v[2], vec_a[0], vec_a[1], vec_a[2]
]
# theta 依存性を可視化
def calculation(vec_r0_box, vec_v0_box, vec_r0, vec_v0, m_box ):
solution = integrate.solve_ivp(
equations, #常微分方程式
method = 'RK45', #計算方法
args = (m_box,), #引数
t_span=( times[0], times[-1]), #積分区間
rtol = 1.e-9, #相対誤差
atol = 1.e-9, #絶対誤差
y0 = np.concatenate([vec_r0_box, vec_v0_box, vec_r0, vec_v0]), #初期値
t_eval = times #結果を取得する時刻
)
#計算結果
vec_r_box_t = np.array( [solution.y[0], solution.y[1], solution.y[2]] )
vec_v_box_t = np.array( [solution.y[3], solution.y[4], solution.y[5]] )
vec_r_t = np.array( [solution.y[6], solution.y[7], solution.y[8]] )
vec_v_t = np.array( [solution.y[9], solution.y[10], solution.y[11]] )
return vec_r_box_t, vec_v_box_t, vec_r_t, vec_v_t
##########################################################
if __name__ == "__main__":
#描画用リスト
xls_box = []
yls_box = []
zls_box = []
vxls_box = []
vyls_box = []
vzls_box = []
xls = []
yls = []
zls = []
vxls = []
vyls = []
vzls = []
lls = []
for i in range( len(m_boxs) ):
m_box = m_boxs[ i ]
vec_v0 = np.array([10, 5, 0])
vec_v0_box = np.array([0, 0, 0])
vec_r_box_t, vec_v_box_t, vec_r_t, vec_v_t = calculation( vec_r0_box, vec_v0_box, vec_r0, vec_v0, m_box )
xls_box.append( vec_r_box_t[0] )
yls_box.append( vec_r_box_t[1] )
zls_box.append( vec_r_box_t[2] )
vxls_box.append( vec_v_box_t[0] )
vyls_box.append( vec_v_box_t[1] )
vzls_box.append( vec_v_box_t[2] )
xls.append( vec_r_t[0] )
yls.append( vec_r_t[1] )
zls.append( vec_r_t[2] )
vxls.append( vec_v_t[0] )
vyls.append( vec_v_t[1] )
vzls.append( vec_v_t[2] )
L = np.sqrt( (vec_r_t[0][-1] - vec_r_box_t[0][-1])**2 + (vec_r_t[1][-1] - vec_r_box_t[1][-1])**2 + (vec_r_t[2][-1] - vec_r_box_t[2][-1])**2 )
print(i, np.abs(L-L0))
if( OUTPUT_FLAG ):
for i in range(len(times)):
path = OUTPUT_PASS + str( i ) + '.txt'
f = open(path, 'w')
for j in range( len(xls) ):
f.write(
str(round( xls_box[ j ][ i ], 2)) + "\t" + str(round( yls_box[ j ][ i ], 2)) + "\t" + str(round( zls_box[ j ][ i ], 2)) + "\t" + str(round( vxls_box[ j ][ i ], 2)) + "\t" + str(round( vyls_box[ j ][ i ], 2)) + "\t" + str(round( vzls_box[ j ][ i ], 2)) + "\t" +
str(round( xls[ j ][ i ], 2)) + "\t" + str(round( yls[ j ][ i ], 2)) + "\t" + str(round( zls[ j ][ i ], 2)) + "\t" + str(round( vxls[ j ][ i ], 2)) + "\t" + str(round( vyls[ j ][ i ], 2)) + "\t" + str(round( vzls[ j ][ i ], 2)) + "\n"
)
f.close()
if( PLOT_FLAG ):
#################################################################
# グラフの設定
#################################################################
fig = plt.figure( figsize=(12, 8) ) #画像サイズ
plt.subplots_adjust(left = 0.085, right = 0.97, bottom = 0.10, top = 0.95) #グラフサイズ
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']
#グラフタイトルと軸ラベルの設定
plt.title( TITLE, fontname="Yu Gothic", fontsize=20, fontweight=500 )
plt.xlabel( r"$t \,[{\rm s}]$", fontname="Yu Gothic", fontsize=20, fontweight=500 )
plt.ylabel( r"$x \,[{\rm m}]$" , fontname="Yu Gothic", fontsize=20, fontweight=500 )
#罫線の描画
plt.grid( which = "major", axis = "x", alpha = 0.7, linewidth = 1 )
plt.grid( which = "major", axis = "y", alpha = 0.7, linewidth = 1 )
#描画範囲を設定
plt.xlim([0, TIME_END])
#plt.ylim([-2.0, 1.0])
plt.plot(times, xls_box[0], linestyle='solid', linewidth = 2, color = colors[0])
plt.plot(times, zls_box[0], linestyle='solid', linewidth = 2, color = colors[0])
plt.plot(times, xls[0], linestyle='solid', linewidth = 2, color = colors[1])
plt.plot(times, zls[0], linestyle='solid', linewidth = 2, color = colors[1])
#グラフの表示
plt.show()
if( ANIMATION_FLAG ):
#################################################################
# アニメーションの設定
#################################################################
fig = plt.figure( figsize=(12, 12) ) #画像サイズ
plt.subplots_adjust(left = 0.085, right = 0.97, bottom = 0.10, top = 0.95) #グラフサイズ
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']
#アニメーション作成用
imgs = []
#アニメーション分割数
AN = len( times )
#アニメーションの各グラフを生成
for an in range( AN ):
t = DT * an
if( an % ANIMATION_INTERVAL != 0): continue
img = plt.plot([-5, 5], [rail( -5 ), rail( 5 )], color = 'gray', linestyle='dotted', linewidth = 2 )
for i in range( len( m_boxs ) ):
c = colorsys.hls_to_rgb( i/len(m_boxs) , 0.5 , 1)
img += plt.plot([xls_box[i][an], xls[i][an] ], [zls_box[i][an], zls[i][an]], color = c, linestyle='solid', linewidth = 2 )
img += plt.plot([xls_box[i][an]], [zls_box[i][an]], color = c, marker = 's', markersize = 5)
img += plt.plot([xls[i][an]], [zls[i][an]], color = c, marker = 'o', markersize = 20)
#時刻をグラフに追加
time = plt.text(13.5, 0.35, "t = " + str("{:.2f}".format(t)), fontsize=18)
img.append( time )
imgs.append( img )
#グラフタイトルと軸ラベルの設定
plt.title( TITLE, fontname="Yu Gothic", fontsize=20, fontweight=500 )
plt.xlabel( r"$x \,[{\rm m}]$", fontname="Yu Gothic", fontsize=20, fontweight=500 )
plt.ylabel( r"$y \,[{\rm m}]$" , fontname="Yu Gothic", fontsize=20, fontweight=500 )
#罫線の描画
plt.grid( which = "major", axis = "x", alpha = 0.7, linewidth = 1 )
plt.grid( which = "major", axis = "y", alpha = 0.7, linewidth = 1 )
#描画範囲を設定
plt.xlim([-5, 5])
plt.ylim([-5, 5])
#アニメーションの生成
ani = animation.ArtistAnimation(fig, imgs, interval=50)
#アニメーションの保存
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=360)
#グラフの表示
plt.show()