📂Dynamical Systems

Thomas Attractor and Labyrinth Chaos

Model¹

$$ \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$: coordinates in three-dimensional space
• $b$: damping coefficient

Explanation²

The Thomas attractor possesses a so-called cyclically symmetricity, so it has governing equations that take the same form no matter how the order of the variables is permuted.

The chaos seen in this system is called labyrinth chaos. The figure above shows the trajectory when $b = 0.1$.

Multistability

For example, looking at the case $b = 0.16$, we can see that even for the same system there exist several limit cycles depending on the initial conditions, so it has multistability.

Bifurcation

The reference paper drew the bifurcation diagram as above, but in reality there is multistability, so one must overlay the individual bifurcation diagrams for various initial conditions.

For instance, for four different initial conditions, several bifurcation diagrams are drawn as above.

For random initial conditions, the result is as above.

Code

The following is the Julia code that reproduces the images in this article.

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")