二重振り子のような簡単な系がなぜカオスのような振舞になるのかは非常に興味深い。
参考の①に非常に興味深い説明があるので参照してほしいが、ここではカオスアートを制御することにより、その振舞を見ていこうと思う。

前回、微分方程式に過誤はあったが、カオスアートがある程度制御できることが分かった。
また、二重振り子には少なくとも以下の微分方程式が導出されている。
①単振り子連結
②棒連結
③単振り子連結でも見た目異なる微分方程式が導出されている
これらの系でどういうカオス領域あるいはその結果としてのカオスアートが生成されるのかを調査しようと思う。

やったこと

（１）微分方程式を無次元量で表す
（２）３つの運動方程式から得られるカオスアート
（３）２つの単振り子連結二重振り子の同値性検証
（４）ちょっと驚いたカオスアート

（１）微分方程式を無次元量で表す

そもそも単振り子の運動方程式は

ml\frac{d^2\theta}{dt^2}=-mg\sin\theta

\frac{d^2\theta}{dt^2}=-\frac{g}{l}\sin\theta

上式は触れ角Θ、角速度ω（=dΘ/dt）、そして周期Tで特徴づけられる力学系である。
すなわち微小角における単振り子の周期Tは

T=2\pi\sqrt{\frac{l}{g}}

固有振動数ωは

\omega=\frac{2\pi}{T}
      =\sqrt{\frac{g}{l}}

この固有振動数を用いて、以下のように表現できる。

\theta=C_1\sin\omega t+C_2\cos\omega t

参考のように、この周期と固有振動数は角度が大きい領域でも系を特徴づける量であることが分かる。

【参考】
・振り子

単振り子と同様に、各系の微分方程式を無次元量で表わした。こうすることにより系の振舞を連結する系の固有な振動をイメージしつつ、その絡み合いとしてのカオス領域をイメージしやすくなる。
ということで、以下の関数を利用した。

単振り子連結二重振り子

\theta_1'' = \frac{-\omega^2(\sin\theta_1 - M\sin\theta_2\cos\theta_{12}) -M\sin \theta_{12}(\theta_2'^2L+\theta_1'^2\cos\theta_{12})}{(1-M\cos^2\theta_{12})}

\theta_2'' = \frac{{-\omega^2(\sin\theta_2-\cos\theta_{12}\sin\theta_1)+(\theta_1'^2+ML\theta'^2_2\cos\theta_{12})\sin\theta_{12}}}{L(1-M{\cos^2\theta_{12})}}

def doublePendulum(y, t, L, omegaJ,M):
    theta1, omega1,theta2, omega2 = y
    theta12=theta1 -theta2
    #M=m2/(m1+m2)
    #L=L2/L1
    #omegaJ=g/L1    
    dydt = [omega1,
            (omegaJ*L*(-np.sin(theta1)+M*np.cos(theta12)*np.sin(theta2))-M*L*(omega1*omega1*np.cos(theta12)+L*omega2*omega2)*np.sin(theta12))/(L-M*L*np.cos(theta12)*np.cos(theta12)),
            omega2,
            (omegaJ*(-np.sin(theta2)+np.cos(theta12)*np.sin(theta1))+(M*L*omega2*omega2*np.cos(theta12)+omega1*omega1)*np.sin(theta12))/(L-M*L*np.cos(theta12)*np.cos(theta12))]
    return dydt

棒連結二重振り子

\theta_1'' = \frac{-\omega^2L(M_1\sin\theta_1 - M\sin\theta_2\cos\theta_{12}) -M\sin \theta_{12}(\theta_2'^2L+2\theta_1'^2\cos\theta_{12})}{(1-2M\cos^2\theta_{12})}

\theta_2'' = \frac{{(2\theta_1'^2+2ML\theta'^2_2\cos\theta_{12})\sin\theta_{12}-\omega^2(\sin\theta_2-2M_1\cos\theta_{12}\sin\theta_1)}}{L(1-2M{\cos^2\theta_{12})}}

def doublePendulum2(y, t, L, omegaJ,M,M1):
    theta1, omega1,theta2, omega2 = y
    #M=2*m2/(m1+4*m2)
    #M1=(m1+2*m2)/(m1+4*m2)
    #L=L2/L1
    #omegaJ=g/L1  
    theta12=theta1 -theta2
    dydt = [omega1,
            (omegaJ*L*(-M1*np.sin(theta1)+M*np.cos(theta12)*np.sin(theta2))-M*(2*omega1*omega1*np.cos(theta12)+L*omega2*omega2)*np.sin(theta12))/(1-2*M*np.cos(theta12)*np.cos(theta12)),
            omega2,
            (omegaJ*(-np.sin(theta2)+2*M1*np.cos(theta12)*np.sin(theta1))+(2*M*L*omega2*omega2*np.cos(theta12)+2*omega1*omega1)*np.sin(theta12))/(L-2*M*L*np.cos(theta12)*np.cos(theta12))]
    return dydt

参考③の二重振り子の運動方程式

\theta_1'' = \frac{-\omega^2(\sin\theta_1 + M\sin(\theta_1-2\theta_2)) -2M\sin \theta_{12}(\theta_2'^2L+\theta_1'^2\cos\theta_{12})}{(1-M\cos2\theta_{12})}

\theta_2'' = \frac{2\sin\theta_{12}(M_1\theta_1'^2+\omega^2M_1\cos\theta_1+\theta_2'^2LM_1\cos\theta_{12})}{L(1-M{\cos2\theta_{12})}}

def dualPendulum(y, t, L, omegaJ,M,M1):
    theta1, omega1,theta2, omega2 = y
    #M=m2/(2*m1+m2)
    #M1=(m1+m2)/(2*m1+m2)
    #omegaJ=g/L1
    #L=L2/L1
    dydt = [omega1, 
            (-omegaJ*(np.sin(theta1)+M*np.sin(theta1-2*theta2))-2*np.sin(theta1-theta2)*M*(omega2*omega2*L+omega1*omega1*np.cos(theta1-theta2)))/(1*(1-M*np.cos(2*(theta1-theta2)))), 
            omega2, 
            (2*np.sin(theta1-theta2)*(omega1*omega1*M1+omegaJ*M1*np.cos(theta1)+omega2*omega2*L*M1*np.cos(theta1-theta2)))/(L*(1 - M*np.cos(2*(theta1-theta2))))]
    return dydt

（２）３つの運動方程式から得られるカオスアート

単振り子連結二重振り子

ｘｐ１

ｘｐ２

棒連結二重振り子

ｘｐ１

ｘｐ２

参考③の二重振り子の運動方程式

ｘｐ１

ｘｐ２

（３）２つの単振り子連結二重振り子の同値性検証

そもそも運動方程式の右辺の振舞を関数的に求めてもできるが、ここはカオスアートの同値性から検証してみた。
以下、同じパラメータに対してのカオスアートを比べてみる。

上記の初期角度はだいたい10^-4～10^-5程度の領域で比較しており、ほぼ一致した振舞になっている。
他のカオスアートも角度が小さい領域はほぼ同じ振舞をする。
一方、上記のダイジェスト見てもわかると思うが、いくつか特徴的な絵が出ていない場合がある。以下のように、特に初期角度が大きい場合に差が異なるカオスアートが出現した。ということで、この微分方程式はほぼ単振り子連結二重振り子の運動方程式だといえるが、何か近似等をやっているようだ。
xp1dual_L11L21m1100m21y010.15601991952892066

xp1pen_L11L21m1100m21y010.15601991952892066

xp1dual_L11L21m1100m21y010.19056316004055146

xp1pen_L11L21m1100m21y010.19056316004055146

xp1dual_L11L21m1100m21y010.23275436927724738

xp1pen_L11L21m1100m21y010.23275436927724738

（４）ちょっと驚いたカオスアート

これ何なんだろうというカオスアートがたくさん出てくるが、こういうの出るんですね。
これって、DLA（Diffusion-limited aggregation）で有名な樹枝状成長に類似しています。

そして、以下なんだか不思議なカオスアートをセレクトしてみました。

案外、七角形や五角形が多いし、いわゆる原子のs状態しかも３sとか高次の状態に見えるのがなんとも。。。今回５０００枚以上のカオスアートを見たが、まだまだ別な領域がありそうで、一年がすぐすぎてしまいそうです。

まとめ

・単振り子連結と棒連結二重振り子のカオスアートを見た
・二つの系では、ほぼ似たようなカオスアートが得られる
・低角度領域においても魅力的なカオスアートが得られることが分かった

・強連結な振動子の系は微小振動とか微小角だからと云って、線形近似するのはダメなようだ
・エネルギー保存則の検証をしたいと思う
・運動方程式について、関数のΘやω依存性から上記結論を検証する

おまけ

【参考】
①二重振り子@カオス＆非線形力学入門
②複雑系の物理学(下)琉球大学理学部物質地球科学科集中講義
③Double Pendulum@MyphysicLab.com

棒連結の微分方程式は以下のとおり

\theta_1'' = \frac{-g(m_1 + 2m_2)\sin\theta_1 + 2m_2g\sin\theta_2\cos\theta_{12} -2m_2\sin \theta_{12}(\theta_2'^2L_2+2\theta_1'^2L_1\cos\theta_{12})}{L_1(m_1+4m_2-4m_2\cos^2\theta_{12})}

\theta_2'' = \frac{{(2L_1\theta_1'^2(m_1+4m_2)+4m_2L_2\theta'^2_2\cos\theta_{12})\sin\theta_{12}-g((m_1+4m_2)\sin\theta_2-2(m_1+2m_2)\cos\theta_{12}\sin\theta_1)}}{L_2(m_1+4m_2-4m_2{\cos^2\theta_{12})}}

\theta_1'' = \frac{-\omega^2(M_1\sin\theta_1 - M\sin\theta_2\cos\theta_{12}) -M\sin \theta_{12}(\theta_2'^2L+2\theta_1'^2\cos\theta_{12})}{(1-2M\cos^2\theta_{12})}

\theta_2'' = \frac{{(2\theta_1'^2+2ML\theta'^2_2\cos\theta_{12})\sin\theta_{12}-\omega^2(\sin\theta_2-2M_1\cos\theta_{12}\sin\theta_1)}}{L(1-2M{\cos^2\theta_{12})}}

以下の無次元量とωを導入した。
無次元量で表したことにより、上式の見通しはよくなり、無次元量で記載されている参考の運動方程式と類似したものとなる。

L=\frac{L_2}{L_1}

\omega^2=\frac{g}{L_1}

M_1=\frac{m_1+2m_2}{m_1+4m_2}

M=\frac{2m_2}{m_1+4m_2}

Lagrangianは以下のKとUからL=K-Uで上記微分方程式が得られます。

K=\frac{1}{2}(m_1+4m_2)L_1^2\theta_1'^2+\frac{1}{2}m_2L_2^2\theta_2'^2+2m_2L_1L_2\theta_1'\theta_2'\cos\theta_{12}

U=(m_1+2m_2)gL_1(1-\cos\theta_1)+m_2gL_2(1-\cos\theta_2)