📂Dynamical Systems

Lorenz 96 System

Model ¹

Lorenz-96 system consists of $N$ variables $x_{1}, x_{2}, \ldots, x_{N}$ and is described by the following ordinary differential equations. $$ {\frac{ d x_{k} }{ dt }} = x_{k-1} \left( x_{k+1} - x_{k-2} \right) - x_{k} + F(t) \qquad k = 1 , \cdots , N \pmod{N} $$

Here, $F(t)$ is the term representing the external forcing.

Description

By contrast, the classical Lorenz attractor is Lorenz-63². Unlike Lorenz-63, Lorenz-96 can be defined without restriction on the natural number $N$, and therefore is widely used as a high-dimensional and chaotic system for benchmarking prediction methods.

Bifurcation

The following is a bifurcation diagram showing how the system changes as the value of $A$ varies, in the case where $N = 6$ and $F(t) = A \sin(2t) + 2$ hold.

Code

function factory_lorenz96(A::Number; ic = [0,0,1,1,1,1,1], tspan = [0., 10.], dt = 1e-4)
    function sys(v::AbstractVector)
        T, x1, x2, x3, x4, x5, x6 = v

        dT = 1
        dx1 = x6*(x2 - x5) - x1
        dx2 = x1*(x3 - x6) - x2
        dx3 = x2*(x4 - x1) - x3
        dx4 = x3*(x5 - x2) - x4
        dx5 = x4*(x6 - x3) - x5
        dx6 = x5*(x1 - x4) - x6
        return [dT; [dx1, dx2, dx3, dx4, dx5, dx6] .+ (A*sin(2T) + 2)]
    end

        
    t_ = first(tspan):dt:last(tspan)
    len_t_ = length(t_)
    
    t, tk = .0, 0
    v = ic; DIM = length(v)
    traj = zeros(len_t_+2, 2DIM)
    while tk ≤ len_t_
        t,_,_ = v
        v, dv = RK4(sys, v, dt)

        if t ≥ first(t_)
            tk += 1
            traj[tk+1,         1:DIM ] =  v
            traj[tk  , DIM .+ (1:DIM)] = dv
        end
    end
    return traj[2:(end-2), :]
end
factory_lorenz96(T::Type, args...; kargs...) =
DataFrame(factory_lorenz96(args...; kargs...), ["t", "x1", "x2", "x3", "x4", "x5", "x6", "dt", "dx1", "dx2", "dx3", "dx4", "dx5", "dx6"])

A_ = 1:0.01:4
@showprogress @threads for k in eachindex(A_)
    A = A_[k]
    traj = factory_lorenz96(DataFrame, A, tspan = [0., 2500], dt = 1e-3)
    _traj = traj[150001:end, :]
    CSV.write("G:/lorenz96/lorenz96_$(lpad(k, 4, '0')).csv", _traj)
end

function arglmax(x)
    bits = circshift(x, 1) .< x .> circshift(x, -1)
    bits[1] = false
    bits[end] = false
    return findall(bits)
end
scatter(repeat([3.2], length(_traj.x1[almx1])), _traj.x1[almx1])

plt_bfcn = plot(legend = :none)
@showprogress for k in eachindex(A_)
    _traj = CSV.read("G:/lorenz96SINDy/lorenz96SINDy_$(lpad(k, 4, '0')).csv", DataFrame)
    x1 = _traj.x1[_traj.t .> 1000]
    almx1 = arglmax(x1)
    scatter!(plt_bfcn, repeat([α_[k]], length(x1[almx1])), x1[almx1], color = :red, ms = 1, ma = 0.05)
end
plot(plt_bfcn, xticks = [-3.5, 0, 1, 4], xflip = true)
png("bifurcation_plot.png")