📂Dynamics

Strange Nonchaotic Attractor (SNA)

Definition ^{1 2}

A Strange Nonchaotic Attractor is an attractor that possesses the structure of fractal geometry but is not chaotic in a dynamical sense.

Explanation

In a dynamical sense, something being chaotic means that the largest number in the Lyapunov spectrum is negative, indicating that, although it is a strange attractor, it does not have the important property of being sensitive to initial conditions.

Consider, for instance, a gear transmission system with a degree of freedom of 1.
$$\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*}$$

The dynamical characteristics of this system are not the focus here. When simplified into a dimensionless equation, it can be represented as a nonsmooth system.
$$\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$ is the reciprocal of the golden ratio $(\sqrt{5} - 1)/2$, and other parameters are set as $\bar{b} = 1$, $\zeta = 0.06$, $\bar{k}_{1} = 0.06$, $H = \omega_{h}^{2} = 1$, $\bar{F}_{n} = 0.3$. Viewing the bifurcation diagram and the largest Lyapunov exponent concerning the bifurcation parameter $\bar{F}_{e} \in [-0.2, 0.2]$ gives the following:

It is said that the system’s characteristics change according to the value of $\bar{F}_{e}$ as follows:

• SNA: $\bar{F}_{e} \in [-0.102, -0.0996] \cup [0.118, 0.122]$
• Chaotic: $\bar{F}_{e} \in [-0.126, -0.102] \cup [0.122, 0.137]$
• Quasiperiodic: Outside the mentioned regions.

In reality, the (a), (b) Poincaré sections, and (c) power spectrum at $\bar{F}_{e} = -0.1$ appear as shown above. The Poincaré sections give a fractal sense, and since discrete peaks do not appear in the power spectrum, it is neither periodic nor quasiperiodic, hence confirming it as a strange attractor.

However, by examining the Lyapunov exponent at $\bar{F}_{e} = -0.1$, it can be observed that it converges below $0$ as seen above. This means that although the orbit has the shape of a strange attractor, it is not actually chaotic.

Code

The following is a Julia code that replicates the bifurcation diagram and Lyapunov spectrum from the reference paper. To draw the bifurcation diagram, specific values need to be stored according to the Poincaré sections, and there were quite a few errors in the reference paper, causing some difficulties. The guidelines in the paper should be changed as follows:

• The paper shows a Poincaré section at $\theta = 0 \pmod{2 \pi}$, and it appears that the bifurcation diagram should also be drawn similarly, but in fact, values at $\Omega = 0 \pmod{2 \pi}$ should be recorded for the bifurcation diagram.
• The figure suggests saving values at $x$, but it’s an obvious typo; instead, values at $v$ should be recorded in the Poincaré section at $\Sigma$.

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)