📂動力学

力学系としての二重振り子

モデル

平面上で振り子^pendulumの先にもう一つの振り子がある構造を二重振り子^{double pendulum}と呼ぶ¹。

$$ \begin{align*} \ddot{\phi} =& \left[ 1 - \mu \cos^{2} \left( \psi - \phi \right) \right]^{-1} \left[ \mu g_{1} \sin \psi \cos \left( \psi - \phi \right) \right. \\ & \left. + \mu \dot{\phi}^{2} \sin \left( \psi - \phi \right) \cos \left( \psi - \phi \right) - g_{1} \sin \phi + {\frac{ \mu }{ \lambda }} \dot{\psi} \sin \left( \psi - \phi \right) \right] \\ \ddot{\psi} =& \left[ 1 - \mu \cos^{2} \left( \psi - \phi \right) \right]^{-1} \left[ g_{2} \sin \phi \cos \left( \psi - \phi \right) \right. \\ & \left. - \mu \dot{\psi}^{2} \sin \left( \psi - \phi \right) \cos \left( \psi - \phi \right) - g_{2} \sin \psi - \lambda \dot{\phi} \sin \left( \psi - \phi \right) \right] \end{align*} $$

変数

• $\phi (t)$: $t$ 時点で原点と第一振り子の間の角度だ。
• $\psi (t)$: $t$ 時点で第一振り子と第二振り子の間の角度だ。

パラメータ

• $g \approx 9.80665$: 重力加速度だ。
• $l_{1}$: 原点と第一振り子の間の距離だ。
• $l_{2}$: 第一振り子と第二振り子の間の距離だ。
• $m_{1}$: 第一振り子の質量だ。
• $m_{2}$: 第二振り子の質量だ。

$\lambda$、$g_{1}$、$g_{2}$、$\mu$ は次のように定義される。 $$ \begin{align*} \lambda =& l_{1} / l_{2} \\ g_{1} =& g / l_{1} \\ g_{2} =& g / l_{2} \\ \mu =& m_{2} / M \end{align*} $$

第一振り子の位置を $\left( x_{1}, y_{1} \right)$ とし、第二振り子の位置を $\left( x_{1}, y_{1} \right)$ とすると次を満たす。 $$ \begin{align*} x_{1} =& + l_{1} \sin \phi \\ y_{1} =& - l_{1} \cos \phi \\ x_{2} =& + l_{1} \sin \phi + l_{2} \sin \psi \\ y_{2} =& - l_{1} \cos \phi - l_{2} \cos \psi \end{align*} $$ この方程式を通じて角度の大きさから振り子の位置を正確に計算する。

説明

混沌

ダイナミクスの最も基本的なシステム、特にカオスを示す例としては、写像で定義されるシステム基準でロジスティック写像、微分方程式で定義される動力学系基準でローレンツアトラクタがある。しかし、これらはその直感的なモチーフを理解するまでにかなりの時間がかかるが、二重振り子は特別な説明なしにそれがなぜカオティック^chaoticであるかが分かる。

コード

軌道

以下はRK4を通じて二重振り子の方程式を解き視覚化するジュリアコードだ。

using Plots, DataFrames

h = 1e-3
function RK4(f, y)
    V1 = f(y)
    V2 = f(y + h*V1/2)
    V3 = f(y + h*V2/2)
    V4 = f(y + h*V3)
    return y + h/6*(V1 + 2V2 + 2V3 + V4)
end

g = 9.80665
l1 = 2
l2 = 1
m1 = 1
m2 = 1
λ = l1 / l2
μ = m2 / (m1 + m2)
g1 = g / l1
g2 = g / l2
function f(pppp)
    ϕ, ψ, Φ, Ψ = pppp
    ϕ̇ = Φ
    ψ̇ = Ψ
    Φ̇ = (μ*g1*sin(ψ)*cos(ψ-ϕ) + μ*(Φ^2)*sin(ψ-ϕ)*cos(ψ-ϕ) - g1*sin(ϕ) + (μ/λ)*(Ψ^2)*sin(ψ-ϕ)) / (1-μ*(cos(ψ-ϕ)^2))
    Ψ̇ = (g2*sin(ϕ)*cos(ψ-ϕ) - μ*(Ψ^2)*sin(ψ-ϕ)*cos(ψ-ϕ) - g2*sin(ψ) - λ*(Φ^2)*sin(ψ-ϕ)) / (1-μ*(cos(ψ-ϕ)^2))
    return [ϕ̇ , ψ̇ , Φ̇ , Ψ̇ ]
end

x_ = [[1.0,1,2,2]]
for t in 0:h:50
    push!(x_, RK4(f, x_[end]))
end
_x_ = stack(x_, dims = 1)
ϕ = _x_[:, 1]
ψ = _x_[:, 2]

x1 = +l1*sin.(ϕ)
y1 = -l1*cos.(ϕ)
x2 = +l1*sin.(ϕ) + l2*sin.(ψ)
y2 = -l1*cos.(ϕ) - l2*cos.(ψ)

df = DataFrame(; ϕ, ψ, x1, x2, y1, y2)

p0 = plot(size = [600, 600], xlims = [-3, 3], ylims = [-4, 2], legend = :none)
@time a1 = @animate for k = 1:20:nrow(df)
    p1 = plot(p0, df.x2[1:10:k], df.y2[1:10:k], color = :black, alpha = (1:10:k)./2k)
    p2 = plot(p1, [0, df.x1[k]], [0, df.y1[k]], lw = 2)
    p3 = plot(p2, [df.x1[k], df.x2[k]], [df.y1[k], df.y2[k]], lw = 2)
end
mp4(a1, "double_pendulum.mp4", fps = 60)