📂동역학

이상하지만 캐어릭하지 않은 어트랙터 SNA

정의 ^{1 2}

이상하지만 캐어릭하지 않은 어트랙터^{SNA, strange nonchaotic attractor}는 프랙털 기하의 구조를 갖추었지만 다이내미컬 센스에서는 캐어릭하지 않은 어트랙터다.

설명

다이내미컬 센스에서 캐어릭하다는 것은 쉽게 말해 랴푸노프 스펙트럼 중에서 가장 큰 수가 음수라는 것으로, 이상한 어트랙터인 것은 맞지만 초기조건에 민감한다는 중요한 성질은 가지지 않는다.

예로써 다음과 같이 자유도가 1인 기어 트랜스미션 시스템^{gear transmission system}을 생각해보자. $$\begin{align*} I_{1} \ddot{\theta}_{1} + c \left( r_{b1} - r_{b2} \dot{\theta}_{2} - \dot{e} (\tau) \right) r_{b1} + k (\tau) f \left( r_{b1} \theta_{1} - r_{b2} \theta_{2} - e (\tau) \right) r_{b1} =& T_{1} \\ I_{2} \ddot{\theta}_{2} + c \left( r_{b1} - r_{b2} \dot{\theta}_{2} - \dot{e} (\tau) \right) r_{b2} + k (\tau) f \left( r_{b1} \theta_{1} - r_{b2} \theta_{2} - e (\tau) \right) r_{b2} =& - T_{2} \end{align*}$$

이 시스템의 역학적인 특징은 지금 관심의 대상이 아니다. 이를 간단한 무차원 방정식^{dimensionless equation}으로 바꾸면 다음과 같이 넌스무스한 시스템으로 요약된다. $$\begin{align*} \dot{x} =& v \\ \dot{v} =& \bar{F}_{n} + \bar{F}_{e} H \left( \cos \theta + \cos \Omega \right) - \left( 1 + \bar{k}_{1} \cos \Omega \right) f(x) - 2 \zeta v \\ \dot{\Omega} =& \omega_{h} \\ \dot{\theta} =& \phi \\ f(x) =& \begin{cases} x - \bar{b} & , \text{if } x > \bar{b} \\ 0 & , \text{if } - \bar{b} \le x \le \bar{b} \\ x + \bar{b} & , \text{if } x < - \bar{b} \end{cases} \end{align*}$$ $\phi$ 는 황금비의 역수 $(\sqrt{5} - 1)/2$ 고, 다른 파라미터는 $\bar{b} = 1$, $\zeta = 0.06$, $\bar{k}_{1} = 0.06$, $H = \omega_{h}^{2} = 1$, $\bar{F}_{n} = 0.3$ 로 두고 바이퍼케이션 파라미터 $\bar{F}_{e} \in [-0.2, 0.2]$ 에 대한 바이퍼케이션 다이어그램과 가장 큰 랴푸노프 지수를 보면 다음과 같다.

이 시스템은 $\bar{F}_{e}$ 의 값에 따라서 다음과 같이 시스템의 특성이 변한다고 한다.

• SNA: $\bar{F}_{e} \in [-0.102, -0.0996] \cup [0.118, 0.122]$
• 캐어릭: $\bar{F}_{e} \in [-0.126, -0.102] \cup [0.122, 0.137]$
• 콰지피리어딕: 위에서 언급된 영역 외에 해당한다.

실제로 $\bar{F}_{e} = -0.1$ 에서의 (a), (b) 푸앙카레 섹션과 (c) 파워 스펙트럼은 위와 같이 나타난다. 푸앙카레 섹션으로 보면 프랙털의 느낌이 나고, 파워스펙트럼 상으로 봤을 땐 이산적인 픽이 나타나지 않기 때문에 피리어딕하거나 콰지피리어딕하지는 않아 이상한 어트랙터는 맞는 것으로 보인다.

그러나 $\bar{F}_{e} = -0.1$ 에서의 랴푸노프 지수를 보면 위와 같이 $0$ 아래에서 수렴하는 것을 확인할 수 있다. 이는 이 때의 오빗이 이상한 어트랙터의 모양을 갖추고 있지만 실제로 캐어릭하지는 않다는 의미가 된다.

코드

다음은 참고논문의 바이퍼케이션 다이어그램과 랴푸노프 스펙트럼을 재현하는 줄리아 코드다. 바이퍼케이션 다이어그램을 그리기 위해서는 푸앙카레 섹션에 따라서 특정 값을 저장해야하는데, 참고논문에 오타가 많아서 고생을 꽤나 했다. 논문에 있는 지침들은 다음과 같이 변경되어야 한다.

• 논문에서는 $\theta = 0 \pmod{2 \pi}$ 에서의 푸앙카레 섹션을 보여준다. 그래서 바이퍼케이션 다이어그램을 그릴 때도 그렇게 하는 것처럼 보이는데, 실제로 바이퍼케이션 다이어그램을 그릴 땐 $\Omega = 0 \pmod{2 \pi}$ 에서의 값을 기록해야 한다.
• 그림에선 $x$ 의 값을 저장하는처럼 보이는데, 명백한 오타고 실제론 푸앙카레 섹션 $\Sigma$ 에서의 $v$ 값을 기록한다.

using DataFrames, CSV, ProgressMeter, LinearAlgebra, Base.Threads

function RK4(f, v::AbstractVector, h=1e-2, nonsmooth=nothing)
    if nonsmooth |> isnothing
        V1 = f(v)
        V2 = f(v + (h/2)*V1)
        V3 = f(v + (h/2)*V2)
        V4 = f(v + h*V3)
    else
        V1 = f(v, nonsmooth)
        V2 = f(v + (h/2)*V1, nonsmooth)
        V3 = f(v + (h/2)*V2, nonsmooth)
        V4 = f(v + h*V3, nonsmooth)
    end
    return v + (h/6)*(V1 + 2V2 + 2V3 + V4), V1
end

function gram_schmidt(J)
    N = size(J, 1)
    U, V = deepcopy(J), deepcopy(J)
    U[:,1] = V[:,1] / norm(V[:,1])
    for j in 2:N
        for jp in 1:(j-1)
            V[:,j] -= (J[:,j]'U[:,jp])*U[:,jp]
        end
        U[:,j] = V[:,j] / norm(V[:,j])
    end
    return U, V
end

function lyapunov_exponent(_data::DataFrame, J_::Function, bf_param;
    U = I(ncol(_data)), T = (last(_data.t) - first(_data.t)))

    λ = zeros(size(U, 1))
    for k = 1:nrow(_data)
        J = J_(collect(_data[k, :])..., bf_param)
        U, V = gram_schmidt(U)
        λ += V |> eachcol .|> norm .|> log
        U = RK4(J, U, dt)
    end
    return sort(λ / T, rev=true)
end

function factory_gear(Fe::Number; ic = [0.1, 0.1, 0.1, 0.0], tspan = [0, 1000], dt = 1e-2)
    __b = 1
    ζ = 0.06
    k1 = 0.06
    ωh = 1
    H = ωh^2
    Fn = 0.3
    φ = (√5 - 1)/2
    function sys(xvΩθ::AbstractVector, nonsmooth::Real)
        x,v,Ω,θ=xvΩθ
        
        ẋ = v
        v̇ = Fn + Fe*H*(cos(θ) + cos(Ω)) - (1 + k1*cos(Ω))*nonsmooth - 2*ζ*v
        Ω̇ = ωh
        θ̇ = φ
        return [ẋ, v̇, Ω̇ , θ̇ ]
    end
    
    t_ = first(tspan):dt:last(tspan)
    len_t_ = length(t_)
    
    t, tk = .0, 0
    u = ic; DIM = length(u)
    traj = zeros(len_t_+2, 2DIM)
    while tk ≤ len_t_
        x,v,Ω,θ = u
        t += dt
        nonsmooth = ifelse(x > __b, x - __b, ifelse(x < -__b, x + __b, 0))
        u, du = RK4(sys, u, dt, nonsmooth)

        if t ≥ first(t_)
            tk += 1
            traj[tk+1,         1:DIM ] =  u
            traj[tk  , DIM .+ (1:DIM)] = du
        end
    end
    return traj[2:(end-2), :]
end
factory_gear(T::Type, args...; kargs...) = 
DataFrame(factory_gear(args...; kargs...), ["x", "v", "Ω", "θ", "dx", "dv", "dΩ", "dθ"])

schedules = DataFrame(idx = 1:801, bp = -0.2:0.0005:0.2)
schedules[!, :λ1] .= .0; schedules[!, :λ2] .= .0; schedules[!, :λ3] .= .0; schedules[!, :λ4] .= .0;
vrbl = [:dx, :dv, :dΩ, :dθ], [:x, :v, :Ω, :θ]
dt = 1e-2; tend = 1000;
function J_(x, v, Ω, θ, Fe)
    dfdx = ifelse(abs(x) > 1, 1, 0)
    return [ 0 1 0 0
             -(1 + k1*cos(Ω))*dfdx -2ζ (-Fe*H*sin(Ω) + k1*sin(Ω)*dfdx) -Fe*H*sin(θ)
             0 0 0 0
             0 0 0 0 ]
end

bfcn = DataFrame(hrzn = [], vrtc = [])
@showprogress @threads for dr = eachrow(schedules)
    data = factory_gear(DataFrame, dr.bp; tspan = [0, tend] .+ 500)
    CSV.write("bifurcation/gear/$(lpad(dr.idx, 5, '0')).csv", data, bom = true)

    λ = lyapunov_exponent(data[:, last(vrbl)], J_, dr.bp, T = tend)
    dr[[:λ1, :λ2, :λ3, :λ4]] .= λ

    data = data[data.Ω .> 5000, :]
    idx_sampled = diff([0; mod.(data.Ω, 2π)]) .< 0
    sampledx = data.v[idx_sampled]
    hrzn, vrtc = fill(dr.bp, length(sampledx)), sampledx
    append!(bfcn, DataFrame(; hrzn, vrtc))
end
CSV.write("lyapunov/gear_bifurcation.csv", bfcn, bom = true)
CSV.write("lyapunov/gear_lyapunov.csv", schedules, bom = true)