【Python】重力＋空気抵抗力による運動
地球上で運動する物体には重力のほかにも空気抵抗力が加わるね。空気抵抗力は速度$${ \boldsymbol{v}(t) }$$）に比例する粘性抵抗力（粘性抵抗係数 $${ \gamma }$$）や速度の２乗に比例する慣性抵抗力（慣性抵抗係数 $${ \beta }$$）が一般的に用いられるね。式で表すと次のとおりだよ。
$$
\boldsymbol{F} =m\boldsymbol{g} -\gamma\boldsymbol{v}(t) - \beta|\boldsymbol{v}(t)|\boldsymbol{v}(t)
$$
２項目が粘性抵抗力、３項目が慣性抵抗力で符号がマイナスなのは、速度方向と反対向きだからだね。$${ \gamma }$$ と $${ \beta }$$ は物体の形状や、媒質（空気の種類や圧力など）の性質で決定される値だよ。力が与えられればニュートンの運動方程式（ $${ \boldsymbol{a} = \boldsymbol{F} / m }$$ ）から、物体の加速度は次のように与えられるね。
$$
\boldsymbol{a}(t) =\boldsymbol{g} -\frac{\gamma}{m}\boldsymbol{v}(t) - \frac{\beta}{m}|\boldsymbol{v}(t)|\boldsymbol{v}(t)
$$
Pythonによる計算結果
初速度と投射角度を一定（（$${ v_0 = 10, \theta=45^\circle }$$））として、粘性抵抗係数を変化（$${ \gamma = 0 \sim 4.0}$$ ）させたときの結果だよ。粘性抵抗係数が大きくなるほど放物線が崩れて投射距離が短くなっていくね。
次のグラフは初速度を一定とした水平投射だよ。慣性抵抗係数をを変化（$${ \gamma = 0.1 \sim 5.0}$$ ）させたときの結果だよ。慣性抵抗力が大きいほど、落下速度が遅くなるね。
プログラムソース（Python）
次のプログラムソースは、Scipyの常微分方程式を解く integrate.solve_ivp を用いて運動方程式を計算しているよ。グラフ描画に加えて、アニメーションも出力できるよ。ぜひ試してみてください。
#################################################################
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 = 2.4 #計算終了時刻
DT = 0.01 #時間間隔
PLOT_FLAG = True #グラフ描画
ANIMATION_FLAG = True #アニメーション作成
ANIMATION_INTERVAL = 4 #アニメーション作成時の描画インターバル
###############################
# 計算グローバル変数
###############################
#時刻
times = np.arange(TIME_START, TIME_END, DT)
#重力加速度
vec_g = np.array([0, 0, -10])
#質量
m = 1.0
#粘性係数
gammas = np.arange(0.0, 4.0, 0.005)
#慣性係数
beta = 0
#投射角度
theta = np.pi / 4
#初期位置
vec_r0 = np.array([0, 0, 0])
#初速度の大きさ
v0 = 10
#初期速度
vec_v0 = [ v0 * np.cos( theta ), 0.0, v0 * np.sin( theta ) ]
###############################
# 関数定義
###############################
#重力+空気抵抗力の微分方程式
def equations(t, X , vec_g, m, gamma, beta):
vec_r = np.array([ X[0], X[1], X[2] ])
vec_v = np.array([ X[3], X[4], X[5] ])
v = np.sqrt( np.sum( np.abs( vec_v**2 ) ) )
vec_a = vec_g - gamma / m * vec_v - beta / m * v * vec_v
return [vec_v[0], vec_v[1], vec_v[2], vec_a[0], vec_a[1], vec_a[2]]
# gamma 依存性を可視化
def calculation( gamma ):
#ルンゲクッタで解く
solution = integrate.solve_ivp(
equations, #常微分方程式
method = 'RK45', #計算方法
args = (vec_g, m, gamma, beta), #引数
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.08, 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"$x \,[{\rm m}]$", fontname="Yu Gothic", fontsize=20, fontweight=500 )
plt.ylabel( r"$h \,[{\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, 10])
plt.ylim([0, 2.5])
#グラフの描画
for i in range(len(gammas)):
c = colorsys.hls_to_rgb( i/len(gammas) , 0.5 , 1)
plt.plot(xls[i], yls[i], linestyle='solid', linewidth = 2, color=c)
#グラフの表示
plt.show()
if( ANIMATION_FLAG ):
#################################################################
# アニメーションの設定
#################################################################
fig = plt.figure( figsize=(12, 8) ) #画像サイズ
plt.subplots_adjust(left = 0.08, 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' #数式用フォント
#アニメーション作成用
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)
if( i == 0 ):
img = plt.plot([xls[i][an]], [yls[i][an]], color = c, marker = 'o', markersize = 10)
else:
img += plt.plot([xls[i][an]], [yls[i][an]], color = c, marker = 'o', markersize = 10)
#軌跡描画用リスト
xl = []
yl = []
for tn in range(len(times)):
if( tn >= an ): break
xl.append(xls[i][tn])
yl.append(yls[i][tn])
img += plt.plot( xl, yl, color = c, linestyle='solid', linewidth = 1)
#時刻をグラフに追加
time = plt.text(8.8, 5.05, "t = " + str("{:.2f}".format(t)), fontsize=18)
img.append( time )
imgs.append( img )
#グラフタイトルと軸ラベルの設定
plt.title( TITLE, fontname="Yu Gothic", fontsize=20, fontweight=1000 )
plt.xlabel( r"$x \,[{\rm m}]$", fontname="Yu Gothic", fontsize=20, fontweight=500 )
plt.ylabel( r"$h \,[{\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, 10.0])
plt.ylim([0, 2.5])
#アニメーションの生成
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()これからも応援よろしくお願いたします！
この記事が気に入ったらサポートをしてみませんか？