Predefined Dynamical Systems

antidots([u]; B = 1.0, d0 = 0.3, c = 0.2)

An antidot "superlattice" is a Hamiltonian system that corresponds to a smoothened periodic Sinai billiard with disk diameter d0 and smooth factor c [1].

This version is the two dimensional classical form of the system, with quadratic dynamical rule and a perpendicular magnetic field. Notice that the dynamical rule is with respect to the velocity instead of momentum, i.e.:

\[\begin{aligned} \dot{x} &= v_x \\ \dot{y} &= v_y \\ \dot{v_x} &= B v_y - U_x \\ \dot{v_y} &= -B v_x - U_y \\ \end{aligned}\]

with $U$ the potential energy:

\[U = \left(\tfrac{1}{c^4}\right) \left[\tfrac{d_0}{2} + c - r_a\right]^4\]

if $r_a = \sqrt{(x \mod 1)^2 + (y \mod 1)^2} < \frac{d_0}{2} + c$ and 0 otherwise. That is, the potential is periodic with period 1 in both $x, y$ and normalized such that for energy value of 1 it is a circle of diameter $d_0$. The magnetic field is also normalized such that for value B = 1 the cyclotron diameter is 1.

[1] : G. Datseris et al, New Journal of Physics 2019

arnoldcat(u0 = [0.001245, 0.00875])

\[f(x,y) = (2x+y,x+y) \mod 1\]

Arnold's cat map. A chaotic map from the torus into itself, used by Vladimir Arnold in the 1960s. [1]

[1] : Arnol'd, V. I., & Avez, A. (1968). Ergodic problems of classical mechanics.

betatransformationmap(u0 = 0.25; β=2.0)-> ds

The beta transformation, also called the generalized Bernoulli map, or the βx map, is described by

\[\begin{aligned} x_{n+1} = \beta x (\mod 1). \end{aligned}\]

The parameter β controls the dynamics of the map. Its Lyapunov exponent can be analytically shown to be λ = ln(β) ^[Ott2002]. At β=2, it becomes the dyadic transformation, also known as the bit shift map, the 2x mod 1 map, the Bernoulli map or the sawtooth map. The typical trajectory for this case is chaotic, though there are countably infinite periodic orbits ^[Ott2002].

chua(u0 = [0.7, 0.0, 0.0]; a = 15.6, b = 25.58, m0 = -8/7, m1 = -5/7)

\[\begin{aligned} \dot{x} &= a [y - h(x)]\\ \dot{y} &= x - y+z \\ \dot{z} &= b y \end{aligned}\]

where $h(x)$ is defined by

\[h(x) = m_1 x + \frac 1 2 (m_0 - m_1)(|x + 1| - |x - 1|)\]

This is a 3D continuous system that exhibits chaos.

Chua designed an electronic circuit with the expressed goal of exhibiting chaotic motion, and this system is obtained by rescaling the circuit units to simplify the form of the equation. [1]

The parameters are $a$, $b$, $m_0$, and $m_1$. Setting $a = 15.6$, $m_0 = -8/7$ and $m_1 = -5/7$, and varying the parameter $b$ from $b = 25$ to $b = 51$, one observes a classic period-doubling bifurcation route to chaos. [2]

The parameter container has the parameters in the same order as stated in this function's documentation string.

[1] : Chua, Leon O. "The genesis of Chua's circuit", 1992.

[2] : Leon O. Chua (2007) "Chua circuit", Scholarpedia, 2(10):1488.

coupled_roessler(u0=[1, -2, 0, 0.11, 0.2, 0.1];
ω1 = 0.18, ω2 = 0.22, a = 0.2, b = 0.2, c = 5.7, k1 = 0.115, k2 = 0.0)

Two coupled Rössler oscillators, used frequently in the study of chaotic synchronization. The parameter container has the parameters in the same order as stated in this function's documentation string. The equations are:

\[\begin{aligned} \dot{x_1} &= -\omega_1 y_1-z_1 \\ \dot{y_1} &= \omega_1 x+ay_1 + k_1(y_2 - y_1) \\ \dot{z_1} &= b + z_1(x_1-c) \\ \dot{x_2} &= -\omega_2 y_2-z_2 \\ \dot{y_2} &= \omega_2 x+ay_2 + k_2(y_1 - y_2) \\ \dot{z_2} &= b + z_2(x_2-c) \\ \end{aligned}\]

coupledstandardmaps(M::Int, u0 = 0.001rand(2M); ks = ones(M), Γ = 1.0)

\[\begin{aligned} \theta_{i}' &= \theta_i + p_{i}' \\ p_{i}' &= p_i + k_i\sin(\theta_i) - \Gamma \left[ \sin(\theta_{i+1} - \theta_{i}) + \sin(\theta_{i-1} - \theta_{i}) \right] \end{aligned}\]

A discrete system of M nonlinearly coupled standard maps, first introduced in [1] to study diffusion and chaos thresholds. The total dimension of the system is 2M. The maps are coupled through Γ and the i-th map has a nonlinear parameter ks[i]. The first M parameters are the ks, the M+1th parameter is Γ.

The first M entries of the state are the angles, the last M are the momenta.

[1] : H. Kantz & P. Grassberger, J. Phys. A 21, pp 127–133 (1988)

double_pendulum(u0 = [π/2, 0, 0, 0.5];
                G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)

Famous chaotic double pendulum system (also used for our logo!). Keywords are gravity (G), lengths of each rod (L1 and L2) and mass of each ball (M1 and M2). Everything is assumed in SI units.

The variables order is $[θ₁, ω₁, θ₂, ω₂]$ and they satisfy:

\[\begin{aligned} θ̇₁ &= ω₁ \\ ω̇₁ &= [M₂ L₁ ω₁² \sin φ \cos φ + M₂ G \sin θ₂ \cos φ + M₂ L₂ ω₂² \sin φ - (M₁ + M₂) G \sin θ₁] / (L₁ Δ) \\ θ̇₂ &= ω₂ \\ ω̇₂ &= [-M₂ L₂ ω₂² \sin φ \cos φ + (M₁ + M₂) G \sin θ₁ \cos φ - (M₁ + M₂) L₁ ω₁² \sin φ - (M₁ + M₂) G \sin Θ₂] / (L₂ Δ) \end{aligned}\]

where $φ = θ₂-θ₁$ and $Δ = (M₁ + M₂) - M₂ \cos² φ$.

Jacobian is created automatically (thus methods that use the Jacobian will be slower)!

(please contribute the Jacobian in LaTeX :smile:)

The parameter container has the parameters in the same order as stated in this function's documentation string.

duffing(u0 = [0.1, 0.25]; ω = 2.2, f = 27.0, d = 0.2, β = 1)

The (forced) duffing oscillator, that satisfies the equation

\[\ddot{x} + d \dot{x} + β x + x^3 = f \cos(\omega t)\]

with f, ω the forcing strength and frequency and d the damping.

The parameter container has the parameters in the same order as stated in this function's documentation string.

fitzhugh_nagumo(u = 0.5ones(2); a=3.0, b=0.2, ε=0.01, I=0.0)

Famous excitable system which emulates the firing of a neuron, with rule

\[\begin{aligned} \dot{v} &= av(v-b)(1-v) - w + I \\ \dot{w} &= \varepsilon(v - w) \end{aligned}\]

More details in the Scholarpedia entry.

gissinger(u0 = [3, 0.5, 1.5]; μ = 0.119, ν = 0.1, Γ = 0.9)

\[\begin{aligned} \dot{Q} &= \mu Q - VD \\ \dot{D} &= -\nu D + VQ \\ \dot{V} &= \Gamma -V + QD \end{aligned}\]

A continuous system that models chaotic reversals due to Gissinger [1], applied to study the reversals of the magnetic field of the Earth.

The parameter container has the parameters in the same order as stated in this function's documentation string.

[1] : C. Gissinger, Eur. Phys. J. B 85, 4, pp 1-12 (2012)

grebogi_map(u0 = [0.2, 0.]; a = 1.32, b=0.9, J₀=0.3)

\[\begin{aligned} \theta_{n+1} &= \theta_n + a\sin 2 \theta_n -b \sin 4 \theta_n -x_n\sin \theta_n\\ x_{n+1} &= -J_0 \cos \theta_n \end{aligned}\]

This map has two fixed point at (0,-J_0) and (π,J_0) which are attracting for |1+2a-4b|<1. There is a chaotic transient dynamics before the dynamical systems settles at a fixed point. This map illustrate the fractalization of the basins boundary and its uncertainty exponent α is roughly 0.2.

Example

ds = Systems.grebogi_map(rand(2))
integ = integrator(ds)
θg = range(0, 2π, length=300)
xg = range(-0.5, 0.5, length=300)
bsn, att = basins_map2D(θg, xg, integ)

henon(u0=zeros(2); a = 1.4, b = 0.3)

\[\begin{aligned} x_{n+1} &= 1 - ax^2_n+y_n \\ y_{n+1} & = bx_n \end{aligned}\]

The Hénon map is a two-dimensional mapping due to Hénon [1] that can display a strange attractor (at the default parameters). In addition, it also displays many other aspects of chaos, like period doubling or intermittency, for other parameters.

According to the author, it is a system displaying all the properties of the Lorentz system (1963) while being as simple as possible. Default values are the ones used in the original paper.

The parameter container has the parameters in the same order as stated in this function's documentation string.

[1] : M. Hénon, Commun.Math. Phys. 50, pp 69 (1976)

henonheiles(u0=[0, -0.25, 0.42081,0])

\[\begin{aligned} \dot{x} &= p_x \\ \dot{y} &= p_y \\ \dot{p}_x &= -x -2 xy \\ \dot{p}_y &= -y - (x^2 - y^2) \end{aligned}\]

The Hénon–Heiles system [1] is a conservative dynamical system and was introduced as a simplification of the motion of a star around a galactic center. It was originally intended to study the existence of a "third integral of motion" (which would make this 4D system integrable). In that search, the authors encountered chaos, as the third integral existed for only but a few initial conditions.

The default initial condition is a typical chaotic orbit. The function Systems.henonheiles_ics(E, n) generates a grid of n×n initial conditions, all having the same energy E.

[1] : Hénon, M. & Heiles, C., The Astronomical Journal 69, pp 73–79 (1964)

hindmarshrose(u0=[-1.0, 0.0, 0.0]; a=1, b=3, c=1, d=5, r=0.001, s=4, xr=-8/5, I=2.0) -> ds

\[\begin{aligned} \dot{x} &= y - ax^3 + bx^2 +I - z, \\ \dot{y} &= c - dx^2 -y, \\ \dot{z} &= r(s(x - x_r) - z) \end{aligned}\]

The Hindmarsh-Rose model reproduces the bursting behavior of a neuron's membrane potential, characterized by a fast sequence of spikes followed by a quiescent period. The x variable describes the membane potential, whose behavior can be controlled by the applied current I; the y variable describes the sodium and potassium ionic currents, and z describes an adaptation current ^{[HindmarshRose1984]}.

The default parameter values are taken from ^{[HindmarshRose1984]}, chosen to lead to periodic bursting.

^{[HindmarshRose1984]} : J. L. Hindmarsh and R. M. Rose (1984) "A model of neuronal bursting using three coupled first order differential equations", Proc. R. Soc. Lond. B 221, 87-102.

hodgkinhuxley(u0=[-60.0, 0.0, 0.0, 0.0]; I = 12.0, Vna = 50.0, Vk = -77.0, Vl = -54.4, gna = 120.0,gk = 36.0, gl = 0.3) -> ds

\[\begin{aligned} C_m \frac{dV_m}{dt} = -\overline{g}_\mathrm{K} n^4 (V_m - V_\mathrm{K}) - \overline{g}_\mathrm{Na} m^3 h(V_m - V_\mathrm{Na}) - \overline{g}_l (V_m - Vl) + I\\ \dot{n} &= \alpha_n(V_m)(1-n) - \beta_n(V_m)n \\ \dot{m} &= \alpha_m(V_m)(1-m) - \beta_m(V_m)m \\ \dot{h} &= \alpha_h(V_m)(1-h) - \beta_h(V_m)h \\ \alpha_n(V_m) = \frac{0.01(V+55)}{1 - \exp(\frac{1V+55}{10})} \quad \alpha_m(V_m) = \frac{0.1(V+40)}{1 - \exp(\frac{V+40}{10})} \quad \alpha_h(V_m) = 0.07 \exp(-\frac{(V+65)}{20}) \\ \beta_n(V_m) = 0.125 \exp(-\frac{V+65}{80}) \quad \beta_m(V_m) = 4 \exp(-\frac{V+65}{18}) \quad \beta_h(V_m) = \frac{1}{1 + \exp(-\frac{V+35}{10})} \end{aligned}\]

The Nobel-winning four-dimensional dynamical system due to Hodgkin and Huxley ^{[HodgkinHuxley1952]}, which describes the electrical spiking activity (action potentials) in neurons. A complete description of all parameters and variables is given in ^{[HodgkinHuxley1952]}, ^{[Ermentrout2010]}, and ^[Abbott2005]. The equations and default parameters used here are taken from ^{[Ermentrout2010][Abbott2005]}. They differ slightly from the original paper ^{[HodgkinHuxley1952]}, since they were changed to shift the resting potential to -65 mV, instead of the 0mV in the original paper.

Varying the injected current I from I = -5 to I = 12 takes the neuron from quiescent to a single spike, and to a tonic (repetitive) spiking. This is due to a subcritical Hopf bifurcation, which occurs close to I = 9.5.

^{[HodgkinHuxley1952]} : A. L. Hodgkin, A.F. Huxley J. Physiol., pp. 500-544 (1952).

^{[Ermentrout2010]} : G. Bard Ermentrout, and David H. Terman, "Mathematical Foundations of Neuroscience", Springer (2010).

^[Abbott2005] : L. F. Abbott, and P. Dayan, "Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems", MIT Press (2005).

ikedamap(u0=[1.0, 1.0]; a=1.0, b=1.0, c=0.4, d =6.0) -> ds

\[\begin{aligned} t &= c - \frac{d}{1 + x_n^2 + y_n^2} \\ x_{n+1} &= a + b(x_n \cos(t) - y\sin(t)) \\ y_{n+1} &= b(x\sin(t) + y \cos(t)) \end{aligned}\]

The Ikeda map was proposed by Ikeda as a model to explain the propagation of light into a ring cavity ^{[Skiadas2008]}. It generates a variety of nice-looking, interesting attractors. The default parameters are chosen to give a unique chaotic attractor. A double attractor can be obtained with parameters [a,b,c,d] = [6, 0.9, 3.1, 6], and a triple attractor can be obtained with [a,b,c,d] = [6, 9, 2.22, 6] ^{[Skiadas2008]}.

^{[Skiadas2008]} : "Chaotic Modelling and Simulation: Analysis of Chaotic Models, Attractors and Forms", CRC Press (2008).

kuramoto(D = 20, u0 = range(0, 2π; length = D);
    K = 0.3, ω = range(-1, 1; length = D)
)

The Kuramoto model^{[Kuramoto1975]} of D coupled oscillators with equation

\dot{\phi}_i = \omega_i + \frac{K}{D}\sum_{j=1}^{D} \sin(\phi_j - \phi_i)

logistic(x0 = 0.4; r = 4.0)

\[x_{n+1} = rx_n(1-x_n)\]

The logistic map is an one dimensional unimodal mapping due to May [1] and is used by many as the archetypal example of how chaos can arise from very simple equations.

Originally intentend to be a discretized model of polulation dynamics, it is now famous for its bifurcation diagram, an immensely complex graph that that was shown be universal by Feigenbaum [2].

The parameter container has the parameters in the same order as stated in this function's documentation string.

[1] : R. M. May, Nature 261, pp 459 (1976)

[2] : M. J. Feigenbaum, J. Stat. Phys. 19, pp 25 (1978)

lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3) -> ds

\[\begin{aligned} \dot{X} &= \sigma(Y-X) \\ \dot{Y} &= -XZ + \rho X -Y \\ \dot{Z} &= XY - \beta Z \end{aligned}\]

The famous three dimensional system due to Lorenz [1], shown to exhibit so-called "deterministic nonperiodic flow". It was originally invented to study a simplified form of atmospheric convection.

Currently, it is most famous for its strange attractor (occuring at the default parameters), which resembles a butterfly. For the same reason it is also associated with the term "butterfly effect" (a term which Lorenz himself disliked) even though the effect applies generally to dynamical systems. Default values are the ones used in the original paper.

The parameter container has the parameters in the same order as stated in this function's documentation string.

[1] : E. N. Lorenz, J. atmos. Sci. 20, pp 130 (1963)

lorenz84(u = [0.1, 0.1, 0.1]; F=6.846, G=1.287, a=0.25, b=4.0)

Lorenz-84's low order atmospheric general circulation model

\[\begin{aligned} \dot x = − y^2 − z^2 − ax + aF, \\ \dot y = xy − y − bxz + G, \\ \dot z = bxy + xz − z. \\ \end{aligned}\]

This system has interesting multistability property in the phase space. For the default parameter set we have four coexisting attractors that gives birth to interesting fractalized phase space as shown in ^[Freire2008]. One can see this by doing:

ds = Systems.lorenz84(rand(3))
xg = yg = range(-1.0, 2.0; length=300)
zg = range(-1.5, 1.5; length=30)
bsn, att = basins_of_attraction((xg, yg, zg), ds; mx_chk_att=4)

lorenz96(N::Int, u0 = rand(M); F=0.01)

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

N is the chain length, F the forcing. Jacobian is created automatically. (parameter container only contains F)

lorenzdl(u = [0.1, 0.1, 0.1]; R=4.7)

Diffusionless Lorenz system: it is probably the simplest rotationnaly invariant chaotic flow.

\[\begin{aligned} \dot x = y − x, \\ \dot y = -xz, \\ \dot z = xy - R. \\ \end{aligned}\]

For R=4.7 this system has two coexisting Malasoma strange attractors that are linked together as shown in ^[Sprott2014]. The fractal boundary between the basins of attractor can be visualized with a Poincaré section at z=0:

ds = Systems.lorenzdl()
xg = yg = range(-10.0, 10.0; length=300)
pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2)
bsn, att = basins_of_attraction((xg, yg), pmap)

lotkavolterra(u0=[10.0, 5.0]; α = 1.5, β = 1, δ=1, γ=3) -> ds

\[\begin{aligned} \dot{x} &= \alpha x - \beta xy, \\ \dot{y} &= \delta xy - \gamma y \end{aligned}\]

The famous Lotka-Volterra model is a simple ecological model describing the interaction between a predator and a prey species (or also parasite and host species). It has been used independently in fields such as epidemics, ecology, and economics ^{[Hoppensteadt2006]}, and is not to be confused with the Competitive Lotka-Volterra model, which describes competitive interactions between species.

The x variable describes the number of prey, while y describes the number of predator. The default parameters are taken from ^[Weisstein], which lead to typical periodic oscillations.

^{[Hoppensteadt2006]} : Frank Hoppensteadt (2006) "Predator-prey model", Scholarpedia, 1(10):1563.

^[Weisstein] : Weisstein, Eric W., "Lotka-Volterra Equations." From MathWorld–A Wolfram Web Resource. https://mathworld.wolfram.com/Lotka-VolterraEquations.html

magnetic_pendulum(u=[0.7,0.7,0,0]; d=0.3, α=0.2, ω=0.5, N=3, γs=fill(1.0,N))

Create a pangetic pendulum with N magnetics, equally distributed along the unit circle, with dynamical rule

\[\begin{aligned} \ddot{x} &= -\omega ^2x - \alpha \dot{x} - \sum_{i=1}^N \frac{\gamma_i (x - x_i)}{D_i^3} \\ \ddot{y} &= -\omega ^2y - \alpha \dot{y} - \sum_{i=1}^N \frac{\gamma_i (y - y_i)}{D_i^3} \\ D_i &= \sqrt{(x-x_i)^2 + (y-y_i)^2 + d^2} \end{aligned}\]

where α is friction, ω is eigenfrequency, d is distance of pendulum from the magnet's plane and γ is the magnetic strength.

manneville_simple(x0 = 0.4; ε = 1.1)

\[x_{n+1} = [ (1+\varepsilon)x_n + (1-\varepsilon)x_n^2 ] \mod 1\]

A simple 1D map due to Mannevile^{[Manneville1980]} that is useful in illustrating the concept and properties of intermittency.

The parameter container has the parameters in the same order as stated in this function's documentation string.

more_chaos_example(u = rand(3))

A three dimensional chaotic system introduced in ^[Sprott2020] with rule

\[\begin{aligned} \dot{x} &= y \\ \dot{y} &= -x - \textrm{sign}(z)y \\ \dot{z} &= y^2 - \exp(-x^2) \end{aligned}\]

It is noteworthy because its strange attractor is multifractal with fractal dimension ≈ 3.

nld_coupled_logistic_maps(D = 4, u0 = range(0, 1; length=D); λ = 1.2, k = 0.08)

A high-dimensional discrete dynamical system that couples D logistic maps with a strongly nonlinear all-to-all coupling. For the default parameters it displays several co-existing attractors. The equations are:

\[u_i' = \lambda - u_i^2 + k \sum_{j\ne i} (u_j^2 - u_i^2)\]

Here the prime $'$ denotes next state.

nosehoover(u0 = [0, 0.1, 0])

\[\begin{aligned} \dot{x} &= y \\ \dot{y} &= yz - x \\ \dot{z} &= 1 - y^2 \end{aligned}\]

Three dimensional conservative continuous system, discovered in 1984 during investigations in thermodynamical chemistry by Nosé and Hoover, then rediscovered by Sprott during an exhaustive search as an extremely simple chaotic system. [1]

See Chapter 4 of "Elegant Chaos" by J. C. Sprott. [2]

[1] : Hoover, W. G. (1995). Remark on ‘‘Some simple chaotic flows’’. Physical Review E, 51(1), 759.

[2] : Sprott, J. C. (2010). Elegant chaos: algebraically simple chaotic flows. World Scientific.

pomaeu_manneville(u0 = 0.2; z = 2.5)

The Pomeau-Manneville map is a one dimensional discrete map which is characteristic for displaying intermittency [1]. Specifically, for z > 2 the average time between chaotic bursts diverges, while for z > 2.5, the map iterates are long range correlated [2].

Notice that here we are providing the "symmetric" version:

\[x_{n+1} = \begin{cases} -4x_n + 3, & \quad x_n \in (0.5, 1] \\ x_n(1 + |2x_n|^{z-1}), & \quad |x_n| \le 0.5 \\ -4x_n - 3, & \quad x_n \in [-1, 0.5) \end{cases}\]

[1] : Manneville & Pomeau, Comm. Math. Phys. 74 (1980)

[2] : Meyer et al., New. J. Phys 20 (2019)

qbh([u0]; A=1.0, B=0.55, D=0.4)

A conservative dynamical system with rule

\[\begin{aligned} \dot{q}_0 &= A p_0 \\ \dot{q}_2 &= A p_2 \\ \dot{p}_0 &= -A q_0 -3 \frac{B}{\sqrt{2}} (q_2^2 - q_1^2) - D q_1 (q_1^2 + q_2^2) \\ \dot{p}_2 &= -q_2 [A + 3\sqrt{2} B q_1 + D (q_1^2 + q_2^2)] (x^2 - y^2) \end{aligned}\]

This dynamical rule corresponds to a Hamiltonian used in nuclear physics to study the quadrupole vibrations of the nuclear surface [1,2].

\[H(p_0, p_2, q_0, q_2) = \frac{A}{2}\left(p_0^2+p_2^2\right)+\frac{A}{2}\left(q_0^2+q_2^2\right) +\frac{B}{\sqrt{2}}q_0\left(3q_2^2-q_0^2\right) +\frac{D}{4}\left(q_0^2+q_2^2\right)^2\]

The Hamiltonian has a similar structure with the Henon-Heiles one, but it has an added fourth order term and presents a nontrivial dependence of chaoticity with the increase of energy [3]. The default initial condition is chaotic.

[1]: Eisenberg, J.M., & Greiner, W., Nuclear theory 2 rev ed. Netherlands: North-Holland pp 80 (1975)

[2]: Baran V. and Raduta A. A., International Journal of Modern Physics E, 7, pp 527–551 (1998)

[3]: Micluta-Campeanu S., Raportaru M.C., Nicolin A.I., Baran V., Rom. Rep. Phys. 70, pp 105 (2018)

rikitake(u0 = [1, 0, 0.6]; μ = 1.0, α = 1.0)

\[\begin{aligned} \dot{x} &= -\mu x +yz \\ \dot{y} &= -\mu y +x(z-\alpha) \\ \dot{z} &= 1 - xz \end{aligned}\]

Rikitake's dynamo is a system that tries to model the magnetic reversal events by means of a double-disk dynamo system.

[1] : T. Rikitake Math. Proc. Camb. Phil. Soc. 54, pp 89–105, (1958)

roessler(u0=[1, -2, 0.1]; a = 0.2, b = 0.2, c = 5.7)

\[\begin{aligned} \dot{x} &= -y-z \\ \dot{y} &= x+ay \\ \dot{z} &= b + z(x-c) \end{aligned}\]

This three-dimensional continuous system is due to Rössler [1]. It is a system that by design behaves similarly to the lorenz system and displays a (fractal) strange attractor. However, it is easier to analyze qualitatively, as for example the attractor is composed of a single manifold. Default values are the same as the original paper.

The parameter container has the parameters in the same order as stated in this function's documentation string.

[1] : O. E. Rössler, Phys. Lett. 57A, pp 397 (1976)

rulkovmap(u0=[1.0, 1.0]; α=4.1, β=0.001, σ=0.001) -> ds

\[\begin{aligned} x_{n+1} &= \frac{\alpha}{1+x_n^2} + y_n \\ y_{n+1} &= y_n - \sigma x_n - \beta \end{aligned}\]

The Rulkov map is a two-dimensional phenomenological model of a neuron capable of describing spikes and bursts. It was described by Rulkov ^[Rulkov2002] and is used in studies of neural networks due to its computational advantages, being fast to run.

The parameters σ and β are generally kept at 0.001, while α is chosen to give the desired dynamics. The dynamics can be quiescent for α ∈ (0,2), spiking for α ∈ (2, 2.58), triangular bursting for α ∈ (2.58, 4), and rectangular bursting for α ∈ (4, 4.62) ^{[Rulkov2001][Cao2013]}. The default parameters are taken from ^[Rulkov2001] to lead to a rectangular bursting.

^[Rulkov2002] : "Modeling of spiking-bursting neural behavior using two-dimensional map", Phys. Rev. E 65, 041922 (2002).

^[Rulkov2001] : "Regularization of Synchronized Chaotic Bursts", Phys. Rev. Lett. 86, 183 (2001).

^[Cao2013] : H. Cao and Y Wu, "Bursting types and stable domains of Rulkov neuron network with mean field coupling", International Journal of Bifurcation and Chaos,23:1330041 (2013).

shinriki(u0 = [-2, 0, 0.2]; R1 = 22.0)

Shinriki oscillator with all other parameters (besides R1) set to constants. This is a stiff problem, be careful when choosing solvers and tolerances.

sprott_dissipative_conservative(u0 = [1.0, 0, 0]; a = 2, b = 1, c = 1)

An interesting system due to Sprott^{[Sprott2014b]} where some initial conditios such as [1.0, 0, 0] lead to quasi periodic motion on a 2-torus, while for [2.0, 0, 0] motion happens on a (dissipative) chaotic attractor.

The equations are:

\[\begin{aligned} \dot{x} &= y + axy + xz \\ \dot{y} &= 1 - 2x^2 + byz \\ \dot{z_1} &= cx - x^2 - y^2 \end{aligned}\]

In the original paper there were no parameters, which are added here for exploration purposes.

standardmap(u0=[0.001245, 0.00875]; k = 0.971635)

\[\begin{aligned} \theta_{n+1} &= \theta_n + p_{n+1} \\ p_{n+1} &= p_n + k\sin(\theta_n) \end{aligned}\]

The standard map (also known as Chirikov standard map) is a two dimensional, area-preserving chaotic mapping due to Chirikov [1]. It is one of the most studied chaotic systems and by far the most studied Hamiltonian (area-preserving) mapping.

The map corresponds to the Poincaré's surface of section of the kicked rotor system. Changing the non-linearity parameter k transitions the system from completely periodic motion, to quasi-periodic, to local chaos (mixed phase-space) and finally to global chaos.

The default parameter k is the critical parameter where the golden-ratio torus is destroyed, as was calculated by Greene [2]. The e.o.m. considers the angle variable θ to be the first, and the angular momentum p to be the second, while both variables are always taken modulo 2π (the mapping is on the [0,2π)² torus).

The parameter container has the parameters in the same order as stated in this function's documentation string.

[1] : B. V. Chirikov, Preprint N. 267, Institute of Nuclear Physics, Novosibirsk (1969)

[2] : J. M. Greene, J. Math. Phys. 20, pp 1183 (1979)

stommel_thermohaline(u = [0.3, 0.2]; η1 = 3.0, η2 = 1, η3 = 0.3)

Stommel's box model for Atlantic thermohaline circulation

\[\begin{aligned} \dot{T} &= \eta_1 - T - |T-S| T \\ \dot{S} &= \eta_2 - \eta_3S - |T-S| S \end{aligned}\]

Here $T, S$ denote the dimensionless temperature and salinity differences respectively between the boxes (polar and equatorial ocean basins) and $\eta_i$ are parameters.

stuartlandau_oscillator(u0=[1.0, 0.0]; μ=1.0, ω=1.0, b=1) -> ds

The Stuart-Landau model describes a nonlinear oscillation near a Hopf bifurcation, and was proposed by Landau in 1944 to explain the transition to turbulence in a fluid ^[Landau1944]. It can be written in cartesian coordinates as ^[Deco2017]

\[\begin{aligned} \dot{x} &= (\mu -x^2 -y^2)x - \omega y - b(x^2+y^2)y \\ \dot{y} &= (\mu -x^2 -y^2)y + \omega x + b(x^2+y^2)x \end{aligned}\]

The dynamical analysis of the system is greatly facilitated by putting it in polar coordinates, where it becomes the normal form of the supercritical Hopf bifurcation) ^{[Strogatz2015]}.

\[\begin{aligned} \dot{r} &= \mu r - r^3, \\ \dot{\theta} &= \omega +br^2 \end{aligned}\]

The parameter \mu serves as the bifurcation parameter, \omega is the frequency of infinitesimal oscillations, and b controls the dependence of the frequency on the amplitude. Increasing \mu from negative to positive generates the supercritical Hopf bifurcation, leading from a stable spiral at the origin to a stable limit cycle with radius \sqrt(\mu). ^[Landau1944] : L. D. Landau, "On the problem of turbulence, In Dokl. Akad. Nauk SSSR (Vol. 44, No. 8, pp. 339-349) (1944). ^[Deco2017] : G. Deco et al "The dynamics of resting fluctuations in the brain: metastability and its dynamical cortical core", Sci Rep 7, 3095 (2017). ^{[Strogatz2015]} : Steven H. Strogatz "Nonlinear dynamics and chaos : with applications to physics, biology, chemistry, and engineering", Boulder, CO :Westview Press, a member of the Perseus Books Group (2015).

tentmap(u0 = 0.2; μ=2) -> ds

The tent map is a piecewise linear, one-dimensional map that exhibits chaotic behavior in the interval [0,1] ^[Ott2002]. Its simplicity allows it to be geometrically interpreted as generating a streching and folding process, necessary for chaos. The equations describing it are:

\[\begin{aligned} x_{n+1} = \begin{cases} \mu x, \quad &x_n < \frac{1}{2} \\ \mu (1-x), \quad &\frac{1}{2} \leq x_n \end{cases} \end{aligned}\]

The parameter μ should be kept in the interval [0,2]. At μ=2, the tent map can be brought to the logistic map with r=4 by a change of coordinates.

^[Ott2002] : E. Ott, "Chaos in Dynamical Systems" (2nd ed.) Cambridge: Cambridge University Press (2010).

thomas_cyclical(u0 = [1.0, 0, 0]; b = 0.2)

\[\begin{aligned} \dot{x} &= \sin(y) - bx\\ \dot{y} &= \sin(z) - by\\ \dot{z} &= \sin(x) - bz \end{aligned}\]

Thomas' cyclically symmetric attractor is a 3D strange attractor originally proposed by René Thomas^[Thomas1999]. It has a simple form which is cyclically symmetric in the x,y, and z variables and can be viewed as the trajectory of a frictionally dampened particle moving in a 3D lattice of forces. For more see the Wikipedia page.

Reduces to the labyrinth system for b=0, see See discussion in Section 4.4.3 of "Elegant Chaos" by J. C. Sprott.

towel(u0 = [0.085, -0.121, 0.075])

\[\begin{aligned} x_{n+1} &= 3.8 x_n (1-x_n) -0.05 (y_n +0.35) (1-2z_n) \\ y_{n+1} &= 0.1 \left[ \left( y_n +0.35 \right)\left( 1+2z_n\right) -1 \right] \left( 1 -1.9 x_n \right) \\ z_{n+1} &= 3.78 z_n (1-z_n) + b y_n \end{aligned}\]

The folded-towel map is a hyperchaotic mapping due to Rössler [1]. It is famous for being a mapping that has the smallest possible dimensions necessary for hyperchaos, having two positive and one negative Lyapunov exponent. The name comes from the fact that when plotted looks like a folded towel, in every projection.

Default values are the ones used in the original paper.

[1] : O. E. Rössler, Phys. Lett. 71A, pp 155 (1979)

ueda(u0 = [3.0, 0]; k = 0.1, B = 12.0)

\[\ddot{x} + k \dot{x} + x^3 = B\cos{t}\]

Nonautonomous Duffing-like forced oscillation system, discovered by Ueda in

1. It is one of the first chaotic systems to be discovered.

The stroboscopic plot in the (x, ̇x) plane with period 2π creates a "broken-egg attractor" for k = 0.1 and B = 12. Figure 5 of [1] is reproduced by

using Plots
ds = Systems.ueda()
a = trajectory(ds, 2π*5e3, dt = 2π)
scatter(a[:, 1], a[:, 2], markersize = 0.5, title="Ueda attractor")

For more forced oscillation systems, see Chapter 2 of "Elegant Chaos" by J. C. Sprott. [2]

[1] : Ruelle, David, ‘Strange Attractors’, The Mathematical Intelligencer, 2.3 (1980), 126–37

[2] : Sprott, J. C. (2010). Elegant chaos: algebraically simple chaotic flows. World Scientific.

vanderpol(u0=[0.5, 0.0]; μ=1.5, F=1.2, T=10) -> ds

\[\begin{aligned} \ddot{x} -\mu (1-x^2) \dot{x} + x = F \cos(\frac{2\pi t}{T}) \end{aligned}\]

The forced van der Pol oscillator is an oscillator with a nonlinear damping term driven by a sinusoidal forcing. It was proposed by Balthasar van der Pol, in his studies of nonlinear electrical circuits used in the first radios ^{[Kanamaru2007][Strogatz2015]}. The unforced oscillator (F = 0) has stable oscillations in the form of a limit cycle with a slow buildup followed by a sudden discharge, which van der Pol called relaxation oscillations ^{[Strogatz2015][vanderpol1926]}. The forced oscillator (F > 0) also has periodic behavior for some parameters, but can additionally have chaotic behavior.

The van der Pol oscillator is a specific case of both the FitzHugh-Nagumo neural model ^{[Kanamaru2007]}. The default damping parameter is taken from ^{[Strogatz2015]} and the forcing parameters are taken from ^{[Kanamaru2007]}, which generate periodic oscillations. Setting \mu=8.53 generates chaotic oscillations.

^{[Kanamaru2007]} : Takashi Kanamaru (2007) "Van der Pol oscillator", Scholarpedia, 2(1):2202.

^{[Strogatz2015]} : Steven H. Strogatz (2015) "Nonlinear dynamics and chaos : with applications to physics, biology, chemistry, and engineering", Boulder, CO :Westview Press, a member of the Perseus Books Group.

^{[vanderpol1926]} : B. Van der Pol (1926), "On relaxation-oscillations", The London, Edinburgh and Dublin Phil. Mag. & J. of Sci., 2(7), 978–992.