2重平面振子

１. モデルの定式化

double pendulum

棒$${l_1}$$及び$${l_2}$$が鉛直となす角を$${\varphi_1,\varphi_2}$$とする.
座標原点を支点におくことにすると, 質点の座標はそれぞれ
$${(l_1\sin\varphi_1, -l_1\cos\varphi_1), \qquad (l_1\sin\varphi_1+l_2\sin\varphi_2, -l_1\cos\varphi_1-l_2\cos\varphi_2)}$$

質点$${m_1}$$に対する運動エネルギー$${T_1}$$,  及びポテンシャル・エネルギー$${U_1}$$は,
$${\displaystyle T_1 = \frac{1}{2}m_1(l_1\dot{\varphi_1})^2\quad,\qquad U_1 = -m_1gl_1\cos\varphi_1}$$

($${\displaystyle T_1=\frac{1}{2}m_1(\dot{x_1}^2+\dot{y_1}^2)}$$を計算しても,同じ結果になる.)

一方質点$${m_2}$$について, そのデカルト座標$${(x_2, y_2)}$$から,
$${\begin{array}{ll}x_2 &= l_1\sin\varphi_1 + l_2\sin\varphi_2 \quad,\quad y_2 = -(l_1\cos\varphi_1 + l_2\cos\varphi_2)\\\\\dot{x_2} &= l_1\dot{\varphi_1}\cos\varphi_1 + l_2\dot{\varphi_2}\cos\varphi_2 \quad,\quad \dot{y}_2 = l_1\dot{\varphi_1}\sin\varphi_1 + l_2\dot{\varphi_2}\sin\varphi_1\\\\\dot{x}_2^2 &= l_1^2\dot{\varphi_1}^2\cos^2\varphi_1 + 2l_1l_2\dot{\varphi_1}\dot{\varphi_2}\cos\varphi_1\cos\varphi_2 + l_2^2\dot{\varphi_2}^2\cos^2\varphi_2 \quad,\\\\\dot{y}_2^2 &= l_1^2\dot{\varphi_1}^2\sin^2\varphi_1 + 2l_1l_2\dot{\varphi_1}\dot{\varphi_2}\sin\varphi_1\sin\varphi_2 + l_2^2\dot{\varphi_2}^2\sin^2\varphi_2\end{array}}$$

$${\cos(\alpha-\beta)=\cos\alpha\cos\beta+\sin\alpha\sin\beta}$$, また $${\sin^2\alpha+\cos^2\alpha=1}$$ だったことを思い出して, 運動エネルギーは,
$${\displaystyle T_2 = \frac{1}{2}m_2(\dot{x_2}^2+\dot{y_2}^2)=\frac{m_2}{2}\left(l_1^2\dot{\varphi_1}^2+l_2^2\dot{\varphi_2}^2+2l_1l_2\cos(\varphi_1-\varphi_2)\dot{\varphi_1}\dot{\varphi_2}\right)}$$
また, ポテンシャル・エネルギーは,
$${U_2 = -m_2gl_1\cos\varphi_1 - m_2gl_2\cos\varphi_2}$$

従って, Lagrangian$${L=T_1+T_2-U_1-U_2}$$は,

$${\begin{array}{ll} L &= \displaystyle\frac{m_1+m_2}{2}l_1^2\dot{\varphi_1}^2 + \frac{m_2}{2}l_2^2\dot{\varphi_2}^2 + m_2l_1l_2\dot{\varphi_1}\dot{\varphi_2}\cos(\varphi_1-\varphi_2)  \\\\&\quad +(m_1+m_2)gl_1\cos\varphi_1 + m_2gl_2\cos\varphi_2\end{array}}$$

になる. 

$${\begin{array}{ll}\displaystyle\frac{\partial L}{\partial\varphi_1}&=-m_2l_1l_2\sin(\varphi_1-\varphi_2)\dot{\varphi_1}\dot{\varphi_2}-(m_1+m_2)gl_1\sin\varphi_1\quad,\\\\\displaystyle\frac{\partial L}{\partial\dot{\varphi_1}}&=(m_1+m_2)l_1^2\dot{\varphi_1}+m_2l_1l_2\cos(\varphi_1-\varphi_2)\dot{\varphi_2}\quad,\\\\\displaystyle\frac{\partial L}{\partial\varphi_2}&=m_2l_1l_2\sin(\varphi_1-\varphi_2)\dot{\varphi_1}\dot{\varphi_2}-m_2gl_2\sin\varphi_2\quad,\\\\\displaystyle\frac{\partial L}{\partial\dot{\varphi_2}}&=m_2l_2^2\dot{\varphi_2}+m_2l_1l_2\cos(\varphi_1-\varphi_2)\dot{\varphi_1}\end{array}}$$

より, Euler-Lagrange eq. $${\displaystyle \frac{\mathrm{d}}{dt}\frac{\partial L}{\partial \dot{\varphi_i}}-\frac{\partial L}{\partial \varphi_i}=0}$$ は,

$${\begin{array}{ll}(m_1+m_2)l_1^2\ddot{\varphi_1}+m_2l_1l_2\cos(\varphi_1-\varphi_2)\ddot{\varphi_2}\\\\\qquad +m_2l_1l_2\sin(\varphi_1-\varphi_2)\dot{\varphi_1}\dot{\varphi_2}+(m_1+m_2)gl_1\sin\varphi_1=0\quad,\\\\m_2l_2^2\ddot{\varphi_2}+m_2l_1l_2\cos(\varphi_1-\varphi_2)\ddot{\varphi_1}\\\\\qquad -m_2l_1l_2\sin(\varphi_1-\varphi_2)\dot{\varphi_1}\dot{\varphi_2} +m_2gl_2\sin\varphi_2=0\end{array}}$$

２. Pythonによる模擬実験

連立の常微分方程式にして, 関数odeを定義する.

$${\begin{cases} \displaystyle\frac{d\varphi_1}{dt}&=\dot{\varphi_1}\\\\\displaystyle\frac{d\varphi_2}{dt}&=\dot{\varphi_2}\\\\\displaystyle\frac{d\dot{\varphi_1}}{dt}&=\frac{Mg\sin\varphi_1+m_2\sin H(l_1\cos H+l_2)\dot{\varphi_1}\dot{\varphi_2}-m_2g\cos H\sin\varphi_2}{-\left(Ml_1-m_2l_1\cos^2H\right)}\\\\\displaystyle\frac{d\dot{\varphi_2}}{dt}&=\frac{Mg\cos H\sin\varphi_1+\sin H(Ml_1+m_2l_2\cos H)\dot{\varphi_1}\dot{\varphi_2}-Mg\sin\varphi_2}{Ml_2-m_2l_2\cos^2H} \end{cases}}$$
ここで, $${H=\varphi_1-\varphi_2 ,\quad M=m_1+m_2}$$である.

単振子との違いは質点が2個あるので, 関数ode の戻り値は,
単振子のときに$${\left[\varphi, \dot{\varphi}\right]}$$としていたのを,$${\left[\varphi_1, \dot{\varphi_1}, \varphi_2, \dot{\varphi_2}\right]}$$とする.

積分は, scipy.integrateのodeint を,
アニメーションは, matplotlib.animationのFuncAnimation を使う.

from numpy import sin, cos, pi
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np
from scipy.integrate import odeint
#
# Constants
#
G = 9.8            # [m/s^2]  acceleration of gravity
PHI1_0 = pi        # [rad] initial theta
PHI2_0 = -pi/6.    # [rad] initial theta
V1_0 = 0.          # [m/s]    initial velocity
V2_0 = 0.          # [m/s]    initial velocity
L1 = 1.            # [m]      length of pendulum
L2 = 1.            # [m]      length of pendulum
M1 = 1.            # [kg]     mass
M2 = 0.4           # [kg]     mass
DURATION = 15.     # [s]      duration time
INTERVAL = 0.05    # [s]      interval time
#
# Differential Equation
#
def ode(f, t):
    phi1, dphi1, phi2, dphi2 = f
    M, H = M1 + M2, phi1 - phi2
    dphi1dt = (M*G*sin(phi1)\
                + M2*(L2+L1*cos(H))*sin(H)*dphi1*dphi2\
                - M2*G*cos(H)*sin(phi2))\
                /(-M*L1+M2*L1*cos(H)*cos(H))
    dphi2dt = (M*G*cos(H)*sin(phi1)
                +(M2*L2*cos(H)+M*L1)*sin(H)*dphi1*dphi2\
                - M*G*sin(phi2))\
                /(M*L2-M2*L2*(cos(H))**2)
    return [dphi1, dphi1dt, dphi2, dphi2dt]
#
# Initial Condition
#
f_0 = [PHI1_0, V1_0/L1, PHI2_0, V2_0/L2]
t = np.arange(0, DURATION + INTERVAL, INTERVAL)
#
# Solve the Equation
#
sol = odeint(ode, f_0, t)
phi1, phi2 = sol[:, 0], sol[:, 2]
x1 = L1 * sin(phi1)
y1 = - L1 * cos(phi1)
x2 = x1 + L2 * sin(phi2)
y2 = y1 - L2 * cos(phi2)
#
# Prepare the Screen to display
#
fig = plt.figure()
ax = fig.add_subplot(111, aspect='equal', autoscale_on=False,
                     xlim = [-L1 - L2, L1 + L2], ylim = [-L1 - L2, L1 + L2])
ax.grid()
markers_on = [1,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, x1[i], x2[i]]
    next_y = [0, y1[i], y2[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 time 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
#ani.save('double_pendulum.mp4', fps=FPS, extra_args=['-vcodec', 'libx264'])
ani.save('double_pendulum.gif', writer='imagemagick', fps=FPS)