tobiom’s diary

備忘録

2重振り子の数値シミュレーション

2重振り子はカオス性を示す物理系として最も有名な例の1つ(だと思っている)。

この2重振り子をシミュレーションしてみた。

解析方針

この2重振り子、微分方程式の数値解法として主流なRunge-Kutta法を適用するとエネルギーが発散してシミュレーションが破綻し易いという事を訊き、FPUのシミュレーションを作る時に利用したSymplectic法を用いる事を考えた。

だが、実装した事のあるSymplectic法はHamiltonianが座標xだけの関数V(x)と運動量pだけの関数T(p)に分離できる前提の陽解法であった。

そこでSymplectic性をもつ陰解法であるGauss-Legendre法を採用。

これによってエネルギーを比較的高精度で保存するシミュレーションを実装する事が可能となった。

しかし、シミュレーションする系がカオスである以上、原理的に現実の系をほぼ寸分のずれも無く再現する事は不可能であるため、あくまでも初期値鋭敏性を分かり易く視覚化するのみのシミュレータになる。

Hamiltonianの導出について

以下、Hamiltonianは運動エネルギー項とポテンシャル項の和となるとする。

2重振り子の系のHamiltonianは形が複雑なため、一旦Lagrangianを導出してからLegendre変換でHamiltonianを求めることにした。

Lagrangian自体は比較的簡単に求められて

L=T-V
T=\displaystyle\frac{1}{2}(m_1+m_2)l_1^2\dot{\theta_1}^2+\displaystyle\frac{1}{2}m_2l_2^2\dot{\theta_2}^2+\displaystyle\frac{1}{2}m_2l_1l_2\dot{\theta_1}\dot{\theta_2}\cos{\left(\theta_1-\theta_2\right)}
V=-(m_1+m_2)gl_1\cos{\theta_1}-m_2gl_2\cos{\theta_2}

である。

ここからLegendre変換p=\displaystyle\frac{\partial L}{\partial \dot{q}}を用いてp_i=\displaystyle\frac{\partial L}{\partial \dot{\theta_i}}を計算する。
このときL=T+Vにおいて\dot{\theta_i}Tの中にのみ現れるのでT偏微分を考えればよろしい。

p_1=(m_1+m_2)l_1^2\dot{\theta_1}+\displaystyle\frac{1}{2}m_2l_1l_2\dot{\theta_2}\cos{\left(\theta_1-\theta_2\right)}

p_2=m_2l_2^2\dot{\theta_2}+\displaystyle\frac{1}{2}m_2l_1l_2\dot{\theta_1}\cos{\left(\theta_1-\theta_2\right)}

この式を用いてL\theta_ip_iの関数に書き換えるとHが得られる。

直接p_i\dot{\theta_i}について解いて代入する方法だと面倒なので、ある程度計算をやり易くするために次の対称行列Aを使う。

A=\begin{pmatrix} (m_1+m_2)l_1^2 & \displaystyle\frac{1}{2}m_2l_1l_2\cos{\left(\theta_1-\theta_2\right)} \\ \displaystyle\frac{1}{2}m_2l_1l_2\cos{\left(\theta_1-\theta_2\right)} & m_2l_2^2 \end{pmatrix}

するとこの行列A\theta_iのベクトル\vec{\bf\theta}=\begin{pmatrix}\theta_1\\\theta_2\end{pmatrix}を用いて\vec{\bf p}=\begin{pmatrix}p_1\\p_2\end{pmatrix}
\vec{\bf p}=A\vec{\bf \dot{\theta}}、即ち\vec{\bf \dot{\theta}}=A^{-1}\vec{\bf p}と表せる。

ここで\theta_iの関数TT=\frac{1}{2}\dot{\vec{\bf\theta}}^tA\vec{\bf\dot{\theta}}であったから、ここに\vec{\bf \dot{\theta}}=A^{-1}\vec{\bf p}を代入すると

T=\frac{1}{2}\left(A^{-1}\vec{\bf p}\right)^tA\left(A^{-1}\vec{\bf p}\right)=\frac{1}{2}\vec{\bf p}^t A^{-1} \vec{\bf p}

こうして得られたp_iの関数としてのTにポテンシャル項Vを足すとHamiltonianが得られる。

運動方程式

以上の結果を正準方程式に適用すると、以下のようになる。

\dot{\vec{\bf p}}=-\displaystyle\frac{\partial (T+V)}{\partial\vec{\bf\theta}}=-\frac{1}{2}\vec{\bf p}^t \displaystyle\frac{\partial A^{-1}}{\partial\vec{\bf\theta}} \vec{\bf p} - \begin{pmatrix}(m_1+m_2)gl_1\sin{\theta_1}\\m_2gl_2\sin{\theta_2}\end{pmatrix}

\dot{\vec{\bf\theta}}=\displaystyle\frac{\partial T}{\partial\vec{\bf p}}=A^{-1} \vec{\bf p}

動画


2重振り子
何故かカクついてる