【Python】重力＋張力＋強制振動による強制振子運動
強制振子運動とは単振子の支点を強制的に振動させることで生じる運動だよ。単振動運動と同様、特定の角振動数（共鳴角振動数）で支点を振動させると、振り子の振幅が増幅されるね。ただし、単振動運動とは異なり、共鳴角振動数が振り子の振幅によって異なるので、単一の振動数で振り子が回転することはないね。下図は強制振動運動に関係するベクトル量の模式図だよ。
単振子運動とは異なり支点も空間に固定ではなく移動させることを考慮して、位置ベクトルを$${\boldsymbol{r}_ {\rm box}(t)}$$と表して、支点を基準とした相対的なおもりの位置ベクトルをひもベクトルと呼んで
$$
\boldsymbol{ L }(t) = \boldsymbol{r}(t) - \boldsymbol{r}_ {\rm box}(t)
$$
と定義するとするよ。このひもベクトルを用いると張力ベクトルは次のようになるね。
$$
\boldsymbol{ S } = -S\, \hat{\boldsymbol{ L }}(t)
$$
ただし、$${\hat{\boldsymbol{ L }}(t) = \boldsymbol{ L }(t)/L}$$だよ。伸び縮みしないひもを実現する条件は単振子のときと同様、次のように表すことができるね。
$$
|\boldsymbol{ L }(t) | = | \boldsymbol{r}(t) - \boldsymbol{r}_ {\rm box}(t)| =L
$$
続いて、単振子のときと同様、ひもベクトルの２乗を時間で微分すると、速度ベクトルに関する条件式が得られるね。
$$
\left[ \boldsymbol{r}(t) - \boldsymbol{r}_ {\rm box}(t)\right]\cdot \left[ \boldsymbol{v}(t) - \boldsymbol{v}_ {\rm box}(t)\right] =0
$$
ただし、$${\boldsymbol{v}_ {\rm box}(t)}$$は支点の速度ベクトルだよ。この関係式は相対位置ベクトルと相対速度ベクトルが直交することを表しているね。さらにこの式の両辺を時間で微分することで、加速度ベクトルに対する拘束条件が得られるね。
$$
\left[\boldsymbol{a}(t) - \boldsymbol{a}_ {\rm box}(t) \right]\cdot\left[\boldsymbol{r}(t) - \boldsymbol{r}_ {\rm box}(t) \right] + \left|\boldsymbol{v}(t) - \boldsymbol{v}_ {\rm box}(t)\right|^2=0
$$
ただし、$${\boldsymbol{a}_ {\rm box}(t)}$$は支点の加速度ベクトルだよ。速度に関する拘束条件と比べると複雑に見えるけれども、$${\boldsymbol{a}(t)}$$の項について解くと
$$
\boldsymbol{a}(t)\cdot\boldsymbol{L}(t) = \boldsymbol{a}_ {\rm box}(t)\cdot\boldsymbol{L}(t)- |\boldsymbol{v}(t) - \boldsymbol{v}_ {\rm box}(t)|^2
$$
おもりの加速度ベクトルに関する条件が得られるね（$${\boldsymbol{ L }(t) = \boldsymbol{r}(t) - \boldsymbol{r}_ {\rm box}(t)}$$）。
一方、張力と重力の合力で運動するおもりの加速度ベクトルはニュートンの運動方程式から
$$
\boldsymbol{a}(t) = \frac{1}{m}\left[\boldsymbol{S}(t)+\boldsymbol{F}_G \right]= \frac{1}{m}\left[-S \hat{\boldsymbol{L}}(t)+m\boldsymbol{g} \right]
$$
となるね。この右辺のうち$${S}$$が未知の量だけれども、さきに導出した加速度ベクトルの条件式と連立することで$${S}$$を決定することができるよ。具体的には上式の両辺に$${\boldsymbol{L}(t)}$$との内積をとると
$$
\boldsymbol{a}(t)\cdot\boldsymbol{L}(t) = \frac{1}{m}\left[-LS+m\boldsymbol{g}_g(t) \cdot\boldsymbol{L}(t) \right]
$$
となることから、先の加速度ベクトルの条件式に代入すると、
$$
S = \frac{m}{L} \left[ \left|\boldsymbol{v}(t) - \boldsymbol{v}_ {\rm box}(t)\right|^2 - \boldsymbol{a}_ {\rm box}(t)\cdot\boldsymbol{L}(t)+ \boldsymbol{g}\cdot\boldsymbol{L}(t) \right]
$$
と張力$${S}$$が得られるね。この式の右辺の支点の位置ベクトル、速度ベクトル、加速度ベクトルが既知で、おもりの位置ベクトルと速度ベクトルが得られている場合には、張力の大きさが得られることを意味するね。つまり、張力ベクトルは次のように得られるね。
ひもによる張力ベクトルの表式
$$
\boldsymbol{S}(t) =- \frac{m}{L^2} \left[ \left|\boldsymbol{v}(t) - \boldsymbol{v}_ {\rm box}(t)\right|^2 - \boldsymbol{a}_ {\rm box}(t)\cdot\boldsymbol{L}(t)+ \boldsymbol{g}\cdot\boldsymbol{L}(t) \right] \boldsymbol{ L }(t)
$$
さらにおもりに加わる力は次のようになるね。
振り子のおもりに加わる合力（張力＋重力）
$$
\boldsymbol{ F }(t) = \boldsymbol{ F }_G+\boldsymbol{ S }(t)= m\boldsymbol{ g } - \frac{m}{L^2} \left[ \left|\boldsymbol{v}(t) - \boldsymbol{v}_ {\rm box}(t)\right|^2 - \boldsymbol{a}_ {\rm box}(t)\cdot\boldsymbol{L}(t)+ \boldsymbol{g}\cdot\boldsymbol{L}(t) \right] \boldsymbol{ L }(t)
$$
最後にニュートンの運動方程式からおもりの加速度ベクトルは次のようになるね。
【計算アルゴリズム】振り子のおもりの加速度ベクトル（張力＋重力）
$$
\boldsymbol{ a }(t) = \boldsymbol{ g } - \frac{1}{L^2} \left[ \left|\boldsymbol{v}(t) - \boldsymbol{v}_ {\rm box}(t)\right|^2 - \boldsymbol{a}_ {\rm box}(t)\cdot\boldsymbol{L}(t)+ \boldsymbol{g}\cdot\boldsymbol{L}(t) \right] \boldsymbol{ L }(t)
$$
支点の位置ベクトル・速度ベクトル・加速度ベクトル
強制振子運動の例として支点の位置をx軸方向に振幅$${l}$$、角振動数$${\omega_{\rm box}}$$で振動させることを考えるよ。この場合、支点の位置ベクトル$${\boldsymbol{r}_{\rm box}(t) }$$、速度ベクトル$${\boldsymbol{v}_{\rm box}(t) }$$、加速度ベクトル$${\boldsymbol{a}_{\rm box}(t) }$$は次のとおりに与えられるね。
$$
\boldsymbol{r}_{\rm box}(t) = \boldsymbol{r}_{\rm box}(0) + l\sin(\omega_{\rm box} t)\, \hat{\boldsymbol{x}}
$$
$$
\boldsymbol{v}_{\rm box}(t) = l\omega_{\rm box} \cos(\omega_{\rm box} t) \, \hat{\boldsymbol{x}}
$$
$$
\boldsymbol{a}_{\rm box}(t) = -l\omega_{\rm box}^2 \sin(\omega_{\rm box} t) \, \hat{\boldsymbol{x}}
$$
計算結果
次の図は長さ（$${L}$$）10~12までのひもにつけられたおもりを、支点の角振動数$${\omega_{\rm box}=1}$$で強制振動させたときの位置（$${z}$$）の時間変化だよ。長さが10のときが一番揺さぶられて、長くなるほど、振幅が小さくなっているね。ちなみに角振動数$${\omega_{\rm box}=1 }$$は微小振動時の共鳴角振動数に対応するよ。
Pythonプログラムソース
次のプログラムソースは、張力＋重力＋強制振動による運動を計算できるよ。もし良かったら試してみてください。
#################################################################
TITLE = "張力＋重力＋強制振動による運動"
#################################################################
import os
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 = 200.0 #計算終了時刻
DT = 0.01 #時間間隔
PLOT_FLAG = True #グラフ描画
ANIMATION_FLAG = False #アニメーション作成
ANIMATION_INTERVAL = 20 #アニメーション作成時の描画インターバル
OUTPUT_FLAG = False
OUTPUT_PASS = "./data_20240925/"
if( OUTPUT_FLAG and os.path.isdir(OUTPUT_PASS) == False) : os.mkdir(OUTPUT_PASS)
###############################
# 計算グローバル変数
###############################
#時刻
times = np.arange(TIME_START, TIME_END, DT)
#重力加速度
vec_g = np.array([0, 0, -10])
#質量
m = 1.0
#粘性係数
gamma = 0.0
#ひもの長さ
L0s = np.arange(10, 13, 0.1)
#おもりの初期速度
vec_v0 = np.array([0, 0, 0])
#支点の振幅と角振動数
l = 1
omega_box = 1
###############################
# 関数定義
###############################
#支点の位置と速度
def box(t):
vec_r_box = np.array([ l * np.sin( omega_box * t ) , 0, 0])
vec_v_box = np.array([ l * omega_box * np.cos( omega_box * t ), 0, 0])
vec_a_box = np.array([ -l * omega_box**2 * np.sin( omega_box * t ), 0, 0])
return vec_r_box, vec_v_box, vec_a_box
#張力＋重力の微分方程式
def equations(t, X, L0, vec_g, m, compensationK, compensationGamma):
vec_r = np.array([ X[0], X[1], X[2] ])
vec_v = np.array([ X[3], X[4], X[5] ])
#支点の位置ベクトル・速度ベクトル・加速度ベクトル
vec_r_box, vec_v_box, vec_a_box = box( t )
#紐ベクトル
vec_L = vec_r - vec_r_box
#紐の長さ
L = np.sqrt( np.dot(vec_L, vec_L) )
#相対速度ベクトル
vec_v_v_box = vec_v - vec_v_box
#相対速度の大きさ
v_v_box = np.sqrt( np.dot(vec_v_v_box, vec_v_v_box) )
#張力
vec_S = - m / L**2 * ( v_v_box**2 - np.dot(vec_a_box, vec_L) + np.dot(vec_g, vec_r) ) * vec_L
#位置ベクトルと速度ベクトルの単位ベクトル
hat_L = vec_L / L
if( v_v_box > 0 ): hat_v = vec_v_v_box / v_v_box
else: hat_v = vec_v_v_box * 0
#補正力
f_c = - compensationK * ( L - L0 ) * hat_L - compensationGamma * np.dot(hat_v, hat_L) * hat_L
#質点に加わる力
F = vec_S + m * vec_g + f_c * np.sqrt( np.dot(vec_S, vec_S) )
#ニュートンの運動方程式
vec_a = F / m
return [vec_v[0], vec_v[1], vec_v[2], vec_a[0], vec_a[1], vec_a[2]]
# theta 依存性を可視化
def calculation(vec_r0, vec_v0, L0, compensationK, compensationGamma ):
solution = integrate.solve_ivp(
equations, #常微分方程式
method = 'RK45', #計算方法
args = (L0, vec_g, m, compensationK, compensationGamma), #引数
t_span=( times[0], times[-1]), #積分区間
rtol = 1.e-5, #相対誤差
atol = 1.e-5, #絶対誤差
y0 = np.concatenate([vec_r0, vec_v0]), #初期値
t_eval = times #結果を取得する時刻
)
#計算結果
return solution.y[ 0 ], solution.y[ 1 ], solution.y[ 2 ], solution.y[ 3 ], solution.y[ 4 ], solution.y[ 5 ]
##########################################################
if __name__ == "__main__":
#描画用リスト
xls = []
yls = []
zls = []
vxls = []
vyls = []
vzls = []
#theta による違い
for i in range(len(L0s)):
print( i )
L0 = L0s[ i ]
#初期位置と速度
vec_r0 = np.array([ 0, 0, -L0])
vec_v0 = np.array([ 0, 0, 0])
x, y, z, vx, vy, vz = calculation( vec_r0, vec_v0, L0, 0.3, 0.5 )
xls.append( x )
yls.append( y )
zls.append( z )
vxls.append( vx )
vyls.append( vy )
vzls.append( vz )
if( OUTPUT_FLAG ):
for i in range(len(times)):
path = OUTPUT_PASS + str( i ) + '.txt'
f = open(path, 'w')
for j in range( len(xls) ):
f.write( str(round( xls[ j ][ i ], 2)) + "\t" + str(round( yls[ j ][ i ], 2)) + "\t" + str(round( zls[ j ][ i ], 2)) + "\t" + str(round( vxls[ j ][ i ], 2)) + "\t" + str(round( vyls[ j ][ i ], 2)) + "\t" + str(round( vzls[ j ][ i ], 2)) + "\n")
f.close()
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([-13.0, 0])
#グラフの描画
for i in range(len(L0s)):
c = colorsys.hls_to_rgb( i/len(L0s) , 0.5 , 1)
plt.plot(times, zls[i], linestyle='solid', linewidth = 2, color = c)
#グラフの表示
plt.show()
if( ANIMATION_FLAG ):
#################################################################
# アニメーションの設定
#################################################################
fig = plt.figure( figsize=(10, 6.5) ) #画像サイズ
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( L0s )
#アニメーションの各グラフを生成
for an in range( AN ):
t = DT * an
vec_r_box, vec_v_box, vec_a_box = box( t )
if( an % ANIMATION_INTERVAL != 0): continue
for i in range( len( L0s ) ):
c = colorsys.hls_to_rgb( i/len(L0s) , 0.5 , 1)
if( i == 0 ):
img = plt.plot([xls[i][an]], [zls[i][an]], color = c, marker = 'o', markersize = 15)
else:
img += plt.plot([xls[i][an]], [zls[i][an]], color = c, marker = 'o', markersize = 15)
img += plt.plot([ vec_r_box[0], xls[i][an] ], [vec_r_box[2], zls[i][an]], color = c, linestyle='solid', linewidth = 2 )
#時刻をグラフに追加
time = plt.text(13.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"$x \,[{\rm m}]$", 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([-10, 10])
plt.ylim([-13.0, 0.1])
#アニメーションの生成
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()