【Python】移動する円経路に束縛された質点運動アニメーション

これまでに導出した経路に束縛された質点運動の計算アルゴリズムは、経路が移動する場合でも成り立つね。導出したアルゴリズムの$${\boldsymbol{r}_0, \boldsymbol{v}_0, \boldsymbol{a}_0}$$を具体的に与えるだけだね。今回は円経路の中心座標が次の式で与えた円運動の場合を質点の運動をシミュレーションした結果を示すよ。

$$
\boldsymbol{r}_{\rm circle}(t) =\left\{ \begin{matrix} r\cos( \omega t ) \cr r\sin( \omega t ) \end{matrix} \right.
$$

$$
\boldsymbol{v}_{\rm circle}(t) =\left\{ \begin{matrix} - r \omega \sin( \omega t ) \cr r\omega\cos( \omega t ) \end{matrix} \right.
$$

$$
\boldsymbol{a}_{\rm circle}(t) =\left\{ \begin{matrix} - r \omega^2 \cos( \omega t ) \cr -r\omega^2\sin( \omega t ) \end{matrix} \right.
$$

シミュレーション結果

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   = 10.0       #計算終了時刻
DT = 0.01               #時間間隔
PLOT_FLAG = False        #グラフ描画
ANIMATION_FLAG = True  #アニメーション作成
ANIMATION_INTERVAL = 10 #アニメーション作成時の描画インターバル
OUTPUT_FLAG = False
OUTPUT_PASS = "./data_20241224/"
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

#初速度の大きさ
v0s = [0.0] #np.arange(0.0, 0.0002+0.000001, 0.000001)
###############################
# 関数定義
###############################
#相対誤差
rtol = 1e-9
#絶対誤差                   
atol = 1e-9                   
#補正ばね弾性係数
compensationK = 0.004
#補正粘性抵抗係数
compensationGamma = 0.4

# 経路の設定
class Path:
	#コンストラクタ
	def __init__(self):
		#円の半径
		self.R = 1
		#円経路の回転半径
		self.RR = 0.5
		#円経路の角速度
		self.omega = 2.0 * np.pi /2

	#円の中心位置ベクトル
	def center( self, time ):
		x = self.RR * np.cos(self.omega * time - np.pi/2)
		y = 0
		z = self.RR * np.sin(self.omega * time - np.pi/2) + self.R + self.RR
		return np.array([x, y, z])
	#円の中心速度ベクトル
	def velocity( self, time ):
		x = -self.RR * self.omega * np.sin(self.omega * time - np.pi/2)
		y = 0
		z = self.RR * self.omega * np.cos(self.omega * time - np.pi/2)
		return np.array([x, y, z])
	#円の中心加速度ベクトル
	def acceleration( self, time ):
		x = -self.RR * self.omega**2 * np.cos(self.omega * time - np.pi/2)
		y = 0
		z = -self.RR * self.omega**2 * np.sin(self.omega * time - np.pi/2)
		return np.array([x, y, z])

	#経路の位置ベクトルを指定する媒介変数関数
	def position( self, time, theta ):
		_center =  self.center( time ); 
		x = self.R * np.cos(theta) + _center[0]
		y = 0                      + _center[1]
		z = self.R * np.sin(theta) + _center[2]
		return np.array([x, y, z])

	#接線ベクトルを指定する媒介変数関数
	def tangent( self, time, theta ):
		x = -np.sin(theta)
		y = 0
		z = np.cos(theta)
		return np.array([x, y, z])
	#曲率ベクトルを指定する媒介変数関数
	def curvature( self, time, theta ):
		x = -np.cos(theta) / self.R
		y = 0
		z = -np.sin(theta) / self.R
		return np.array([x, y, z])
	#媒介変数の取得
	def getTheta( self, time, r ):
		_center =  self.center( time ); 
		#相対位置ベクトル
		bar_r = r - _center
		bar_r_length = np.sqrt(np.sum( bar_r**2 ))
		sinTheta = bar_r[2] / bar_r_length
		if( sinTheta > 0 ):
			theta = np.arccos( bar_r[0] / bar_r_length )
		else:
			theta = 2.0 * np.pi - np.arccos( bar_r[0] / bar_r_length )

		return theta


#経路拘束力の微分方程式
def equations(time, X, path, compensationK, compensationGamma):
	vec_r =  np.array([ X[0], X[1], X[2] ])
	vec_v =  np.array([ X[3], X[4], X[5] ])

	#媒介変数の取得
	theta = path.getTheta ( time, vec_r )

	#媒介変数に対する位置ベクトル、接線ベクトル、曲率ベクトルの計算
	position = path.position( time, theta )
	tangent = path.tangent( time, theta )
	curvature = path.curvature( time, theta )

	#経路中心の位置ベクトル・速度ベクトル・加速度ベクトル
	vec_r_box = path.center( time )
	vec_v_box = path.velocity( time )
	vec_a_box = path.acceleration( time )

	#相対速度
	bar_v = vec_v - vec_v_box

	#微係数dl/dtとd^2l/dt^2を計算
	dl_dt = np.dot( bar_v, tangent )
	d2l_dt2 = np.dot( vec_g, tangent ) - np.dot( vec_a_box, tangent) 

	#加速度を計算
	vec_a = vec_a_box + d2l_dt2 * tangent + dl_dt * dl_dt * curvature

	#ズレ位置ベクトル
	bar_r = vec_r - position
	#補正ばね弾性力
	fk = -compensationK * bar_r
	#法線ベクトル
	n1 = curvature / np.sqrt(np.sum( curvature**2 ))
	n2 = np.cross(n1, tangent )
	#補正粘性抵抗力
	fgamma1 = -compensationGamma * np.dot( bar_v, n1 ) * n1 
	fgamma2 = -compensationGamma * np.dot( bar_v, n2 ) * n2
	#経路補正力を加える
	vec_a += (fk + fgamma1 + fgamma2)/m * np.abs(m*vec_a)

	return [vec_v[0], vec_v[1], vec_v[2], vec_a[0], vec_a[1], vec_a[2]]

# 初速度依存性を可視化
def calculation(vec_r0, vec_v0, path, compensationK, compensationGamma):

	solution = integrate.solve_ivp(
		equations,                     #常微分方程式
		method = 'RK45',                #計算方法
		args = (path, compensationK, compensationGamma),              #引数
		t_span=( times[0], times[-1]),  #積分区間
		rtol = rtol,                   #相対誤差
		atol = atol,                   #絶対誤差
		y0 = np.concatenate([vec_r0, vec_v0]), #初期値
		t_eval = times                  #結果を取得する時刻
	)
	#計算結果
	return solution.y[ 0 ], solution.y[ 1 ], solution.y[ 2 ], solution.y[ 3 ], solution.y[ 4 ], solution.y[ 5 ]

##########################################################
if __name__ == "__main__":

	#経路の生成
	path = Path()
	#描画用リスト
	xls = []
	yls = []
	zls = []
	vxls = []
	vyls = []
	vzls = []
	#vo による違い 
	for i in range(len(v0s)):
		print( i )

		#初期位置と速度
		vec_r0 = np.array([ 0, 0, 0 ])
		vec_v0 = np.array([ v0s[i], 0, 0])

		x, y, z, vx, vy, vz = calculation( vec_r0, vec_v0, path, compensationK, compensationGamma )
		xls.append( x )
		yls.append( y )
		zls.append( z )
		vxls.append( vx )
		vyls.append( vy )
		vzls.append( vz )

		vec_r = np.array([ x[-1], y[-1], z[-1] ])
		vec_v = np.array([ vx[-1],vy[-1],vz[-1] ])

		_center = path.center( times[-1] )
		vecL = vec_r - _center
		L = np.sqrt( vecL[0]**2 + vecL[1]**2 + vecL[2]**2 )
		deltaL_abs = np.abs( L - path.R ) 

		print( "|L-L0| = ", deltaL_abs )
		

	if( OUTPUT_FLAG ):
		for i in range(len(times)):
			f_path = OUTPUT_PASS + str( i ) + '.txt'
			f = open(f_path, 'w')

			for j in range( len(xls) ):
				f.write( str(round( xls[ j ][ i ], 4)) + "\t" +  str(round( yls[ j ][ i ], 4)) + "\t" + str(round( zls[ j ][ i ], 4)) + "\t" +  str(round( vxls[ j ][ i ], 4)) + "\t" + str(round( vyls[ j ][ i ], 4)) + "\t" +  str(round( vzls[ j ][ i ], 4)) + "\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"$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([0, TIME_END])
		plt.ylim([0, 2])

		#グラフの描画
		for i in range(len(v0s)):
			c = colorsys.hls_to_rgb( i/len(v0s) , 0.5 , 1)
			plt.plot(times, zls[i], linestyle='solid', linewidth = 4)#, color = c)

		#グラフの表示
		plt.show()

	if( ANIMATION_FLAG ):

		#################################################################
		# アニメーションの設定
		#################################################################
		fig = plt.figure(figsize=(8, 8)) #画像サイズ
		plt.subplots_adjust(left = 0.07, right = 0.98, bottom = 0.07, 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 )

		#軌跡リスト
		r_trajectorys = [ 0.0 ] * len( v0s )

		#アニメーションの各グラフを生成
		for an in range( AN ):
			t = DT * an

			if( an % ANIMATION_INTERVAL != 0): continue


			#円座標点配列の生成
			thetas = np.linspace(0, 2*np.pi, 200)

			_center = path.center( t )
			cirsx = np.cos( thetas ) + _center[0]
			cirsy = np.sin( thetas ) + _center[2]


			#軌跡の描画
			img = plt.plot(cirsx, cirsy, color = "black", linestyle='dotted', linewidth = 1 )

			for i in range( len( v0s ) ):
				c = colorsys.hls_to_rgb( i/len(v0s) , 0.5 , 1)
				if( i == 0 ): 
					img += plt.plot([xls[i][an]], [zls[i][an]], color = c, marker = 'o', markersize = 20)
				else:
					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"$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([-1.6, 1.6])
		plt.ylim([-0.2, 3.2])
		#アニメーションの生成
		ani = animation.ArtistAnimation(fig, imgs, interval=50)
		#アニメーションの保存
		#ani.save("output.html", writer=animation.HTMLWriter())
		#ani.save("output.gif", writer="imagemagick")
		#ani.save("output_30.mp4", writer="ffmpeg", dpi=360)
		#グラフの表示
		plt.show()