Tutorial

Let's attempt to detect a periodic orbit of the Lorenz system. First we define the Lorenz system itself.

using PeriodicOrbits

function lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
    return CoupledODEs(lorenz_rule, u0, [σ, ρ, β]; diffeq=(abstol=1e-10, reltol=1e-10))
end
@inbounds function lorenz_rule(u, p, t)
    du1 = p[1]*(u[2]-u[1])
    du2 = u[1]*(p[2]-u[3]) - u[2]
    du3 = u[1]*u[2] - p[3]*u[3]
    return SVector{3}(du1, du2, du3)
end

ds = lorenz()

3-dimensional CoupledODEs
 deterministic: true
 discrete time: false
 in-place:      false
 dynamic rule:  lorenz_rule
 ODE solver:    Tsit5
 ODE kwargs:    (abstol = 1.0e-10, reltol = 1.0e-10)
 parameters:    [10.0, 28.0, 2.6666666666666665]
 time:          0.0
 state:         [0.0, 10.0, 0.0]

Next, we give initial guess of the location of the periodic orbit and its period.

u0_guess = SVector(1.0, 2.0, 5.0)
T_guess = 4.2
ig = InitialGuess(u0_guess, T_guess)

InitialGuess{SVector{3, Float64}, Float64}
 u0: [1.0, 2.0, 5.0]
 T:  4.2

Then we pick an appropriate algorithm that will detect the PO. In this case we can use any algorithm intended for continuous-time dynamical systems. We choose Optimized Shooting algorithm, for more information see OptimizedShooting.

alg = OptimizedShooting(Δt=0.01, n=3)

OptimizedShooting{@NamedTuple{reltol::Float64, abstol::Float64, maxiters::Int64}}(0.01, 3, (reltol = 1.0e-6, abstol = 1.0e-6, maxiters = 1000))

Finally, we combine all the pieces to find the periodic orbit.

po = periodic_orbit(ds, alg, ig)
po

PeriodicOrbit{3, Float64, Float64}
 u0:       [0.20152, 0.43088, 13.35185]
 T:        4.41777
 stable:   missing
 discrete: false
 length:   44

The closed curve of the periodic orbit can be visualized using plotting library such as Makie.

using CairoMakie

u0 = po.points[1]
T = po.T
traj, t = trajectory(ds, T, u0; Dt = 0.01)
fig = Figure()
ax = Axis3(fig[1,1], azimuth = 0.6pi, elevation= 0.1pi)
lines!(ax, traj[:, 1], traj[:, 2], traj[:, 3], color = :blue, linewidth=1.7)
scatter!(ax, u0)
fig

To ensure that the detected period is minimal, eg. it is not a multiple of the minimal period, we can use minimal_period.

minT_po = minimal_period(ds, po)
minT_po

PeriodicOrbit{3, Float64, Float64}
 u0:       [0.20152, 0.43088, 13.35185]
 T:        4.41777
 stable:   missing
 discrete: false
 length:   100

Whether two periodic orbits are equivalent up to some tolerance. Function poequal can be used.

equal = poequal(po, minT_po; dthres=1e-3, Tthres=1e-3)
"Detected periodic orbit had minimal period: $(equal)"

"Detected periodic orbit had minimal period: true"

To determine whether found periodic orbit is stable or unstable, we can apply isstable function.

po = isstable(ds, po)
"Detected periodic orbit is $(po.stable ? "stable" : "unstable")."

"Detected periodic orbit is unstable."