【Python】ばね弾性力＋重力＋粘性抵抗力による運動（減衰振動運動）
今日も古典力学の復習を進めていくよ。今回は前回の「ばね弾性力＋重力」に加えて粘性抵抗力が加わった場合の振動運動だよ。重力＋ばね弾性力＋粘性抵抗力が加わる物体の力は次のとおりだね。
$$
\boldsymbol{F} =m\boldsymbol{g} -k\boldsymbol{r}(t) -\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}$$、初速度$${0}$$、初期位置を原点として、粘性抵抗係数を $${\gamma=0.0\sim0.1}$$まで変化させたときの時間発展の様子だよ。運動の様子は包絡関数と三角関数の積で表すことのできる減衰振動運動となるね。
プログラムソース（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 = True #グラフ描画
ANIMATION_FLAG = False #アニメーション作成
ANIMATION_INTERVAL = 20 #アニメーション作成時の描画インターバル
###############################
# 計算グローバル変数
###############################
#時刻
times = np.arange(TIME_START, TIME_END, DT)
#重力加速度
vec_g = np.array([0, 0, -10])
#質量
m = 1.0
#ばね定数
k = 1
#粘性係数
gammas = np.arange(0.0, 0.1, 0.001)
#初期位置
vec_r0 = np.array([0, 0, 0])
#初期速度
vec_v0 = np.array([0, 0, 0])
###############################
# 関数定義
###############################
#ばね弾性力＋重力の微分方程式
def equations(t, X , vec_g, m, gamma, k):
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 - gamma / m * vec_v
return [vec_v[0], vec_v[1], vec_v[2], vec_a[0], vec_a[1], vec_a[2]]
# k 依存性を可視化
def calculation( gamma ):
solution = integrate.solve_ivp(
equations, #常微分方程式
method = 'RK45', #計算方法
args = (vec_g, m, gamma, k), #引数
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 = []
#gamma による違い
for gamma in gammas:
x, y, z = calculation( gamma )
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(gammas)):
c = colorsys.hls_to_rgb( i/len(gammas) , 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( gammas )
#アニメーションの各グラフを生成
for an in range( AN ):
t = DT * an
if( an % ANIMATION_INTERVAL != 0): continue
for i in range( len( gammas ) ):
c = colorsys.hls_to_rgb( i/len(gammas) , 0.5 , 1)
gamma = gammas[ i ]
if( i == 0 ):
img = plt.plot([gamma], [yls[i][an]], color = c, marker = 'o', markersize = 15)
else:
img += plt.plot([gamma], [yls[i][an]], color = c, marker = 'o', markersize = 15)
img += plt.plot([ gamma, gamma ], [0, 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"$\gamma \,[{\rm Pa \cdot 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.0, 0.1])
plt.ylim([-20.0, 0.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()