2重振り子

(問題) 2つのおもりによって構成された2重振り子(Double Pendulum)がある。 これらは重力と逆向きの方向を$y$軸とする$xy$平面上で回転運動するものとする。 おもり1は質量$m_1$で、 原点$O$に固定された回転軸$O_1$から長さ$l_1$の糸でつるされている。 もうひとつのおもり2は質量$m_2$で、 おもり1に固定された回転軸$O_2$から長さ$l_2$の糸でつるされている。 この2重振り子のラグランジュ運動方程式を記述せよ。 また、計算機シミュレーションなどをしやすいように、 求められた2階常微分方程式を連立の1階常微分方程式に変換せよ。 また、この振り子が各々の角速度に比例する抵抗 や外力 を受けている場合においても同様に求めよ。

(解答) おもり1,2のy軸に対する角度を $\theta_1,\theta_2$とする。 おもり1の位置$G_1$の座標 $(x_{G_1},y_{G_1})$およびその速度$V_{G_1}$の2乗は、

$$x_{G_1} = l_1\sin\theta_1$$ (1)

$$y_{G_1} = -l_1\cos\theta_1$$ (2)

$$V_{G_1}^2 = \dot x_{G_1}^2 + \dot y_{G_1}^2$$

$$= (l_1\dot\theta_1\cos\theta_1)^2 + (l_1\dot\theta_1\sin\theta_1)^2$$

$$= l_1^2\dot\theta_1^2$$ (3)

となる。 おもり2の位置$G_2$の座標 $(x_{G_2},y_{G_2})$およびその速度$V_{G_2}$の2乗は、

$$x_{G_2} = l_1\sin\theta_1+l_2\sin\theta_2$$ (4)

$$y_{G_2} = -l_1\cos\theta_1-l_2\cos\theta_2$$ (5)

$$V_{G_2}^2 = \dot x_{G_2}^2 + \dot y_{G_2}^2$$

$$= (l_1\dot\theta_1\cos\theta_1+l_2\dot\theta_2\cos\theta_2)^2 + (l_1\dot\theta_1\sin\theta_1+l_2\dot\theta_2\sin\theta_2)^2$$

$$= l_1^2\dot\theta_1^2 + l_2^2\dot\theta_2^2 + 2l_1l_2\dot\theta_1\dot\theta_2\cos(\theta_1-\theta_2)$$ (6)

となる。 ここで、おもり1の運動エネルギー$T_1$および位置エネルギー$U_1$は、

$$T_1 = \frac12 m_1 V_{G_1}^2$$

$$= \frac12 m_1 l_1^2 \dot\theta_1^2$$ (7)

$$U_1 = m_1 g y_{G_1}$$

$$= -m_1 g l_1\cos\theta_1$$ (8)

また、おもり2の運動エネルギー$T_2$および位置エネルギー$U_2$は、

$$T_2 = \frac12 m_2 V_{G_2}^2$$

$$= \frac12 m_2 \left\{ l_1^2\dot\theta_1^2 + l_2^2\dot\theta_2^2 + 2l_1l_2\dot\theta_1\dot\theta_2\cos(\theta_1-\theta_2) \right\}$$ (9)

$$U_2 = m_2 g y_{G_2}$$

$$= -m_2 g \left( l_1\cos\theta_1 + l_2\cos\theta_2 \right)$$ (10)

となる。 この2重振り子のラグランジュ関数$L$は、

$$L = T - U$$

$$= (T_1+T_2) - (U_1+U_2)$$

$$= \frac12 m_1 l_1^2 \dot\theta_1^2 + \frac12 m_2 l_1^2\dot\theta_1^2 + \frac12 m_2 l_2^2\dot\theta_2^2$$

$$+ m_2 l_1l_2\dot\theta_1\dot\theta_2\cos(\theta_1-\theta_2) + m_1 g l_1\cos\theta_1 + m_2 g \left( l_1\cos\theta_1 + l_2\cos\theta_2 \right)$$

$$= \frac12 \left( m_1 + m_2 \right) l_1^2 \dot\theta_1^2 + \frac12 m_2 l_2^2 \dot\theta_2^2$$

$$+ m_2 l_1l_2\dot\theta_1\dot\theta_2\cos(\theta_1-\theta_2) + (m_1 + m_2) g l_1\cos\theta_1 + m_2 g l_2\cos\theta_2$$ (11)

となる。 このラグランジュ関数$L$から $\theta_1,\theta_2$ に関するラグランジュ運動方程式をたてるには、 以下の関係式を用いればよい。

$$\begin{eqnarray*} \frac{d}{dt} \frac{\partial L}{\partial\dot\theta_1} - \frac... ...rtial\dot\theta_2} - \frac{\partial L}{\partial\theta_2} &=& 0 \end{eqnarray*}$$

これを具体的に求めると、

$$\frac{\partial L}{\partial\dot\theta_1} = \left( m_1 + m_2 \right) l_1^2 \dot\theta_1 + m_2 l_1l_2\dot\theta_2\cos(\theta_1-\theta_2)$$

$$\Rightarrow \frac{d}{dt} \frac{\partial L}{\partial\dot\theta_1} = \left( m_1 + m_2 \right) l_1^2 \ddot\theta_1 + m_2 l_1l_2\left\{ ... ...a_2) - \dot\theta_2(\dot\theta_1-\dot\theta_2) \sin(\theta_1-\theta_2) \right\}$$ (12)

$$- \frac{\partial L}{\partial\theta_1} = m_2 l_1l_2\dot\theta_1\dot\theta_2\sin(\theta_1-\theta_2) + \left( m_1 + m_2 \right) g l_1\sin\theta_1$$ (13)

$$\frac{\partial L}{\partial\dot\theta_2} = m_2 l_2^2 \dot\theta_2 + m_2 l_1l_2\dot\theta_1\cos(\theta_1-\theta_2)$$

$$\Rightarrow \frac{d}{dt} \frac{\partial L}{\partial\dot\theta_2} = m_2 l_2^2 \ddot\theta_2 + m_2 l_1l_2\left\{ \ddot\theta_1\cos(\th... ...a_2) - \dot\theta_1(\dot\theta_1-\dot\theta_2) \sin(\theta_1-\theta_2) \right\}$$ (14)

$$- \frac{\partial L}{\partial\theta_2} = -m_2 l_1l_2\dot\theta_1\dot\theta_2\sin(\theta_1-\theta_2) + m_2 g l_2\sin\theta_2$$ (15)

ゆえに、 $\theta_1,\theta_2$に関するラグランジュ運動方程式は、

$$\left( m_1 + m_2 \right) l_1^2 \ddot\theta_1 + \left( m_1 + m_2 \... ...theta_2\cos(\theta_1-\theta_2) + \dot\theta_2^2\sin(\theta_1-\theta_2) \right\} = 0$$ (16)

$$m_2 l_2^2 \ddot\theta_2 + m_2 g l_2\sin\theta_2 + m_2 l_1l_2\left... ...theta_1\cos(\theta_1-\theta_2) - \dot\theta_1^2\sin(\theta_1-\theta_2) \right\} = 0$$ (17)

となる。$\Box$

ここで、求められた2階の微分方程式を連立1階微分方程式に直すため、 便宜上以下のようにおく。

$$\begin{eqnarray*} C &=& \cos(\theta_1-\theta_2) \\ S &=& \sin(\theta_1-\theta_2) \\ m_{1,2} &=& m_1 + m_2 \end{eqnarray*}$$

$$m_{1,2} l_1^2 \ddot\theta_1 + m_{1,2} g l_1\sin\theta_1 + m_2 l_1l_2 \left( C\ddot\theta_2 + S\dot\theta_2^2 \right) = 0$$ (18)

$$m_2 l_2^2 \ddot\theta_2 + m_2 g l_2\sin\theta_2 + m_2 l_1l_2 \left( C\ddot\theta_1 - S\dot\theta_1^2 \right) = 0$$ (19)

上式を整理すると,

$$\begin{eqnarray*} m_{1,2} l_1^2\ddot\theta_1 + m_2 l_1l_2 C\ddot\theta_2 &=& -m... ...\theta_1 &=& m_2 l_1l_2 S\dot\theta_1^2 - m_2 g l_2\sin\theta_2 \end{eqnarray*}$$

$$\begin{eqnarray*} m_{1,2} l_1^2 m_2 l_2^2\ddot\theta_1 + m_2 l_1l_2 C m_2 l_2^2... ..._1^2 l_2^2 CS\dot\theta_1^2 - m_2 l_1l_2 Cm_2 g l_2\sin\theta_2 \end{eqnarray*}$$

$$\Downarrow$$

$$( m_{1,2} l_1^2 m_2 l_2^2- m_2^2 l_1^2 l_2^2 C^2)\ddot\thet... ... l_1\sin\theta_1 m_2 l_2^2+ m_2 l_1l_2 Cm_2 g l_2\sin\theta_2$$

$$\begin{eqnarray*} m_{1,2} l_1^2 m_2 l_1l_2 C\ddot\theta_1 + m_2^2 l_1^2 l_2^2 C... ...m_2 l_1l_2 S\dot\theta_1^2 - m_{1,2} l_1^2m_2 g l_2\sin\theta_2 \end{eqnarray*}$$

$$\Downarrow$$

$$( m_2^2 l_1^2 l_2^2 C^2- m_{1,2} l_1^2 m_2 l_2^2)\ddot\thet... ...\sin\theta_1 m_2 l_1l_2 C+ m_{1,2} l_1^2m_2 g l_2\sin\theta_2$$

$$\ddot\theta_1 = \frac{-m_2 l_1l_2 S m_2 l_2^2\dot\theta_2^... ...sin\theta_2 }{ m_{1,2} l_1^2 m_2 l_2^2- m_2^2 l_1^2 l_2^2 C^2}$$ (20)

$$\ddot\theta_2 = \frac{-m_2^2 l_1^2 l_2^2 CS\dot\theta_2^2 ... ...sin\theta_2 }{ m_2^2 l_1^2 l_2^2 C^2- m_{1,2} l_1^2 m_2 l_2^2}$$ (21)

となる。さらに、

$$\vartheta_1 = \dot\theta_1, \qquad \vartheta_2 = \dot\theta_2$$

とおくことにより以下の連立1階微分方程式を得る。

$$\dot\theta_1 = \vartheta_1$$ (22)

$$\dot\theta_2 = \vartheta_2$$ (23)

$$\dot\vartheta_1 = \frac{-m_2 l_1l_2 S m_2 l_2^2\vartheta_2^2 - m_2^2 l_1^2 l_2^2 CS... ...2^2 g l_1l_2^2 C\sin\theta_2 }{ m_{1,2} l_1^2 m_2 l_2^2- m_2^2 l_1^2 l_2^2 C^2}$$ (24)

$$\dot\vartheta_2 = \frac{-m_2^2 l_1^2 l_2^2 CS\vartheta_2^2 - m_{1,2} l_1^2 m_2 l_1l... ...} l_1^2m_2 g l_2\sin\theta_2 }{ m_2^2 l_1^2 l_2^2 C^2- m_{1,2} l_1^2 m_2 l_2^2}$$ (25)

$\Box$

ちなみに、式(24), (25)は冗長であるので整理する ためにそれぞれ分母分子を $m_2^2 l_1 l_2^2$, $m_2^2 l_1^2 l_2$で割ると、

$$\begin{eqnarray*} \dot\vartheta_1 &=& \frac{g (\sin\theta_2 C - \frac{m_{1,2}}... ...a_1^2 + l_2 \vartheta_2^2 C) S}{l_2 (\frac{m_{1,2}}{m_2} - C^2)} \end{eqnarray*}$$

となる¹。 また、各々の角速度に比例する抵抗などの減衰 や強制振動のような外力 がある場合は、 $\dot\theta_1,\dot\theta_2$それぞれに対する減衰定数を $\lambda_1,\lambda_2$、外力を $\sigma_1,\sigma_2$ とすると、 $\theta_1,\theta_2$に関するラグランジュ運動方程式は 式(16),(17)から、

$$\left( m_1 + m_2 \right) l_1^2 \ddot\theta_1 + \left( m_1 + m_2 \... ...theta_2\cos(\theta_1-\theta_2) + \dot\theta_2^2\sin(\theta_1-\theta_2) \right\} = \sigma_1 - \lambda_1\dot\theta_1$$ (26)

$$m_2 l_2^2 \ddot\theta_2 + m_2 g l_2\sin\theta_2 + m_2 l_1l_2\left... ...theta_1\cos(\theta_1-\theta_2) - \dot\theta_1^2\sin(\theta_1-\theta_2) \right\} = \sigma_2 - \lambda_2\dot\theta_2$$ (27)

となる。 $\Box$

よって、式(22)-(25)に対応する減衰 や外力 がある場合の連立一階微分方程式は以下のようになる²。

$$\dot\theta_1 = \vartheta_1$$ (28)

$$\dot\theta_2 = \vartheta_2$$ (29)

$$\dot\vartheta_1 = \frac{\displaystyle{\left\{ \begin{array}{l} \displaystyle{ -m_2 ... ...rray}\right\}}}{\displaystyle{ m_{1,2} l_1^2 m_2 l_2^2- m_2^2 l_1^2 l_2^2 C^2}}$$ (30)

$$\dot\vartheta_2 = \frac{\displaystyle{\left\{ \begin{array}{l} \displaystyle{ -m_2^... ...rray}\right\}}}{\displaystyle{ m_2^2 l_1^2 l_2^2 C^2- m_{1,2} l_1^2 m_2 l_2^2}}$$ (31)

$\Box$ ［補足］ ハミルトンの正準方程式の導出

ここで、運動エネルギー$T$を以下のように表す。

$$T = \frac12 \dot{\bm\theta}^\top {R} \dot{\bm\theta}$$ (32)

但しここで、 $\bm\theta\in\mathbb{R}^2$、

$${R} = \left(\begin{array}{cc} m_{1,2} l_1^2& m_2 l_1l_2 C\\ m_2 l_1l_2 C& m_2 l_2^2 \end{array} \right)$$ (33)

である。すると、角運動量$\bm{p}$は、 と表される。ハミルトニアン$H$は

$$H = \frac12 \bm{p}^\top {R}^{-1} \bm{p} + U$$ (34)

であり、ハミルトニアン$H$からハミルトンの正準方程式を得るには以下

$$\dot{\bm\theta} = \frac{\partial H}{\partial \bm{p}} = \frac{\partial}{\partial \bm{p}} \left( \frac12 \bm{p}^\top {R}^{-1} \bm{p} + U \right)$$

$$= {R}^{-1} \bm{p}$$ (35)

$$\dot{\bm{p}} = - \frac{\partial H}{\partial \bm\theta} + \bm{F} = - \frac{\parti... ...tial \bm\theta} \left( \frac12 \bm{p}^\top {R}^{-1} \bm{p} + U \right) + \bm{F}$$

$$= - \frac12 \frac{\partial}{\partial \bm\theta} \left( \bm{p}^\top {R}^{-1} \bm{p} \right) - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= \frac12 \frac{\partial}{\partial \bm\theta} \left( \dot{\bm\theta... ...op {R} \dot{\bm\theta} \right) - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$ (36)

の関係式を用いればよい³。 ちなみに$\bm{F}$は、前述の抵抗や外力がある場合は、

$$\bm{F} = \left[\begin{array}{c} \sigma_1 - \lambda_1\dot\... ...1 \\ \sigma_2 - \lambda_2\dot\theta_2 \end{array} \right]$$ (37)

さもなくば${F}_d=0$($d=1,2$)である。

具体的に求めると、ハミルトンの正準方程式は以下のようになる。

$$\dot{\bm\theta} = {R}^{-1} \bm{p}$$

$$= \frac1{ m_{1,2} l_1^2 m_2 l_2^2- m_2^2 l_1^2 l_2^2 C^2}\left[\beg... ..._2 l_1l_2 C{p}_2\\ - m_2 l_1l_2 C{p}_1 + m_{1,2} l_1^2{p}_2 \end{array}\right]$$ (38)

$$\dot{\bm{p}} = \frac12 \frac{\partial}{\partial \bm\theta} \left( \dot{\bm\theta... ...bm\theta}, \ldots \right)^\top - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= \left[\begin{array}{c} -m_2 l_1l_2 S\dot\theta_1\dot\theta_2\\ m... ...m_{1,2} g l_1\sin\theta_1 \\ m_2 g l_2\sin\theta_2 \end{array}\right] + \bm{F}$$ (39)

$\Box$ ［発展］ 動座標系におけるラグランジュ運動方程式及びハミルトンの正準方程式の導出

$O_1$が時変の位置 $\bm{\Gamma}=(\Gamma_x, \Gamma_y)$によって動くことに より、振り子が強制振動される場合を考える⁴。 $i$番目の質点の位置ベクトルを $\bm{r}_i=(x_{G_i},y_{G_i})$として、この 動座標系における運動エネルギー$T$は、

$$T = \frac12\sum_{i=1}^n m_i\vert\dot{\bm{r}}_i\vert^2 + ... ...bm{r}}_i + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2$$ (43)

となり、位置エネルギー$U$は、

$$U = \sum_{i=1}^n m_i g y_{G_i} + \sum_{i=1}^n m_i g \Gamma_y$$ (44)

となる。ここで、

$$\sum_{i=1}^n m_i\dot{\bm\Gamma}^\top\dot{\bm{r}}_i = m_1 (\dot\Gamma_x, \dot\Gamma_y)\left[\begin{array}{c} l_1 \dot\t... ... l_1 \dot\theta_1\sin\theta_1 + l_2 \dot\theta_2\sin\theta_2 \end{array}\right]$$

$$= \left[\begin{array}{c} m_{1,2} l_1 \left(\dot\Gamma_x\cos\theta_1... ... \dot\theta_1\\ \dot\theta_2 \end{array}\right] = \bm{p}_O^\top\dot{\bm\theta}$$ (45)

とおくと、ラグランジュ関数$L$は$L=T-U$より、

$$L = \frac12 \dot{\bm\theta}^\top {R} \dot{\bm\theta} + \bm{p}_O^\top\dot{\bm\theta} + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2 - U$$ (46)

となり、ラグランジュ運動方程式は、

$${R} \ddot{\bm\theta} + \dot{R} \dot{\bm\theta} + \frac{\p... ...heta}\right) + \frac{\partial U}{\partial\bm\theta} = \bm{F}$$ (47)

となる。但しここで、

$$\frac{\partial \bm{p}_O}{\partial t} = \left[\begin{array... ...s\theta_2+\ddot\Gamma_y\sin\theta_2\right) \end{array}\right]$$ (48)

である。 この場合の角運動量$\bm{p}$は、 であるので、 より、ハミルトニアン$H$は、

$$H = \frac12 \left( \bm{p} - \bm{p}_O \right)^\top {R}^{-1} ... ...ight) + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2 + U$$ (49)

であり、ハミルトニアン$H$からハミルトンの正準方程式を得るには以下

$$\dot{\bm\theta} = \frac{\partial H}{\partial \bm{p}} = \frac{\partial}{\partial \bm... ...{p}_O \right) + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2 + U \right\}$$

$$= {R}^{-1} \left( \bm{p} - \bm{p}_O \right)$$ (50)

$$\dot{\bm{p}} = - \frac{\partial H}{\partial \bm\theta} + \bm{F} = - \frac{\parti... ...ght) + \frac12\sum_{i=1}^n m_i\vert\dot{\bm\Gamma}\vert^2 + U \right\} + \bm{F}$$

$$= - \frac12 \frac{\partial}{\partial \bm\theta} \left\{ \left( \bm{... ...p} - \bm{p}_O \right) \right\} - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= - \left( \frac12 \left( \bm{p} - \bm{p}_O \right)^\top \frac{\par... ..._O \right) ,\ldots\right)^\top - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= - \left( - \frac12 \dot{\bm\theta}^\top \frac{\partial {R}}{\part... ...\bm\theta} ,\ldots\right)^\top - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$

$$= \frac12 \frac{\partial}{\partial \bm\theta} \left( \dot{\bm\theta... ...O^\top \dot{\bm\theta} \right) - \frac{\partial U}{\partial \bm\theta} + \bm{F}$$ (51)

の関係式を用いればよい。但しここで、

$$\frac{\partial}{\partial \bm\theta}\left( \bm{p}_O^\top \dot... ...ot\Gamma_y\cos\theta_2\right) \dot\theta_2 \end{array}\right]$$ (52)

である。

... となる¹

これは明らかに文献[5]の式(B5), (B6)、 文献[12]の式(1.3)、文献[13]の式(1.12)と等し い。文献[9]の式(4C-4)、文献[11] p.30 の式とは異なっており、おそらくこれらは誤植であると思われる。

... がある場合の連立一階微分方程式は以下のようになる²

ソース で1箇所符合が反転しており、その結果、式(20), (24), (30)の一部符合が間違っておりました。 お詫びして訂正致します(2003年10月10日)。

... の関係式を用いればよい³

この導出には、 の関係から、

$$\frac{\partial}{\partial \theta_i} \left( {R} \cdot {R}^{-1}... ...{R}^{-1} + {R} \frac{\partial {R}^{-1}}{\partial \theta_i} = 0$$

であるので、これに $\bm{p}^\top {R}^{-1}$を左から、$\bm{p}$を右からかけて、

$$\bm{p}^\top {R}^{-1} \frac{\partial {R}}{\partial \theta_i} ... ...R}^{-1} {R} \frac{\partial {R}^{-1}}{\partial \theta_i} \bm{p}$$

ゆえに、

$$\dot{\bm\theta}^\top \frac{\partial {R}}{\partial \theta_i} ... ...\bm{p}^\top \frac{\partial {R}^{-1}}{\partial \theta_i} \bm{p}$$

が成立することを利用している[14]。

... より、振り子が強制振動される場合を考える⁴

以下は強制振動の一例 とその速度、加速度。

$$\bm{\Gamma} = \left(\Delta_x\sin(\Omega_x t), \Delta_y\cos(\Omega_y t)\right)$$ (40)

$$\dot{\bm\Gamma} = \left(\Delta_x\Omega_x\cos(\Omega_x t), -\Delta_y\Omega_y\sin(\Omega_y t)\right)$$ (41)

$$\ddot{\bm\Gamma} = \left(-\Delta_x\Omega_x^2\sin(\Omega_x t), -\Delta_y\Omega_y^2\cos(\Omega_y t)\right)$$ (42)

このサイトについてのご意見ご要望は taiji@aihara.co.jp までお願いします。 Copyright (C) 1996-1998,2003-2004 Taiji Yamada, AIHARA Electrical Engineering Co., Ltd.
All rights reserved.