【Python】ばね弾性力＋重力＋粘性抵抗力＋強制振動による運動（強制振動運動）
今日も古典力学の復習を進めていくよ。今回は前回の「ばね弾性力＋重力＋粘性抵抗力」に加わった物体を強制的に振動させるの強制振動運動だよ。$${\boldsymbol{r}_{\rm box}(t)}$$は支点の位置で、支点を振動させることでばね弾性力を通じて強制振動させることを考えるよ。重力＋ばね弾性力＋粘性抵抗力＋強制振動が加わる物体の力は次のとおりだね。
$$
\boldsymbol{F} =m\boldsymbol{g} -k\left[\boldsymbol{r}(t) - \boldsymbol{r}_{\rm box}(t)\right] -\gamma \boldsymbol{v}(t)
$$
力が与えられればニュートンの運動方程式（ \( \boldsymbol{a} = \boldsymbol{F} / m \) ）から、物体の加速度は次のように与えられるね。
$$
\boldsymbol{a}(t) =\boldsymbol{g} -\frac{k}{m}\,\boldsymbol{r}(t) -\frac{\gamma}{m}\, \boldsymbol{v}(t)
$$
Pythonによる計算結果
次のグラフは重力加速度$${\boldsymbol{g}=(0,0,-10)}$$、質量$${m=1}$$、ばね定数$${k=1}$$、粘性抵抗係数$${\gamma = 0.1}$$、初速度$${0}$$、初期位置を重力との釣り合いの位置（$${\boldsymbol{r}=(0,0,-10)}$$）として、角振動数を $${\omega=1.0\sim1.5}$$まで変化させたときの時間発展の様子だよ。粘性抵抗力が存在するので共鳴振動数に対応する角振動数（$${\omega=\sqrt{k/m}}$$）であっても振幅は無限に増幅するわけではなく一定値に収束するね。
プログラムソース（Python）
#################################################################
TITLE = "ばね弾性力＋重力＋粘性抵抗力＋強制振動による運動"
#################################################################
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 = 100.0 #計算終了時刻
DT = 0.01 #時間間隔
PLOT_FLAG = False #グラフ描画
ANIMATION_FLAG = True #アニメーション作成
ANIMATION_INTERVAL = 20 #アニメーション作成時の描画インターバル
###############################
# 計算グローバル変数
###############################
#時刻
times = np.arange(TIME_START, TIME_END, DT)
#重力加速度
vec_g = np.array([0, 0, -10])
#質量
m = 1.0
#ばね定数
k = 1
#粘性係数
gamma = 0.1
#角振動数
omegas = np.arange(1.0, 1.2, 0.01) * np.sqrt( k / m )
#初期位置
vec_r0 = vec_g / k
#初期速度
vec_v0 = np.array([0, 0, 0])
###############################
# 関数定義
###############################
#ばね弾性力＋重力の微分方程式
def equations(t, X , vec_g, m, gamma, k, omega ):
vec_box = np.array([ 0, 0, 1.0 * np.sin( omega * t) ])
vec_r = np.array([ X[0], X[1], X[2] ])
vec_v = np.array([ X[3], X[4], X[5] ])
vec_a = vec_g - k * (vec_r - vec_box) - gamma / m * vec_v
return [vec_v[0], vec_v[1], vec_v[2], vec_a[0], vec_a[1], vec_a[2]]
# omega 依存性を可視化
def calculation( omega ):
solution = integrate.solve_ivp(
equations, #常微分方程式
method = 'RK45', #計算方法
args = (vec_g, m, gamma, k, omega), #引数
t_span=( times[0], times[-1]), #積分区間
y0 = np.concatenate([vec_r0, vec_v0]), #初期値
t_eval = times #結果を取得する時刻
)
#計算結果
return solution.y[ 0 ], solution.y[ 1 ], solution.y[ 2 ]
##########################################################
if __name__ == "__main__":
#描画用リスト
xls = []
yls = []
#omega による違い
for omega in omegas:
x, y, z = calculation( omega )
xls.append( x )
yls.append( z )
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([-20.0, 0.2])
#グラフの描画
for i in range(len(omegas)):
c = colorsys.hls_to_rgb( i/len(omegas) , 0.5 , 1)
plt.plot(times, yls[i], linestyle='solid', linewidth = 2, color = c)
#グラフの表示
plt.show()
if( ANIMATION_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']
#アニメーション作成用
imgs = []
#アニメーション分割数
AN = len( times )
#軌跡リスト
r_trajectorys = [ 0.0 ] * len( omegas )
#アニメーションの各グラフを生成
for an in range( AN ):
t = DT * an
if( an % ANIMATION_INTERVAL != 0): continue
for i in range( len( omegas ) ):
c = colorsys.hls_to_rgb( i/len(omegas) , 0.5 , 1)
omega = omegas[ i ]
if(i == 0):
img = plt.plot([ omega ], [np.sin( omega * t )], color = c, marker = 's', markersize = 10)
else:
img += plt.plot([ omega ], [np.sin( omega * t )], color = c, marker = 's', markersize = 10)
img += plt.plot([omega], [yls[i][an]], color = c, marker = 'o', markersize = 15)
img += plt.plot([ omega, omega ], [np.sin( omega * t ), yls[i][an]], color = c, linestyle='solid', linewidth = 2 )
#時刻をグラフに追加
time = plt.text(26.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"$\omega \,[{\rm rad / 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([omegas[0], omegas[-1]])
plt.ylim([-20.0, 1.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()