📂動力学

ストレンジノンカオティックアトラクター (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)