############################################################################
# 強制振り子運動シミュレーション
############################################################################
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)
######################################
# ４次のルンゲ・クッタ クラス
######################################
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):
#１段目
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
#################################################################
## 計算開始
#インスタンス
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()
物理シミュレーションが趣味なシミュ君のページ