using DynamicalSystems
using CairoMakie

Example: Lorenz96

Inspired by juliadynamics official doc

A Lorenz 96 model is a dynamical system that shows chaotic behaviours:

\[ \frac{dx_i}{dt} = (x_{i+1} - x_{i-2})x_{i-1} - x_i + F, \]

for $i \in \{1, \cdots, N\}$ and $N + j = j$.

We have $N$ coupled ODEs, and therefore we use thr CoupledODEs function provided by Julia.

# We first make the dynamic rule function
# use ! to define in-place, usually suggested for high dimensional systems for efficiency. 
function lorenz96_rule!(du, u, p, t) # note that t is not explicitly used in this function, but will be in CoupledODEs
    F = p[1]; N = length(u) #u is the state, an array container that contains the variables and the starting state of the system

    # boundary conditions (for the closed loop)
    du[1] = (u[2] - u[N - 1]) * u[N] - u[1] + F 
    du[2] = (u[3] - u[N]) * u[1] - u[2] + F 
    du[N] = (u[1] - u[N - 2]) * u[N - 1] - u[N] + F

    # the general condition
    for n in 3:(N - 1)
        du[n] = (u[n + 1] - u[n - 2]) * u[n - 1] - u[n] + F
        end
        return nothing # always return nothing for in place 
    end

lorenz96_rule! (generic function with 1 method)

# We define the initial state and parameters
N = 6 # number of variables 
u0 = range(0.1, 1; length = N) # the state space
p0 = [17.0] # the force F 

lorenz96_1 = CoupledODEs(lorenz96_rule!, u0, p0)
lorenz96_2 = CoupledODEs(lorenz96_rule!, u0, 3) # another example with p0 = 3

6-dimensional CoupledODEs
 deterministic: true
 discrete time: false
 in-place:      true
 dynamic rule:  lorenz96_rule!
 ODE solver:    Tsit5
 ODE kwargs:    (abstol = 1.0e-6, reltol = 1.0e-6)
 parameters:    3
 time:          0.0
 state:         [0.1, 0.28, 0.46, 0.64, 0.82, 1.0]

# We then get the data points (trajectory) 
total_time = 12.5 
sampling_time = 0.02
Y, t = trajectory(lorenz96_1, total_time; Ttr = 2.2, Δt = sampling_time)
Y2, t2 = trajectory(lorenz96_2, total_time; Ttr = 2.2, Δt = sampling_time)

(6-dimensional StateSpaceSet{Float64} with 626 points, 2.2:0.02:14.7)

fig = Figure(size = (1000, 200))
F1 = 17
F2 = 3

ax1 = Axis(fig[1, 1]; xlabel = "time", ylabel = "variable", title = L"F = %$F1")
ax2 = Axis(fig[1, 2]; xlabel = "time", title = L"F = %$F2")

for var in columns(Y)
    lines!(ax1, t, var)
end
for var in columns(Y2)
    lines!(ax2, t2, var)
end

fig

---