📂동역학

토마스 어트랙터와 미궁 혼돈

모델 ¹

$$ \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$ 일 때를 보면 같은 시스템이라도 초기 조건에 따라 리미트 사이클이 여러가지 존재해서, 다중안정성을 가진다는 것을 알 수 있다.

바이퍼케이션

참고논문에서는 위와 같이 바이퍼케이션 다이어그램을 그렸지만, 실제로는 다중안정성이 있어서 여러가지 초기조건에 대해서 각각의 바이퍼케이션 다이어그램을 중첩해서 그려야 한다.

가령 네 가지의 초기조건에 대해서는 위와 같이 여러 가지의 바이퍼케이션 다이어그램이 그려진다.

랜덤한 초기조건에 대해서는 위와 같다.

코드

다음은 본문의 이미지들을 재현하는 줄리아 코드다.

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