📂동역학

아이자와 어트랙터

모델

$$ \begin{align*} \dot{x} =& (z - b)x - d y \\ \dot{y} =& d x + (z - b)y \\ \dot{z} =& c + a z - \frac{z^{3}}{3} - (x^{2} + y^{2})(1 + e z) + f z x^{3} \end{align*} $$

파라미터

• $a = 0.95$
• $b = 0.7$
• $c = 3.5$
• $d = 0.25$
• $e = 0.1$
• $f = 0.1$

설명

아이자와 어트랙터^{Aizawa attractor}는 우변에 최대 4차항까지 포함해서 복잡한 혼돈을 보여주는 동역학계로써 시계열 예측 기법 등에서 벤치마크로써 등장하기도 한다.

그러나 필자는 이 시스템의 기원과 유래를 살펴보던 중 아이자와 본인의 논문이 존재하지 않으며, 그나마 아이자와 어트랙터가 등장하는 문헌들에서도 방정식의 표현이 제각각인 것을 알게 되었다^{1 2 3 4 5}. 이미 명성을 얻어서 인용되고 있는 어트랙터치고는 꽤나 기묘한 일이다. 본문에는 상미분방정식은 그들 중 가장 보편적인 형태라고 생각되는 표현을 모아서 정리하고 남겼다.

어트랙터

이상한 어트랙터 중에서도 상당히 이상하게 생겼다. 기본적으로 구의 형태를 취하면서도 $z$-축을 따라 회전하며 오르내리는 구조를 가졌다.

바이퍼케이션

적어도 이 포스트를 작성하는 지금까지도, 필자는 아이자와 어트랙터와 관련된 바이퍼케이션을 언급한 연구 결과를 본 적이 없다. 위의 바이퍼케이션 다이어그램은 $b$ 를 $0.6$ 에서 $0.9$ 까지 변화시키면서, $z$ 의 최대값을 기록한 것이다. 보다시피 흔히 보여야할 루트 투 카오스^{route to chaos} 등이 보이지 않고 굉장히 기괴한 형태를 가지고 있다. 내가 지금은 아이자와 어트랙터를 연구할 시간이 없지만, 언젠가는 이 시스템에 대해 깊게 파고들어보고 싶다.

코드

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

using DataFrames, DifferentialEquations, OrdinaryDiffEqLowOrderRK

function factory_aizawa(b::Number; ic = [0.1, 0.1, 0.1], saveat = 0:1e-2:100)
      a,  c,   d,   e,  f = (
    .95, .6, 3.5, .25, .1)
    function sys(du, u, p, t)
        x, y, z = u
        b = p[1]

        du[1] = (z - b)*x - d*y
        du[2] = d*x + (z - b)*y
        du[3] = c + a*z - (z^3)/3 - (x^2 + y^2)*(1 + e*z) + f*z*x^3
        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_aizawa(T::Type, args...; kargs...) =
DataFrame(factory_aizawa(args...; kargs...), ["t", "x", "y", "z", "dx", "dy", "dz"])

sol = factory_aizawa(DataFrame, 0.7)
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.6, 0.9; p0, p1 = 0.72, 0.74;
p_ = range(pm, pM, length = 1001)
bfcn = callbfcn("G:/BF/aizawa/bfcnA.jld2")
@showprogress @threads for k in eachindex(p_)
    sol = factory_aizawa(DataFrame, p_[k], saveat = 0:1e-3:3000)
    z_ = sol.z[sol.t .≥ 2000]
    bfcn[p_[k]] = z_[arglmax(z_)]
end