Juliaで二重振り子シミュレーションを作ってみた

最近、懐かしい力学系を動かしたくなってJulia入門しました。

書き方はPythonやRに近く、思ったより早く動くモノが完成しました。勉強に使った書籍は↓です。

Juliaではじめる数値計算入門

題材としては、力学系のカオス入門である二重振り子を選びました。最終的にグワングワン動く絵が描けて楽しかったです！なので、今回はこちらを紹介していきます。

例のごとく、レポジトリを公開しています。良かったらCloneして遊んでください！

二重振り子

二重振り子は力学系のカオス解を得られる典型例としてよく使われます。ここでのカオスとは、ほんの少しの初期条件の違いでも全く異なる複雑な動きが得られることを指します。

二重振り子は二つの振り子を縦に連結させたもので、実世界でもシミュレーションすることができます。でんじろう先生の素晴らしい動画をみてください！

僅かな角度の違いで、全く違う動きがでてきて面白いですね。一見これをコンピュータで計算するのはとても難しいように思います。

しかし、微分方程式を利用することで比較的簡単にこの二重振り子の解を求めることができます。

二重振り子の運動方程式

詳しい導出は他の記事に丸投げしますが、二重振り子の運動方程式は以下のようになります。パラメータがいっぱいなので慣れていない方は流してもらって大丈夫ですw。

$$
\begin{aligned}
&\theta_{1}: 内側の振り子の角度 \\
&\theta_{2}: 外側の振り子の角度 \\
&\omega_{1}: 内側の振り子の角速度 \\
&\omega_{2}: 外側の振り子の角速度 \\
&l_{1}: 内側の振り子の長さ \\
&l_{2}: 外側の振り子の長さ \\
&w_{1}: 内側の振り子の質量 \\
&w_{2}: 外側の振り子の質量 \\
&g: 重力加速度 \coloneqq 9.8 \\
\\
&C = cos(\theta_{1} - \theta_{2}) \\
&S = sin(\theta_{1} - \theta_{2}) \\
&m = w_{1} + w_{2} \\
\\
&\dot{\theta_{1}} = \omega_{1} \\
&\dot{\theta_{2}} = \omega_{2} \\
&\dot{\omega_{1}} = \frac{g(sin(\theta_{2}) C - \frac{m}{w_{2}} sin(\theta_{1})) - (l_{1} \omega_{1}^{2} C + l_{2} \omega_{2}^{2})S}{l_{1} (\frac{m}{w_{2}} - C^{2})} \\
&\dot{\omega_{2}} = \frac{g \frac{m}{w_{2}}(sin(\theta_{1}) C - sin(\theta_{2})) + (\frac{m}{w_{2}} l_{1} \omega_{1}^{2} + l_{2} \omega_{2}^{2} C)S}{l_{2} (\frac{m}{w_{2}} - C^{2})} \\
\end{aligned}
$$

角度と角速度以外のパラメータは最初に指定する定数です。そして、角度と角速度は連立一次微分方程式で与えられます。この方程式は解くことができません（一般解を得ることができない）が、角度と角速度の初期値（=時間0のときの値）を与えることで、逐次的にこれらの数値を計算することができます。

※厳密には逐次計算で得られるのは近似解です。

Juliaによる実装

初期条件

まず、シミュレーションで使うパラメータを入力します。上で書いた運動方程式で使用するパラメータの他に、時間長(t_max)と時間刻み(t_delta)を指定します。皆さんが使うときには、このパラメータを色々変えてみて下さい。

特に、シミュレーションにおけるパラメータではこの時間刻みが重要です。大きすぎると計算は早く終わりますが動きがカクカクしてしまい、詳細な動きを追うことができませんが、細かすぎると計算が長くなり解が発散してしまう可能性があります。コンピューターサイエンス理論ではこのシミュレーションをなるべく細かく行うための様々な工夫がなされています。

Base.@kwdef struct DoublePendulumParameters
    pendulum_length::NTuple{2, Float16} = (1.0, 0.5)
    pendulum_weight::NTuple{2, Float16} = (0.5, 0.2)
    gravity::Float16 = 9.8

    theta_init::Vector{Float64} = [3 * π / 4, 3 * π / 4]
    omega_init::Vector{Float64} = [0.0, 0.0]

    t_max::Float16 = 60.0
    t_delta::Float16 = 0.05

    locus_max::Float16 = 2.0
end

運動方程式

上で書いた運動方程式をプログラムに書きました。（パラメータが多すぎて、2回くらい間違えましたw。）

function motion_equation(
    theta::Vector{Float64},
    omega::Vector{Float64},
    params::DoublePendulumParameters,
)
    len = params.pendulum_length
    weight = params.pendulum_weight
    gravity = params.gravity

    c_const = cos(theta[1] - theta[2])
    s_const = sin(theta[1] - theta[2])
    m_const = (weight[1] + weight[2]) / weight[2]

    values::Vector{Float64} = zeros(Float64, 2)
    values[1] = (
        gravity * (c_const * sin(theta[2]) - m_const * sin(theta[1]))
        - (len[1] * omega[1] ^ 2 * c_const + len[2] * omega[2] ^ 2) * s_const
    ) / (len[1] * (m_const - c_const ^ 2))
    values[2] = (
        gravity * m_const * (c_const * sin(theta[1]) - sin(theta[2]))
        + (m_const * len[1] * omega[1] ^ 2 + len[2] * omega[2] ^ 2 * c_const) * s_const
    ) / (len[2] * (m_const - c_const ^ 2))

    return values
end

数値計算

常微分方程式を初期値問題で解くときに、上の関数を逐次更新していっても近似解は得られません。上の方程式式で求まるのはあくまで微分値なので、一工夫加える必要があります。そのための方法をしてルンゲクッタ法があります。（シミュレーションをするくらいであれば、これを使うくらいの理解で大丈夫です。）

function runge_kutta(
    func,
    theta::Vector{Float64},
    omega::Vector{Float64},
    params::DoublePendulumParameters,
)
    k1_theta = omega * params.t_delta
    k1_omega = func(theta, omega, params) * params.t_delta

    k2_theta = (omega + k1_omega / 2) * params.t_delta
    k2_omega = func(theta + k1_theta / 2, omega + k1_omega / 2, params) * params.t_delta

    k3_theta = (omega + k2_omega / 2) * params.t_delta
    k3_omega = func(theta + k2_theta / 2, omega + k2_omega / 2, params) * params.t_delta

    k4_theta = (omega + k3_omega) * params.t_delta
    k4_omega = func(theta + k3_theta, omega + k3_omega / 2, params) * params.t_delta

    new_theta = theta + (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta) / 6
    new_omega = omega + (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega) / 6

    return new_theta, new_omega
end

これまでに定義したパラメータや関数を使って、パラメータの解を計算していきます。加えて、描画用に得られた角度から振り子の先端の空間座標も計算しておきます。

function solution(params::DoublePendulumParameters)
    times = [params.t_delta:params.t_delta:params.t_max;]
    points = PendulumPoints([], [], [])

    theta = params.theta_init
    omega = params.omega_init
    pendulum_points = trans_point(theta, params)

    push!(points.time, 0.0)
    push!(points.pendulum_1_point, pendulum_points[1])
    push!(points.pendulum_2_point, pendulum_points[2])

    for now_t in times
        theta, omega = runge_kutta(motion_equation, theta, omega, params)
        pendulum_points = trans_point(theta, params)

        push!(points.time, now_t)
        push!(points.pendulum_1_point, pendulum_points[1])
        push!(points.pendulum_2_point, pendulum_points[2])
    end

    return points
end

function trans_point(theta::Vector{Float64}, params::DoublePendulumParameters)
    pendulum_1::Vector{Float64} = [0, 0]
    pendulum_2::Vector{Float64} = [0, 0]

    pendulum_1[1] = params.pendulum_length[1] * sin(theta[1])
    pendulum_1[2] = - params.pendulum_length[1] * cos(theta[1])
    pendulum_2[1] = pendulum_1[1] + params.pendulum_length[2] * sin(theta[2])
    pendulum_2[2] = pendulum_1[2] - params.pendulum_length[2] * cos(theta[2])

    return pendulum_1, pendulum_2
end

グラフ化

得られた座標をいい感じに描画してGIF形式で保存します。

function create_pendulum_gif(points::PendulumPoints, params::DoublePendulumParameters)
    times = points.time
    pendulum_1 = points.pendulum_1_point
    pendulum_2 = points.pendulum_2_point

    len = (params.pendulum_length[1] + params.pendulum_length[2]) * 1.1

    locus_x::Vector{Float64} = []
    locus_y::Vector{Float64} = []
    locus_t::Vector{Float64} = []
    all_locus_x::Vector{Float64} = []
    all_locus_y::Vector{Float64} = []
    all_locus_t::Vector{Float64} = []

    anim = @animate for idx in 1:length(times)
        push!(locus_x, pendulum_2[idx][1])
        push!(locus_y, pendulum_2[idx][2])
        push!(locus_t, times[idx])
        push!(all_locus_x, pendulum_2[idx][1])
        push!(all_locus_y, pendulum_2[idx][2])
        push!(all_locus_t, times[idx])

        if length(locus_x) > params.locus_max / params.t_delta
           popfirst!(locus_x) 
           popfirst!(locus_y) 
           popfirst!(locus_t) 
        end

        dinamics = plot(
            [0, pendulum_1[idx][1], pendulum_2[idx][1]],
            [0, pendulum_1[idx][2], pendulum_2[idx][2]],
            title = "Pendulum dinamics",
            lw = 3,
            lc = "black",
            markershape = :circle,
            ms = 5,
            mc = "black",
        )
        if length(locus_x) >= 2
            plot!(
                dinamics,
                locus_x,
                locus_y,
                linez=locus_t,
                lw = 5,
                c = :binary,
                legend = false,
                colorbar=false,
            )
        end

        locus = plot(
            title = "Pendulum locus",
            all_locus_x,
            all_locus_y,
            linez=all_locus_t,
            lw = 5,
            c = :jet,
            la = 0.5,
            legend = false,
        )

        plot(
            dinamics,
            locus,
            layout = (1, 2),
            plot_title = "time: $(times[idx])",
            legend = false,
            xaxis = false,
            yaxis = false,
            aspect_ratio = :equal,
            size = (1000, 500),
        )
        xlims!(-len, len)
        ylims!(-len, len)
    end

    mkpath("output/double_pendulum")
    gif(anim, "output/double_pendulum/simulation.gif", fps = Int.(1 / params.t_delta))
end

実際を描画を貼り付けてみました。左側が振り子の動きで、右側が外側の振り子の軌跡です。

めちゃくちゃ複雑に動いていますね！完璧だと思います。

二重振り子シミュレーション（note用に小さくしています）

ということで、Julia入門した+力学系したい欲を満たしたの記事でした。

Juliaは思ったより良い言語だと思いました。特に、for文を気兼ねなく書ける点が素晴らしいと思います。これからもっと流行るかもですね！

また、微分方程式は本当に色々なことができて楽しいので、ものづくりが好きな人ははまると思います！是非入門してみてください！（しかも、工作と比べてお金がかからない！w）

最後までお読みいただきありがとうございました。