📂동역학

부알리 어트랙터

모델 ¹

$$ \begin{align*} \dot{x} =& x ( a - y ) + \alpha z \\ \dot{y} =& - y ( b - x^2 ) \\ \dot{z} =& - x ( c - s z ) - \beta x \end{align*} $$

파라미터

• $a = 4$
• $\alpha = 0.3$
• $b = 1$
• $c = 1.5$
• $\beta = 0.05$
• $s = 1$

설명

$$ \begin{align*} \dot{x} =& x ( a - y ) \\ \dot{y} =& - y ( b - x^2 ) \end{align*} $$

부알리 어트랙터^{Bouali attractor}는 위와 같은 롯카-볼테라 모델 타입의 시스템에 $z$ 축을 추가하여 혼돈화^{chaotification} 시킨 모델이다. $\dot{z}$ 의 형태가 꼭 저렇게 되어야하는 것은 아니고, 자연스러운 진동^oscillation을 가지고 피식자인 $x$ 에 피드백을 주기만 하면 된다.

실제로 부알리 방정식이라고 하면 또 다른 유사한 형태로써 다음이 알려져 있고, 전형적인 파라미터는 $\alpha = 3$, $\beta = 2.2$, $\gamma = 1$, $\mu = 0.001$ 이다².

$$ \begin{align*} \dot{x} =& a x (1 - y) - \beta z \\ \dot{y} =& - \gamma y (x - 1) \\ \dot{z} =& \mu x \end{align*} $$

어트랙터

$s = 1$ 일 때의 이상한 어트랙터다. 로렌츠 어트랙터와 비교하면 날개 아래로 뻗어나가는 무언가가 생겨난 느낌을 받을 수 있다.

바이퍼케이션

$s$ 를 0.5부터 1.2까지 변화시키면서 그린 바이퍼케이션 다이어그램이다. 최대차항이 3차라서 그런지, 아이자와 어트랙터 못지 않게 기괴한 형태를 보여준다.

코드

다음은 부알리 어트랙터를 시뮬레이션 할 수 있는 줄리아 코드다.

using DataFrames, DifferentialEquations, OrdinaryDiffEqLowOrderRK

function factory_bouali(s::Number; ic = [1.0, 1.0, 0.1], saveat = 0:1e-2:100)
      a,   α, b,   c,    β = (
      4, 0.3, 1, 1.5, 0.05)
    function sys(du, u, p, t)
        x, y, z = u
        s = p[1]

        du[1] = x*(a - y) + α*z
        du[2] = -y*(b - x^2)
        du[3] = -x*(c - s*z) - β*z
        return du
    end
    sol = solve(ODEProblem(sys, ic, (0, last(saveat)), [s]), RK4(), dt = saveat.step.hi, adaptive=false, maxiters = Inf)
    dsol = [sol.k[2][1] stack([solk[2] for solk in sol.k[2:end]])]
    matrix = [sol.t sol[:, :]' dsol']
    return matrix[sol.t .≥ first(saveat), :][1:end-1, :]
end
factory_bouali(T::Type, args...; kargs...) =
DataFrame(factory_bouali(args...; kargs...), ["t", "x", "y", "z", "dx", "dy", "dz"])

sol = factory_bouali(DataFrame, 1, ic = [2,1,0], saveat = 4500:1e-2:5000)
plot(sol.x, sol.y, sol.z, alpha = .5, color = :black, xlabel = L"x", ylabel = L"y", zlabel = L"z", size = [400, 400])

pm, pM = 0.5, 1.2; p0, p1 = 0.7, 0.8;
p_ = range(pm, pM, length = 1001)
bfcn = callbfcn("G:/BF/bouali/bfcnA.jld2")
@showprogress @threads for k in eachindex(p_)
    sol = factory_bouali(DataFrame, p_[k], saveat = 2500:1e-2:3000)
    bfcn[p_[k]] = sol.x[arglmax(sol.x)]
end
scatter(dict2bifurcation(bfcn)..., xticks = [pm, p0, p1, pM], ms = .5, ma = .5, msw = 0, color = :black); png("temp1")