📂力学系

トーマス・アトラクターとラビリンス・カオス

モデル¹

$$ \begin{align*} \dot{x} =& - b x + \sin y \\ \dot{y} =& - b y + \sin z \\ \dot{z} =& - b z + \sin x \end{align*} $$

• $x, y, z$: 3次元空間における座標
• $b$: 減衰係数^{damping coefficient}

説明²

トーマス・アトラクターはいわゆる循環対称性^{cyclically symmetricity}を持ち、変数の順序がどのように入れ替わっても同じ形になる支配方程式を持っている。

このシステムで見られるカオスは迷宮カオス^{labyrinth chaos}と呼ばれる。上の図は$b = 0.1$のときのトラジェクトリを表したものである。

多重安定性

例えば$b = 0.16$のときを見ると、同じシステムであっても初期条件に応じてリミットサイクルが複数存在し、多重安定性を持つことがわかる。

分岐

参考論文では上のように分岐図を描いているが、実際には多重安定性があるため、さまざまな初期条件についてそれぞれの分岐図を重ね合わせて描く必要がある。

例えば4つの初期条件については、上のようにさまざまな分岐図が描かれる。

ランダムな初期条件については上のようになる。

コード

以下は本文の画像を再現するJuliaコードである。

using JLD2, ProgressMeter, DataFrames, DifferentialEquations, Plots, StatsBase

function factory_thomas(b::Number; ic = rand(3), saveat = 0:1e-2:2000)
    function sys(du, u, p, t)
        x, y, z = u
        b = p[1]

        du[1] = sin(y) - b*x
        du[2] = sin(z) - b*y
        du[3] = sin(x) - b*z
        return du
    end
    sol = solve(ODEProblem(sys, ic, (0, last(saveat)), [b]), RK4(), dt = saveat.step.hi, adaptive=false, maxiters = 1e+7)
    matrix = Matrix([sol.t'; sol[:, :]; stack([sys(zeros(3), u, [b], 0) for u in sol.u])]')
    return matrix[sol.t .≥ first(saveat), :][1:end-1, :]
end
factory_thomas(T::Type, args...; kargs...) =
DataFrame(factory_thomas(args...; kargs...), ["t", "x", "y", "z", "dx", "dy", "dz"])

data = factory_thomas(DataFrame, 0.10, saveat = 1000:1e-2:2000)
plt_0 = plot(data.x, data.y, data.z, color = :black, size = [400, 400])

data = factory_thomas(DataFrame, 0.16, ic = [.1, .2, .3], saveat = 1000:1e-2:2000)
plt_1 = plot(data.x, data.y, data.z, color = :black, size = [400, 400], title = "ic = (.1, .2, .3)")
data = factory_thomas(DataFrame, 0.16, ic = [.1, .2, .5], saveat = 1000:1e-2:2000)
plt_2 = plot(data.x, data.y, data.z, color = :black, size = [400, 400], title = "ic = (.1, .2, .5)")
data = factory_thomas(DataFrame, 0.16, ic = [.1, .2, .8], saveat = 1000:1e-2:2000)
plt_3 = plot(data.x, data.y, data.z, color = :black, size = [400, 400], title = "ic = (.1, .2, .8)")
plot(plt_1, plt_2, plt_3, size = [600, 200], layout = (1, 3))

b_ = .10:2e-5:.24
if !isfile("bifurcation_thomas.jld2")
    @info "Calculating bifurcation data for Thomas..."
    bfcn = Dict{Float64, Vector{Float64}}()
    @showprogress @threads for k in eachindex(b_)
        sol = factory_thomas(DataFrame, b_[k], ic = ic0)
        x_ = sol.x[sol.t .≥ 1000]
        bfcn[b_[k]] = x_[arglmax(x_)]
    end
    JLD2.@save "bifurcation_thomas.jld2" bfcn
else
    @info "Loading bifurcation data for Thomas from file..."
    JLD2.@load "bifurcation_thomas.jld2" bfcn
end
scatter(dict2bifurcation(bfcn)..., ms = .5, ma = .5, msw = 0, color = :black); png("temp")