【Python】張力＋重力による単振子運動
単振り子とは伸び縮みしないひもに結びつけたおもりによる振り子だよ。ばねの長さに比例した復元力を与えれば良い単振動と違って、伸び縮みしないひもに加わる張力の計算方法は簡単じゃないね。まずは単振り子に関係する物理量を図示するね。
おもりには重力の他に張力が加わるのだけれども、この張力は先にも述べたけれどもおもりの位置だけでは決めることはできないね。支点の位置を原点としておもりの位置ベクトルを$${\boldsymbol{r}(t)}$$とすると、張力$${\boldsymbol{S}}$$は
$$
\boldsymbol{S} = - S \hat{\boldsymbol{r}}(t)
$$
と表すことができるね。ただし、$${\hat{\boldsymbol{r}} = \boldsymbol{r}/|{\boldsymbol{r}}| }$$だよ。重力を含めておもりに加わる力は次のとおりだね。
$$
\boldsymbol{F} = \boldsymbol{S}+\boldsymbol{F}_G = - S\, \hat{\boldsymbol{r}}(t) + m\boldsymbol{g}
$$
$${S}$$はこのおもりの運動に関する拘束条件で決定することができるので導出するよ。拘束条件はひもの長さがいつでも一定という条件
$$
|\boldsymbol{r}(t)| =L
$$
だよ。これ拘束条件は単に長さに関する条件のみを含んでいるわけではなく、位置ベクトルの内積
$$
\boldsymbol{r}(t)\cdot\boldsymbol{r}(t)=L^2
$$
を取っておいて、両辺を時間で微分すると
$$
\boldsymbol{v}(t)\cdot\boldsymbol{r}(t)=0
$$
が得られるね。この条件は位置ベクトルと速度ベクトルが必ず直交することを表しているよ。さらにもう一度時間で微分すると
$$
\boldsymbol{a}(t)\cdot\boldsymbol{r}(t) + |\boldsymbol{v}(t)|^2 =0
$$
が得られるね。これはおもりの加速度に関する拘束条件にもなっているね。一方、ニュートンの運動方程式から加速度ベクトルは
$$
\boldsymbol{a}(t) = \frac{\boldsymbol{F}}{m} =-\frac{S}{m}\, \hat{\boldsymbol{r}}(t) + \boldsymbol{g}
$$
で与えられるので、この両辺に位置ベクトルをとると
$$
\boldsymbol{a}(t)\cdot \boldsymbol{r}(t) = -\frac{S L }{m}+\boldsymbol{g}\cdot\boldsymbol{r}(t)
$$
が得られるので、さきの加速度ベクトルに対する拘束条件と比較すると、張力の大きさは
$$
S = \frac{m}{L} \left[ |\boldsymbol{v}(t)|^2 +\boldsymbol{g}\cdot\boldsymbol{r}(t) \right]
$$
となるね。つまり、張力ベクトルは
$$
\boldsymbol{S} =- \frac{m}{L^2} \left[|\boldsymbol{v}(t)|^2 + \boldsymbol{g}\cdot\boldsymbol{r}(t) \right] \boldsymbol{r}(t)
$$
となるね。
Pythonによる計算結果
次のグラフは重力加速度$${\boldsymbol{g}=(0,0,-10}}$$、質量$${m=1}$$、初速度$${0}$$として、ひもの長さを0.01から1まで変化させたときのおもりの位置（$${z}$$成分）の時間変化だよ。初期位置としておもりが丁度水平になるように配置しているよ。ひもの長さが長いほど周期も長くなる様子がわかるね。
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 = 10.0 #計算終了時刻
DT = 0.01 #時間間隔
PLOT_FLAG = True #グラフ描画
ANIMATION_FLAG = False #アニメーション作成
ANIMATION_INTERVAL = 10 #アニメーション作成時の描画インターバル
###############################
# 計算グローバル変数
###############################
#時刻
times = np.arange(TIME_START, TIME_END, DT)
#重力加速度
vec_g = np.array([0, 0, -10])
#質量
m = 1.0
#粘性係数
gamma = 0.0
#ひもの長さ
L0s = np.arange(0.01, 1+0.05, 0.05)
#初期速度
vec_v0 = np.array([0, 0, 0])
#初期角度
theta = np.pi / 2.0
###############################
# 関数定義
###############################
#張力＋重力の微分方程式
def equations(t, X , vec_g, m):
vec_r = np.array([ X[0], X[1], X[2] ])
vec_v = np.array([ X[3], X[4], X[5] ])
#紐の長さ
L = np.sqrt( np.dot(vec_r, vec_r) )
#速度の大きさ
v = np.sqrt( np.dot(vec_v, vec_v) )
#張力
vec_S = - m / L**2 * ( v**2 + np.dot(vec_g, vec_r) ) * vec_r
#質点に加わる力
F = vec_S + m * vec_g
#ニュートンの運動方程式
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( L0 ):
#初期位置
vec_r0 = np.array([ L0 * np.sin( theta ), 0, L0 * np.cos( theta )])
solution = integrate.solve_ivp(
equations, #常微分方程式
method = 'RK45', #計算方法
args = (vec_g, m), #引数
t_span=( times[0], times[-1]), #積分区間
rtol = 1.e-10, #相対誤差
atol = 1.e-10, #絶対誤差
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 = []
#theta による違い
for L0 in L0s:
x, y, z = calculation( L0 )
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([-1.0, 0])
#グラフの描画
for i in range(len(L0s)):
c = colorsys.hls_to_rgb( i/len(L0s) , 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( L0s )
#アニメーションの各グラフを生成
for an in range( AN ):
t = DT * an
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]], [yls[i][an]], color = c, marker = 'o', markersize = 15)
else:
img += plt.plot([xls[i][an]], [yls[i][an]], color = c, marker = 'o', markersize = 15)
img += plt.plot([ 0, xls[i][an] ], [0, yls[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([-1, 1])
plt.ylim([-1, 0])
#アニメーションの生成
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()