二重振り子のシミュレーション


二重振り子を例にカオス理論について解説します(オタ芸ではありません。似てるけど)

H﨑　 2018/6/28　 力学, 二重振り子

二重振り子とは？？

普通の振り子の先にもう一つ振り子をつないだだけの簡単な構造なのですが、これがまた興味深い。いわゆるカオスな運動をします。

カオス理論は、力学系の一部に見られる、数的誤差により予測できないとされている複雑な様子を示す現象を扱う理論である。ここで言う予測できないとは、決してランダムということではない。その振る舞いは決定論的法則に従うものの、初期値鋭敏性ゆえに、数値解析の過程での誤差によっても、得られる値と真の値とのずれが増幅される。そのため予測が事実上不可能という意味である。(Wikipediaより引用)

もう少し噛み砕いて説明すると、振り子の動きを運動方程式として書くことができる、つまり振り子はある法則に従って動いているのですが、その運動方程式を使って二重振り子の動きを予測することは不可能であると言うことです。

どれだけ高性能なスパコンを使っても、数十秒後には誤差が指数関数的に増加して、バラバラの動きとなってしまいます。

そのとっても無意味な運動方程式を立ててやる！！というのが今回のお題です。

どうでもいいですけど、二重振り子の動きとオタ芸の腕の動きってなんか似てません？はい、どうでもいいですね。

ひたすら計算

Drawing

二つのおもりで構成される二重振り子について考えてみたいと思います。

おもり1は質量$m_1$で紐の長さは$l_1$で角度を$\theta_1$、座標を$(x,y)$
おもり2は質量$m_2$で紐の長さは$l_2$で角度を$\theta_2$、座標を$(X,Y)$とします。
そうすると

$$x = l_1\sin\theta_1　,　y = -l_1\cos\theta_1$$ $$X = l_1\sin\theta_1+l_2\sin\theta_2　,　Y = -l_1\cos\theta_1-l_2\cos\theta_2$$

各座標の速度は時間tで微分して

$$\dot{x} = l_1\dot{\theta_1}\cos\theta_1　,　\dot{y} = l_1\dot{\theta_1}\sin\theta_1$$ $$\dot{X} = l_1\dot{\theta_1}\cos\theta_1+l_2\dot{\theta_2}\sin\theta_2　,　\dot{Y} = l_1\dot{\theta_1}\sin\theta_1+l_2\dot{\theta_2}\sin\theta_2$$

おもり１の速度を$v$、おもり２の速度を$V$とすると

$$v^2 = \dot{x}^2+\dot{y}^2 = l_1^2\dot{\theta_1}^2$$
\begin{align} V^2 &= \dot{X}^2+\dot{Y}^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\cos\theta_2+\sin\theta_1\sin\theta_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) \end{align}

よって運動エネルギーは

\begin{align} T &= \frac{1}{2}m_1v^2+\frac{1}{2}m_2V^2\\
&=\frac{1}{2}(m_1+m_2)l_1^2\dot(\theta_1)^2+\frac{1}{2}m_2l_2^2\dot{\theta_2}^2+m_2l_1l_2\dot{\theta_1}\dot{\theta_2}cos(\theta_1-\theta_2) \end{align}

位置エネルギーは

\begin{align} U &= m_1gl_1(1-\cos\theta_1)+m_2g\Bigl(l_1(1-\cos\theta_1)+l_2(1-\cos\theta_2)\Bigr)\\
&= (m_1+m_2)gl_1(1-cos\theta_1)+m_2gl_2(1-cos\theta_2) \end{align}

ラグランジアンは

\begin{align} L &= \frac{1}{2}(m_1+m_2)l_1^2\dot{\theta_1}^2+\frac{1}{2}m_2l_2^2\dot{\theta_2}^2+m_2l_1l_2\dot{\theta_1}\dot{\theta_2}cos(\theta_1-\theta_2)\\
&-(m_1+m_2)gl_1(1-\cos\theta_1)-m_2gl_2(1-\cos\theta_2) \end{align}

簡便のため$m_1=m_2=m$、$l_1=l_2=l$という場合を考えると

\begin{align} L &= ml^2\dot{\theta_1}^2+\frac{m}{2}l^2\dot{\theta_2}+ml^2\dot{\theta_2}\dot{\theta_2}\cos(\theta_1-\theta_2)\\
&-2mgl(1-cos\theta_1)-mgl(2\cos\theta_1+\cos\theta_2-3)\\
&= ml^2\Bigl(\dot{\theta_1}^2+\frac{1}{2}\dot{\theta_2}^2+\dot{\theta_1}\dot{\theta_2}cos(\theta_1-\theta_2)\Bigr)\\
&+mgl(2\cos\theta_1+\cos\theta_2-3) \end{align}

保存力場中でのラグランジュ運動方程式は

$$\frac{d}{dt}(\frac{\partial L }{\partial \dot{\theta_1}})-\frac{\partial L }{\partial \theta_1} = 0$$ $$\frac{d}{dt}(\frac{\partial L }{\partial \dot{\theta_2}})-\frac{\partial L }{\partial \theta_2} = 0$$

それぞれの項を計算すると

$$\frac{d}{dt}\Bigl(\frac{\partial L }{\partial \dot{\theta_1}}\Bigl)= 2ml^2\ddot{\theta_1}+ml^2\Bigl(\ddot{\theta_2}\cos(\theta_1-\theta_2)-\dot{\theta_2}(\dot{\theta_1}-\dot{\theta_2})sin(\theta_1-\theta_2)\Bigl)$$ $$\frac{\partial L }{\partial \theta_1} = -ml^2\dot{\theta_1}\dot{\theta_2}\sin(\theta_1-\theta_2)-2mgl\sin\theta_1$$ $$\frac{d}{dt}\Bigl(\frac{\partial L }{\partial \dot{\theta_2}}\Bigl)= ml^2\ddot{\theta_2}+ml^2\Bigl(\ddot{\theta_1}\cos(\theta_1-\theta_2)-\dot{\theta_1}(\dot{\theta_1}-\dot{\theta_2})sin(\theta_1-\theta_2)\Bigl)$$ $$\frac{\partial L }{\partial \theta_2} = ml^2\dot{\theta_2}\dot{\theta_2}\sin(\theta_1-\theta_2)-mgl\sin\theta_2$$

$\ddot{\theta_1}$、$\ddot{\theta_2}$について整理すると

$$\ddot{\theta_1} = \frac{1}{2}\Bigl(-\ddot{\theta_2}\cos(\theta_1-\theta_2)-\dot{\theta_2}^2\sin(\theta_1-\theta_2)-2\frac{g}{l}\sin\theta_1\Bigl)$$ $$\ddot{\theta_2} = -\ddot{\theta_1}\cos(\theta_1-\theta_2)+\dot{\theta_1}^2\sin(\theta_1-\theta_2)-\frac{g}{l}\sin\theta_2$$

まとめ

$\ddot{\theta_1}$、$\ddot{\theta_2}$に関しての常微分方程式までは導出できたので、次回はルンゲ・クッタ法を用いて数値計算してみたいと思います。