【第4話】強制振り子の計算アルゴリズム【ルンゲクッタで行こう!④】

############################################################################
# 強制振り子運動シミュレーション
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
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 = 1.0

#支点の初期位置
box_r0 = np.array([0.0, 0.0, 0.0])

#初期位置と速度
r0 = np.array([0.0, 0.0, -1.0])
v0 = np.array([0.0, 0.0, 0.0])

#位置ベクトル
L0 = r0 - box_r0
#ひもの長さ
L0_abs = np.sqrt( np.sum( L0**2 ) )

#強制振動の大きさと角振動数
box_l = 0.2
box_omega = 0.5 * (2.0 * np.pi) 

#支点の運動関数
def boxFunction( t ):
	box_rx = box_r0[0] + box_l * np.sin( box_omega * t)
	box_vx = box_l * box_omega * np.cos( box_omega * t)
	box_ax = - box_l * box_omega**2 * np.sin( box_omega * t)
	box_r = np.array([box_rx, 0.0, 0.0])
	box_v = np.array([box_vx, 0.0, 0.0])
	box_a = np.array([box_ax, 0.0, 0.0])
	return box_r, box_v, box_a


#重力加速度
g = np.array([0.0, 0.0, -10.0])
#補正パラメータ
compensationK = 0.0
compensationGamma = 0.0

#張力ベクトル
def S ( t, r_t, v_t):
	#支点の位置ベクトル
	box_r, box_v, box_a = boxFunction(t)	
	#おもり位置ベクトル
	L_t = r_t - box_r
	#おもり速度ベクトル
	v_bar = v_t - box_v
	#ひもの長さ
	L_abs = np.sqrt( np.sum(L_t**2) )

	return -m / L_abs ** 2 * ( np.sum(v_bar**2) - np.dot(box_a, L_t) + np.dot(g, L_t) ) * L_t

def F ( t, r_t, v_t ):
	#支点の位置ベクトル
	box_r, box_v, box_a = boxFunction(t)

	#おもり位置ベクトル
	L_t = r_t - box_r
	#おもり速度ベクトル
	v_bar = v_t - box_v
	#ひもの長さ
	L_abs = np.sqrt( np.sum( L_t**2 ) )
	#おもり速度の大きさ
	v_bar_abs = np.sqrt( np.sum( v_bar**2 ) )
	#単位方向ベクトル
	L_hat = L_t.copy() / L_abs
	#張力ベクトル
	S_ = S( t, r_t, v_t ) 
	#張力の大きさ
	S_coefficient = np.sqrt( np.sum( S_**2 ) )

	#補正力ベクトル
	compensationK = S_coefficient 
	compensationGamma = v_bar_abs * 10
	fc1 = - compensationK * ( L_abs - L0_abs ) * L_hat
	fc2 = - compensationGamma * np.dot(v_bar, L_hat) * L_hat

	#合力ベクトル
	return  S_ + m * g + fc1 + fc2


#時間区間
t_min = 0
t_max = 40.0
#時間区間数
NT = 4000
skip =  2
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)

######################################
# 4次のルンゲ・クッタ クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		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.__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( t, r_t, v_t ) / m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#1段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#2段目
		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			
		)
		#3段目
		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			
		)
		#4段目
		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 

#################################################################
## 計算開始

#インスタンス
rk4 = RK4(dt, r0, v0)

#アニメーション作成用
ims=[]
def plot( t, r ):
	box_r, box_v, box_a = boxFunction(t)
	print( "{:.2f}".format(t),  np.sqrt((r[0]-box_r[0])**2 + (r[2]-box_r[2])**2) - 1  )

	#各コマを描画
	img  = plt.plot([box_r[0]], [box_r[2]], colors[0], marker = 's', markersize = 10, linestyle='solid', linewidth = 0 )
	img  += plt.plot([r[0]], [r[2]], colors[1], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
	img  += plt.plot([box_r[0], r[0]], [box_r[2], r[2]], colors[1], linestyle='solid', linewidth = 1 )
	time = plt.text(0.8, 1.25, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )

#各時刻における運動方程式の計算
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t, rk4.r )
	#ルンゲ・クッタ法による時間発展
	rk4.timeEvolution( t )
	#位置と速度を更新
	rk4.r += rk4.dr
	rk4.v += rk4.dv

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([-1.2, 1.2])
plt.ylim([-1.2, 1.2])

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.mp4", writer="ffmpeg", dpi=300)

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

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です