📂동역학

동역학계로써의 이중 진자

모델

평면 상에서 진자^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)