任意結合張力による運動：９連結子による角運動量保存運動
前々回に導出した任意結合張力による運動の例として、重力なしで９個の質点を繋いで角運動量が一定になるように初速度を与えたときの運動をシミュレーションするよ。
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
INF = float('inf')
###############################
# 実験室グローバル変数
###############################
TIME_START = 0.0 #計算開始時刻
TIME_END = 1.0 #計算終了時刻
DT = 0.01 #時間間隔
PLOT_FLAG = True #グラフ描画
ANIMATION_FLAG = True #アニメーション作成
ANIMATION_INTERVAL = 1 #アニメーション作成時の描画インターバル
OUTPUT_FLAG = False
OUTPUT_PASS = "./data_20241004/"
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, 0])
#質量
ml = [1.0, 3.0, 1.0, 5.0, 1.0, 1.0, 4.0, 1.0, 3.0]
NumberOfFix = 0
#結合リスト
connections = []
connections.append( [ 0, 1 ] )
connections.append( [ 1, 2 ] )
connections.append( [ 2, 3 ] )
connections.append( [ 3, 4 ] )
connections.append( [ 4, 5 ] )
connections.append( [ 5, 6 ] )
connections.append( [ 6, 7 ] )
connections.append( [ 7, 8 ] )
#初期位置と速度
vec_r = [
np.array([ 0, 0, -4.0 ]),
np.array([ 0, 0, -3.0 ]),
np.array([ 0, 0, -2.0 ]),
np.array([ 0, 0, -1.0 ]),
np.array([ 0, 0, 0.0 ]),
np.array([ 0, 0, 1.0 ]),
np.array([ 0, 0, 2.0 ]),
np.array([ 0, 0, 3.0 ]),
np.array([ 0, 0, 4.0 ])
]
vec_v = [
np.array([ 0, 0, 0]),
np.array([ 0, 0, 0]),
np.array([ 0, 0, 0]),
np.array([ 0, 0, 0]),
np.array([ 0, 0, 0]),
np.array([ 0, 0, 0]),
np.array([ 0, 0, 0]),
np.array([ 0, 0, 0]),
np.array([ 0, 0, 0])
]
#ひもの長さ（初期値）
L0 = []
for n in range(len(connections)):
i = connections[ n ][ 0 ]
j = connections[ n ][ 1 ]
ri = vec_r[ i - NumberOfFix ]
rj = vec_r[ j - NumberOfFix ]
rij = ri - rj
L0.append( np.sqrt( np.dot( rij, rij ) ) )
def box():
return
###############################
# 関数定義
###############################
#相対誤差
rtol = 1e-5
#絶対誤差
atol = 1e-5
#補正ばね弾性係数
compensationK = 0.06
#補正粘性抵抗係数
compensationGamma = 0.5
#接続番号を取得
def getConectionNumber( i, j ):
if( i == j ): return -1
for m in range( len(connections) ):
connection = connections[ m ]
if( i in connection and j in connection ): return m
return -1
#張力＋重力の微分方程式
def equations(t, X, compensationK, compensationGamma):
vec_r = [0] * len(ml)
vec_v = [0] * len(ml)
vec_a = [0] * len(ml)
j = 0
for i in range(len(ml)):
if(ml[ i ] == INF):
#支点の位置ベクトル・速度ベクトル・加速度ベクトル
vec_r[ i ], vec_v[ i ], vec_a[ i ] = box( t )
else:
vec_a[ i ] = vec_g
vec_r[ i ] = np.array([ X[ 6 * j + 0 ], X[ 6 * j + 1 ], X[ 6 * j + 2 ] ])
vec_v[ i ] = np.array([ X[ 6 * j + 3 ], X[ 6 * j + 4 ], X[ 6 * j + 5 ] ])
j += 1
N = len(connections)
vec_rij = [ 0 ] * N
vec_vij = [ 0 ] * N
vec_aij = [ 0 ] * N
vec_nij = [ 0 ] * N
Lij = [ 0 ] * N
vij = [ 0 ] * N
for n in range(N):
i = connections[ n ][ 0 ]
j = connections[ n ][ 1 ]
#相対位置ベクトル
vec_rij[ n ] = vec_r[ i ] - vec_r[ j ]
#ひもの長さ
Lij[ n ] = np.sqrt( np.dot(vec_rij[n], vec_rij[n]) )
#相対単位ベクトル
vec_nij[ n ] = vec_rij[ n ] / Lij[ n ]
#相対速度ベクトル
vec_vij[ n ] = vec_v[ i ] - vec_v[ j ]
#相対速度の大きさ
vij[ n ] = np.sqrt( np.dot(vec_vij[n], vec_vij[n]))
#相対加速度ベクトル
vec_aij[ n ] = vec_a[ i ] - vec_a[ j ]
#a係数（右辺）
A = np.array( [ 0.0 ] * N )
for n in range( N ):
A[ n ] = -1.0 / Lij[ n ] * ( vij[ n ]**2 + np.dot(vec_aij[ n ], vec_rij[ n ]) )
#f成分（左辺）
F = np.zeros(( N, N ))
for n in range( N ):
i = connections[ n ][ 0 ]
j = connections[ n ][ 1 ]
for k in range( len( ml ) ):
nki = getConectionNumber( k, i )
nkj = getConectionNumber( k, j )
if( nki > -1):
vec_n = vec_nij[ nki ]
if( i < k ) : vec_n = -vec_n
F[ n ][ nki ] += np.dot( vec_n, vec_nij[ n ]) / ml[ i ]
if( nkj > -1):
vec_n = vec_nij[ nkj ]
if( j < k ) : vec_n = -vec_n
F[ n ][ nkj ] -= np.dot( vec_n, vec_nij[ n ]) / ml[ j ]
#print( np.linalg.det(F) )
# 求解
vec_f = np.linalg.solve(F, A)
#補正力
vec_fc = [0.0] * N
for n in range(N):
i = connections[ n ][ 0 ]
j = connections[ n ][ 1 ]
vec_fc[ n ] = (compensationK * ( Lij[ n ] - L0[ n ] ) + compensationGamma * np.dot(vec_vij[ n ], vec_nij[ n ]) ) * vec_nij[ n ]
if( i > j ): vec_fc[ n ] = -vec_fc[ n ]
#補正因子
vec_fc[ n ] *= np.abs(vec_f[ n ])
#質点に加わる力
vec_F = [0.0] * len(ml)
for i in range(len(ml)):
if(ml[ i ] == INF): continue
vec_F[ i ] = ml[ i ] * vec_g
for j in range( len( ml )):
nn = getConectionNumber( i, j )
if( nn > -1):
f = vec_f[ nn ] * vec_nij[ nn ] + vec_fc[ nn ]
if( i < j ): f = -f
vec_F[ i ] += f
#ニュートンの運動方程式
vec_a = [0.0] * len(ml)
for i in range(len(ml)):
if(ml[ i ] == INF): continue
vec_a[ i ] = vec_F[ i ] / ml[ i ]
H = []
for i in range(len(ml)):
if(ml[ i ] == INF): continue
H.append( vec_v[ i ][ 0 ])
H.append( vec_v[ i ][ 1 ])
H.append( vec_v[ i ][ 2 ])
H.append( vec_a[ i ][ 0 ])
H.append( vec_a[ i ][ 1 ])
H.append( vec_a[ i ][ 2 ])
return H
def calculation(vec_r, vec_v, compensationK, compensationGamma):
init = []
for i in range(len(vec_r)):
init.append( vec_r[i] )
init.append( vec_v[i] )
solution = integrate.solve_ivp(
equations, #常微分方程式
method = 'RK45', #計算方法
args = (compensationK, compensationGamma), #引数
t_span=( times[0], times[-1]), #積分区間
rtol = rtol, #相対誤差
atol = atol, #絶対誤差
y0 = np.concatenate(init), #初期値
t_eval = times #結果を取得する時刻
)
#計算結果
vec_r_t = []
vec_v_t = []
for i in range(len(vec_r)):
vec_r_t.append( np.array( [solution.y[ 6 * i + 0 ], solution.y[ 6 * i + 1 ], solution.y[ 6 * i + 2 ]] ) )
vec_v_t.append( np.array( [solution.y[ 6 * i + 3 ], solution.y[ 6 * i + 4 ], solution.y[ 6 * i + 5 ]] ) )
return vec_r_t, vec_v_t
##########################################################
if __name__ == "__main__":
vec_v[0][0] = 40 * 3
vec_v[8][0] = -40
vec_r_t, vec_v_t = calculation( vec_r, vec_v, compensationK, compensationGamma)
DeltaLs = []
ts = []
for tn in range(len(times)):
t = tn * DT
Lij = [0.0] * len( connections )
for cn in range(len( connections )) :
i = connections[ cn ][ 0 ]
j = connections[ cn ][ 1 ]
ri = np.array( [ vec_r_t[ i - NumberOfFix ][ 0 ][ tn ], vec_r_t[ i - NumberOfFix ][ 1 ][ tn ], vec_r_t[ i - NumberOfFix ][ 2 ][ tn ] ])
rj = np.array( [ vec_r_t[ j - NumberOfFix ][ 0 ][ tn ], vec_r_t[ j - NumberOfFix ][ 1 ][ tn ], vec_r_t[ j - NumberOfFix ][ 2 ][ tn ] ])
rij = ri - rj
Lij[ cn ] = np.sqrt( np.dot( rij, rij ) )
DeltaL = 0
for cn in range(len( connections )) :
DeltaL += np.abs( Lij[ cn ] - L0[ cn ] ) / len( connections )
if(tn%10==9):
print( DeltaL )
ts.append( t )
DeltaLs.append( DeltaL )
if( OUTPUT_FLAG ):
for tn in range(len(times)):
t = tn * DT
vec_r_box, vec_v_box, vec_a_box = box( t )
path = OUTPUT_PASS + str( tn ) + '.txt'
f = open(path, 'w')
for j in range( len(ml) ):
if( ml[j] == INF ):
for xyz in range(3):
f.write( str( round( vec_r_box[ xyz ], 2)) + "\t")
for xyz in range(3):
f.write( str( round( vec_v_box[ xyz ], 2)) + "\t")
f.write("\n")
else:
for xyz in range(3):
f.write( str( round( vec_r_t[ j - NumberOfFix][ xyz ][ tn ], 2)) + "\t")
for xyz in range(3):
f.write( str( round( vec_v_t[ j - NumberOfFix][ xyz ][ tn ], 2)) + "\t")
f.write("\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"$|\Delta L| \,[{\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([0, 2.0])
plt.yscale("log")
#グラフの描画
plt.plot(ts, DeltaLs, linestyle='solid', linewidth = 4)
#グラフの表示
plt.show()
if( ANIMATION_FLAG ):
#################################################################
# アニメーションの設定
#################################################################
fig = plt.figure( figsize=(8, 6) ) #画像サイズ
plt.subplots_adjust(left = 0.095, right = 0.97, bottom = 0.10, top = 0.92) #グラフサイズ
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([ 0 ], [ 0 ], marker = 'o', markersize = 0.01, color = colors[0])
c = colorsys.hls_to_rgb( 0.2 , 0.5 , 1)
NumberOfFix = 0
for cn in range(len( connections )):
i = connections[ cn ][ 0 ]
j = connections[ cn ][ 1 ]
img += plt.plot([ vec_r_t[ i - NumberOfFix ][ 0 ][ an ], vec_r_t[ j - NumberOfFix ][ 0 ][ an ] ], [vec_r_t[ i - NumberOfFix ][ 2 ][ an ], vec_r_t[ j - NumberOfFix ][ 2 ][ an ]], linestyle='solid', linewidth = 1, color = colors[0])
NumberOfFix = 0
for i in range(len(ml)):
m = ml[i]
img += plt.plot([ vec_r_t[ i - NumberOfFix ][ 0 ][ an ] ], [vec_r_t[ i - NumberOfFix ][ 2 ][ an ]], marker = 'o', markersize = 5 * m, color = colors[i])
#時刻をグラフに追加
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"$z \,[{\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([-4.2, 4.2])
plt.ylim([-4.2, 4.2])
#アニメーションの生成
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()
この記事が気に入ったらサポートをしてみませんか？