支点が円周上を周回している単振子

１．モデルの定式化

支点が鉛直平面内の円周上を一定の角速度$${\gamma}$$で一様に動いている単振子.
($${\gamma<0}$$でCW, $${\gamma>0}$$でCCW, $${\gamma}$$が$${0}$$なら支点は動かない)

円周上の点に注目して（ $${sin(-\gamma t)  = -sin(\gamma t)}$$ だから）CWとCCWの違い

右の図は,CCW方向に回転する場合のモデル化になっている.
質点$${m}$$の位置$${(x,y)}$$は,
$${\begin{cases}x&=a\cos\gamma t + l\sin\varphi\\y&=a\sin\gamma t - l\cos\varphi\end{cases}}$$

質点の速度,その自乗を作って,
$${\begin{aligned}\dot{x}&=-a\gamma\sin\gamma t + l\dot{\varphi}\cos\varphi\\\dot{y}&=a\gamma\cos\gamma t + l\dot{\varphi}\sin\varphi\end{aligned}}$$
$${\begin{aligned}\dot{x}^2&=a^2\gamma^2\sin^2\gamma t + l^2\dot{\varphi}^2\cos^2\varphi - 2al\gamma\dot{\varphi}\sin\gamma t\cos\varphi\\\dot{y}^2&=a^2\gamma^2\cos^2\gamma t + l^2\dot{\varphi}^2\sin^2\varphi + 2al\gamma\dot{\varphi}\sin\varphi\cos\gamma t\end{aligned}}$$

運動エネルギー$${T}$$とポテンシャル・エネルギー$${U}$$は,
$${\begin{aligned}T=\displaystyle\frac{m}{2}\left(\dot{x}^2+\dot{y}^2\right)&=\frac{m}{2}\{a^2\gamma^2 + l^2\dot{\varphi}^2 + 2al\gamma\dot{\varphi}\left(\sin\varphi\cos\gamma t - \cos\varphi\sin\gamma t\right)\}\\&=\frac{m}{2}\{a^2\gamma^2 + l^2\dot{\varphi}^2 + 2al\gamma\dot{\varphi}\sin(\varphi-\gamma t)\}\end{aligned}}$$
$${U=mgy=mga\sin\gamma t - mgl\cos\varphi}$$

Lagrangian $${L=T-U}$$は,
$${\begin{aligned}L&=T-U\\&=\displaystyle\frac{m}{2}a^2\gamma^2 + \frac{m}{2}l^2\dot{\varphi}^2 + alm\gamma\dot{\varphi}\sin(\varphi-\gamma t) - mga\sin\gamma t + mgl\cos\varphi\\&=\frac{m}{2}l^2\dot{\varphi}^2 + alm\gamma\dot{\varphi}\sin(\varphi-\gamma t) + mgl\cos\varphi +\left(\frac{m}{2}a^2\gamma^2 - mga\sin\gamma t\right)\end{aligned}}$$
右辺の最終項は, 次の関数$${f(t)}$$の時間に関する完全導関数になっている.
$${f(t)=\displaystyle\frac{m}{2}a^2\gamma^2t + mga\frac{1}{\gamma}\cos\gamma t,\quad\frac{df(t)}{dt}=\frac{m}{2}a^2\gamma^2-mga\sin\gamma t}$$
時間だけに依存する最終項を除けば, Lagrangianは最終的に次の様になる.
(＊座標と時間の任意の関数$${f}$$の,時間に関する完全導関数$${\dot{f}}$$をLagrangianの中に含んでいても,作用積分によってその変分は消えてしまう項になるので,取り除いても運動方程式としては変わらない.$${\rightarrow}$$ランダウ「力学」p.5)

$${\therefore L=\displaystyle\frac{m}{2}l^2\dot{\varphi}^2 + alm\gamma\dot{\varphi}\sin(\varphi-\gamma t)+mgl\cos\varphi}$$

次の計算をして,
$${\begin{aligned}\displaystyle\frac{\partial L}{\partial\dot{\varphi}}&=ml^2\dot{\varphi} + alm\gamma\sin(\varphi-\gamma t)\\\displaystyle\frac{\partial L}{\partial\varphi}&=alm\gamma\dot{\varphi}\cos(\varphi-\gamma t)- mgl\sin\varphi\end{aligned}}$$
もう少し準備して,
$${\displaystyle\frac{d}{dt}\frac{\partial L}{\partial\dot{\varphi}}=ml^2\ddot{\varphi}-alm\gamma^2\cos(\varphi-\gamma t)}$$

Eular-Langange eq. $${\left(\displaystyle\frac{d}{dt}\frac{\partial L}{\partial\dot{\varphi}}-\frac{\partial L}{\partial\varphi}=0\right)}$$ は,

$${\therefore \ddot{\varphi}=\displaystyle\frac{a\gamma^2}{l}\cos(\varphi-\gamma t)+\frac{a\gamma}{l}\dot{\varphi}\cos(\varphi-\gamma t)-\frac{g}{l}\sin\varphi}$$

一方, 左の図はCW方向に回転する場合のモデル化になっている.
(実際, 右の図のモデルで$${\gamma}$$を$${-\gamma}$$で置き換えても同じ事になる.)
質点$${m}$$の位置$${(x,y)}$$は,
$${\begin{cases}x&=a\cos\gamma t + l\sin\varphi\\y&=-a\sin\gamma t - l\cos\varphi\end{cases}}$$
質点の速度, その自乗と計算していくと,
$${\begin{aligned}\dot{x}&=-a\gamma\sin\gamma t + l\dot{\varphi}\cos\varphi\\\dot{y}&=-a\gamma\cos\gamma t + l\dot{\varphi}\sin\varphi\end{aligned}}$$
$${\begin{aligned}\dot{x}^2&=a^2\gamma^2\sin^2\gamma t + l^2\dot{\varphi}^2\cos^2\varphi - 2al\gamma\dot{\varphi}\sin\gamma t\cos\varphi\\\dot{y}^2&=a^2\gamma^2\cos^2\gamma t + l^2\dot{\varphi}^2\sin^2\varphi - 2al\gamma\dot{\varphi}\sin\varphi\cos\gamma t\end{aligned}}$$

運動エネルギー$${T}$$とポテンシャル・エネルギー$${U}$$は,
$${\begin{aligned}T=\displaystyle\frac{m}{2}\left(\dot{x}^2+\dot{y}^2\right)&=\frac{m}{2}\{a^2\gamma^2 + l^2\dot{\varphi}^2 - 2al\gamma\dot{\varphi}\left(\sin\varphi\cos\gamma t + \cos\varphi\sin\gamma t\right)\}\\&=\frac{m}{2}\{a^2\gamma^2 + l^2\dot{\varphi}^2 - 2al\gamma\dot{\varphi}\sin(\varphi+\gamma t)\}\end{aligned}}$$
$${U=mgy= - mga\sin\gamma t - mgl\cos\varphi}$$

Lagrangian $${L=T-U}$$は,
$${\begin{aligned}L&=T-U\\&=\displaystyle\frac{m}{2}a^2\gamma^2 + \frac{m}{2}l^2\dot{\varphi}^2 - alm\gamma\dot{\varphi}\sin(\varphi+\gamma t) + mga\sin\gamma t + mgl\cos\varphi\\&=\frac{m}{2}l^2\dot{\varphi}^2 - alm\gamma\dot{\varphi}\sin(\varphi+\gamma t) + mgl\cos\varphi +\left(\frac{m}{2}a^2\gamma^2 + mga\sin\gamma t\right)\end{aligned}}$$
右辺の最終項は, 次の関数$${f(t)}$$の時間に関する完全導関数になっている.
$${f(t)=\displaystyle\frac{m}{2}a^2\gamma^2t - mga\frac{1}{\gamma}\cos\gamma t,\quad\frac{df(t)}{dt}=\frac{m}{2}a^2\gamma^2+mga\sin\gamma t}$$
(＊座標と時間の任意の関数$${f}$$の,時間に関する完全導関数$${\dot{f}}$$をLagrangianの中に含んでいても,作用積分によってその変分は消えてしまう項になるので,取り除いても運動方程式としては変わらない.$${\rightarrow}$$ランダウ「力学」p.5)
時間だけに依存する最終項を除けば, Lagrangianは最終的に次の様になる.

$${\therefore L=\displaystyle\frac{m}{2}l^2\dot{\varphi}^2 - alm\gamma\dot{\varphi}\sin(\varphi+\gamma t)+mgl\cos\varphi}$$

次の計算をして,
$${\begin{aligned}\displaystyle\frac{\partial L}{\partial\dot{\varphi}}&=ml^2\dot{\varphi} - alm\gamma\sin(\varphi+\gamma t)\\\displaystyle\frac{\partial L}{\partial\varphi}&=-alm\gamma\dot{\varphi}\cos(\varphi+\gamma t)- mgl\sin\varphi\end{aligned}}$$
もう少し準備して,
$${\displaystyle\frac{d}{dt}\frac{\partial L}{\partial\dot{\varphi}}=ml^2\ddot{\varphi}-alm\gamma^2\cos(\varphi+\gamma t)}$$

Eular-Langange eq. $${\left(\displaystyle\frac{d}{dt}\frac{\partial L}{\partial\dot{\varphi}}-\frac{\partial L}{\partial\varphi}=0\right)}$$ は,

$${\therefore \ddot{\varphi}=\displaystyle\frac{a\gamma^2}{l}\cos(\varphi+\gamma t)-\frac{a\gamma}{l}\dot{\varphi}\cos(\varphi+\gamma t)-\frac{g}{l}\sin\varphi}$$

このCWの式は, CCWの場合のEular-Lgrange eq.の$${\gamma}$$を$${-\gamma}$$で置き換えて得られる式と同じ格好をしている.

２．Pythonによる模擬実験

連立の一階微分方程式に直して,関数ode定義する.
右の図（CCWモデル）の場合のEular-Lagrange eq.からは,
$${\begin{cases}\displaystyle\frac{\mathrm{d}\varphi}{\mathrm{d}t}&=\dot{\varphi}\\\\\displaystyle\frac{\mathrm{d}\dot{\varphi}}{\mathrm{d}t}&=\displaystyle\frac{a\gamma^2}{l}\cos(\varphi-\gamma t) +\frac{a\gamma}{l}\dot{\varphi}\cos(\varphi-\gamma t)- \frac{g}{l}\sin\varphi\end{cases}}$$
円周上の点の位置$${(x_0,y_0)}$$は、
$${x_0=a\cos\gamma t \quad,\quad y_0=a\sin\gamma t}$$
であり、質点の位置$${(x,y)}$$は、
$${x=x_0 + l\sin\varphi \qquad,\quad y=y_0 - l\cos\varphi}$$
である.

左の図（CWモデル）の場合のEular-Lagrange eq.からは,
$${\begin{cases}\displaystyle\frac{\mathrm{d}\varphi}{\mathrm{d}t}&=\dot{\varphi}\\\\\displaystyle\frac{\mathrm{d}\dot{\varphi}}{\mathrm{d}t}&=\displaystyle\frac{a\gamma^2}{l}\cos(\varphi+\gamma t) -\frac{a\gamma}{l}\dot{\varphi}\cos(\varphi+\gamma t)- \frac{g}{l}\sin\varphi\end{cases}}$$
円周上の点の位置$${(x_0,y_0)}$$は、
$${x_0=a\cos\gamma t \quad,\quad y_0=-a\sin\gamma t}$$
であり、質点の位置$${(x,y)}$$は、
$${x=x_0 + l\sin\varphi \qquad,\quad y=y_0 - l\cos\varphi}$$
である.

模擬実験によって, CWとCCWの動きを確認することができる.
わざわざ場合分けしてまでの記述をしているが, （当然の事ながら,） 一方の定式化だけにしておいて、定数GAMMAの符号を反転させても同じ現象になるよね.

import matplotlib.pyplot as plt
import numpy as np
from numpy import sin, cos, pi
from matplotlib.animation import FuncAnimation
from scipy.integrate import odeint
#
# Constants
#
CCW = False       # CW = not CCW
A = 0.8           # [m]       radius
F = 0.1           # [Hz]      frequency
GAMMA = 2.0*pi*F  # [rad/s]   angular velocity
G  = 9.8          # [m/s^2]   acceleration of gravity
PHI0 = pi/6.      # [rad]     initial angle
V0 = 0.05         # [m/s]     initial velocity
L = 1.            # [m]       length of the pendulum
DURATION = 20.    # [s]       duration time
INTERVAL = 0.05   # [s]       interval time
#
# Differential equation
#
def ode(f, t):
    phi, dphi = f
    if CCW:
        # CCW at the right figure
        dfdt = [dphi, (A*GAMMA*GAMMA/L)*cos(phi-GAMMA*t)\
                          + (A*GAMMA/L)*cos(phi-GAMMA*t)*dphi\
                          - (G/L)*sin(phi)]
    else:
        # CW at the left figure
        dfdt = [dphi, (A*GAMMA*GAMMA/L)*cos(phi+GAMMA*t)\
                          - (A*GAMMA/L)*cos(phi+GAMMA*t)*dphi\
                          - (G/L)*sin(phi)]
    return dfdt
#
# Initial condition
#
f0 = [PHI0, V0/L]    # [theta, v] at t = 0
t = np.arange(0, DURATION + INTERVAL, INTERVAL)
#
# Solve the equation
#
sol = odeint(ode, f0, t)
phi = sol[:, 0]
if CCW:
    # CCW at the right sided figure
    circ_x = A * cos(GAMMA*t)
    circ_y = A * sin(GAMMA*t)
    x = circ_x + L * sin(phi)
    y = circ_y - L * cos(phi)
else:
    # CW at the left sided figure
    circ_x = A * cos(GAMMA*t)
    circ_y = - A * sin(GAMMA*t)
    x = circ_x + L * sin(phi)
    y = circ_y - L * cos(phi)
#
# Prepare the screen to display
#
fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False,\
                     xlim=(-2*L, 2*L), ylim=(-2*L, 2*L))
ax.grid()
markers_on = [2]
line, = plt.plot([], [], 'ro-', markevery=markers_on,  animated = True)
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
#
# Simulate
#
def init():
    time_text.set_text('')
    return line, time_text

def update(i):
    next_x = [0, circ_x[i], x[i]]
    next_y = [0, circ_y[i], y[i]]
    line.set_data(next_x, next_y)
    time_text.set_text(time_template % (i*INTERVAL))
    return line, time_text

FRAME_INTERVAL = 1000 * INTERVAL # [msec] interval between frames
ani = FuncAnimation(fig, update, frames=np.arange(0, len(t)),
                    interval=FRAME_INTERVAL, init_func=init, blit=True)
#
# Show and save the results
#
plt.show()
FPS = 1000/FRAME_INTERVAL # frames per second
if CCW:
    fname = 'pendulum4CCW'
else:
    fname = 'pendulum4CW'
#ani.save(fname+'.mp4', fps=FPS, extra_args=['-vcodec', 'libx264'])
ani.save(fname+'.gif', writer='imagemagick', fps=FPS)