PredefinedDynamicalSystems.jl

PredefinedDynamicalSystems.jl

Module which contains pre-defined dynamical systems that can be used by the DynamicalSystems.jl library. To install it, run import Pkg; Pkg.add("PredefinedDynamicalSystems").

Predefined systems exist as functions that return a DynamicalSystem instance. They are accessed like:

ds = PredefinedDynamicalSystems.lorenz(u0; ρ = 32.0)

The alias Systems is also exported as a deprecation.

This module is provided purely as a convenience. It does not have any actual tests, and it is not guaranteed to be stable in future versions. It is not recommended to use this module for anything else besides on-the-spot demonstrative examples.

For some systems, a Jacobian function is also defined. The naming convention for the Jacobian function is \$(name)_jacob. So, for the above example we have J = Systems.lorenz_jacob.

All available systems are provided in the documentation, which you can either find online or build locally by running the docs/make.jl file.

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 ^{[Datseris2019]}.

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.

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. ^[Chua1992]

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. ^[Chua2007]

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

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.

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

The standard forced damped pendulum with a sine response force. duffing oscillator, that satisfies the equation

\[\ddot{x} + d \dot{x} + \sin(x) = 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.

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 ^{[Gissinger2012]}, 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.

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.

guckenheimer_holmes(u0=[-0.55582369,0.05181624,0.37766104];
    a = 0.4,
    b = 20.25,
    c = 3,
    d = 1.6,
    e = 1.7,
    f = 0.44)

\[\begin{aligned} \dot{x} &= ax - by + czx + dz(x^2 + y^2)\\ \dot{y} &= ay + bx + czy\\ \dot{z} &= e - z^2 - f(x^2 + y^2) - az^3 \end{aligned}\]

A nonlinear oscillator ^{[GuckenheimerHolmes1983]}.

halvorsen(u0=[-8.6807408,-2.4741399,0.070775762]; a = 1.4, b = 4.0)

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

An algebraically-simple chaotic system with quadratic nonlinearity ^[Sprott2010].

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 ^{[HénonHeiles1964]} 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.

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.

hindmarshrose_two_coupled(u0=[0.1, 0.2, 0.3, 0.4, 0.5, 0.6];
a = 1.0, b = 3.0, d = 5.0, r = 0.001, s = 4.0, xr = -1.6, I = 4.0,
k1 = -0.17, k2 = -0.17, k_el = 0.0, xv = 2.0)

\[\begin{aligned} \dot x_{i} = y_{i} + bx^{2}_{i} - ax^{3}_{i} - z_{i} + I - k_{i}(x_{i} - v_{s})\Gamma(x_{j}) + k(x_{j} - x_{i})\\ \dot y_{i} = c - d x^{2}_{i} - y_{i}\\ \dot z_{i} = r[s(x_{i} - x_{R}) - z_{i}]\\ \i,j=1,2 (i\neq j).\\ \end{aligned}\]

The two coupled Hindmarsh Rose element by chemical and electrical synapse. it is modelling the dynamics of a neuron's membrane potential. The default parameter values are taken from article "Dragon-king-like extreme events in coupled bursting neurons", DOI:https://doi.org/10.1103/PhysRevE.97.062311.

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.

function hyper_bao(u0 = [5.0, 8.0, 12.0, 21.0];
    a = 36.0,
    b = 3.0,
    c = 20.5,
    d = 0.1,
    k = 21.0)

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

A system showchasing hyperchaos obtained from the Lu system[^Bo-Cheng2008].

[^Bo-Cheng2008]: Bo-Cheng, B., & Zhong, L. (2008). A hyperchaotic attractor coined from chaotic Lü system. Chinese Physics Letters, 25(7), 2396.

function hyper_cai(u0 = [1.0, 1.0, 20.0, 10.0];
    a = 27.5,
    b = 3.0,
    c = 19.3,
    d = 2.9,
    e = 3.3)

\[\begin{aligned} \dot{x} &= a (y - x)\\ \dot{y} &= b x + c y - x z + w\\ \dot{z} &= -d z + y^2\\ \dot{w} &= -e x \end{aligned}\]

A system showchasing hyperchaos obtained from the Finance system^[Cai2007].

function hyper_jha(u0 = [0.1, 0.1, 0.1, 0.1];
    a = 10.0,
    b = 28.0,
    c = 8/3,
    d = 1.3)

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

An extension of the Lorenz system showchasing hyperchaos^{[Hussain2015]}.

function hyper_lorenz(u0 = [-10.0, -6.0, 0.0, 10.0];
    a = 10.0,
    b = 28.0,
    c = 8/3,
    d = -1.0)

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

An extension of the Lorenz system showchasing hyperchaos^[Wang2008]. An hyperchaotic system is characterized by two positive Lyapunov exponents.

function hyper_lu(u0 = [5.0, 8.0, 12.0, 21.0];
    a = 36,
    b = 3.0,
    c = 20.0,
    d = 1.3)

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

A system showchasing hyperchaos obtained from the Lu system^[Chen2006].

function hyper_pang(u0 = [1.0, 1.0, 10.0, 1.0];
    a = 36,
    b = 3.0,
    c = 20.0,
    d = 2.0)

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

A system showchasing hyperchaos obtained from the Lu system^[Pang2011].

function hyper_qi(u0 = [10.0, 15.0, 20.0, 22.0];
    a = 50.0,
    b = 24.0,
    c = 13,
    d = 8,
    e = 33,
    f = 30)

\[\begin{aligned} \dot{x} &= a*(y - x) + y*z\\ \dot{y} &= b*(x + y) - xz\\ \dot{z} &= - c*z - e*w + x*y\\ \dot{w} &= -d*w + f*z +x*y \end{aligned}\]

A hyperchaotic dynamical systems, showcasing a wide range of different behaviors, including rich bifurcations in different directions^[Qi2008].

hyper_roessler(u0 = [-10.0, -6.0, 0.0, 10.0];
    a = 0.25,
    b = 3.0,
    c = 0.5,
    d = 0.05)

\[\begin{aligned} \dot{x} &= -y - z\\ \dot{y} &= x + a*y + w\\ \dot{z} &= b + x*z\\ \dot{w} &= -c*z + d*w \end{aligned}\]

An extension of the Rössler system showchasing hyperchaos^{[Rossler1979]}. An hyperchaotic system is characterized by two positive Lyapunov exponents.

function hyper_wang(u0 = [5.0, 1.0, 30.0, 1.0];
    a = 10.0,
    b = 40.0,
    c = 2.5,
    d = 10.6,
    e = 4.0)

\[\begin{aligned} \dot{x} &= a*(y - x)\\ \dot{y} &= -x*z + b*x + w\\ \dot{z} &= e*x^2 - c*z\\ \dot{w} &= -d*x \end{aligned}\]

An extension of the Wang system showchasing hyperchaos^[Wang2009].

function hyper_xu(u0 = [2.0, -1.0, -2.0, -10.0];
    a = 10.0,
    b = 40.0,
    c = 2.5,
    d = 2.0,
    e = 16.0)

\[\begin{aligned} \dot{x} &= a*(y - x) + w\\ \dot{y} &= b*x + e*x*z\\ \dot{z} &= - c*z - x*y\\ \dot{w} &= x*z - d*y \end{aligned}\]

A system showchasing hyperchaos^{[Letellier2007]}.

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 ^[Lorenz1963], 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.

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)

lorenz_bounded(u0=[-13.284881, -12.444334, 34.188198];
    beta = 2.667,
    r = 64.0,
    rho = 28.0,
    sigma = 10.0
)

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

Lorenz system bounded by a confining potential ^{[SprottXiong2015]}.

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.

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.

morris_lecar(u0=[0.1, 0.1]; I = 0.15, V3 = 0.1, V1 = -0.00, V2 = 0.15, V4 = 0.1,
    VCa = 1, VL = -0.5, VK = -0.7, gCa = 1.2, gK = 2, gL = 0.5, τ = 3) -> ds

The Morris-Lecar model is ubiquitously used in computational neuroscience as a simplified model for neuronal dynamics (2D), and can also be in general as an excitable system ^{[IzhikevichBook]}. It uses the formalism of the more complete Hodgkin-Huxley model (4D), and can be viewed as a simplification thereof, with variables V for the membrane potential and N for the recovery of the Potassium current. Its original parameters were obtained from experimental studies of the giant muscle fiber in the Pacific barnacle ^{[MorrisLecar1981]}. Its evolution is given by:

\[\begin{aligned} \dot{V} &= -g_{Ca} M(V) (V - V_{Ca}) - g_K N (V - V_K) - g_L (V - V_L) + I \\ \dot{N} &= (-N + G(V)) / au \\ \end{aligned}\]

with

\[\begin{aligned} M(V) = 0.5 (1 + \tanh((x-V1)/V2)) \\ G(V) = 0.5 (1 + \tanh((x-V3)/V4)) \\\]

multispecies_competition(option = 1)

A model of competition dynamics between multiple species from Huisman and Weissing^{[Huisman2001]}. It highlights fundamental unpredictability by having extreme multistability, fractal basin boundaries and transient chaos.

TODO: write here equations when we have access to the paper (not open access).

TODO: Describe initial conditions and what option 1 means.

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. ^[Hoover1995]

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

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_0^2) - D q_0 (q_0^2 + q_2^2) \\ \dot{p}_2 &= -q_2 [A + 3\sqrt{2} B q_0 + D (q_0^2 + q_2^2)] \end{aligned}\]

This dynamical rule corresponds to a Hamiltonian used in nuclear physics to study the quadrupole vibrations of the nuclear surface ^{[Eisenberg1975][Baran1998]}.

\[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 [^Micluta-Campeanu2018]. The default initial condition is chaotic.

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

riddled_basins(u0=[0.5, 0.6, 0, 0]; γ=0.05, x̄ = 1.9, f₀=2.3, ω =3.5, x₀=1, y₀=0) → ds

\[\begin{aligned} \dot{x} &= v_x, \quad \dot{y} = v_z \\ \dot{v}_x &= -\gamma v_x - [ -4x(1-x^2) +y^2] + f_0 \sin(\omega t)x_0 \\ \dot{v}_y &= -\gamma v_y - 2y (x+\bar{x}) + f_0 \sin(\omega t)y_0 \end{aligned}\]

This 5 dimensional (time-forced) dynamical system was used by Ott et al ^{[OttRiddled2014]} to analyze riddled basins of attraction. This means nearby any point of a basin of attraction of an attractor A there is a point of the basin of attraction of another attractor B.

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 ^{[Rikitake1958]} is a system that tries to model the magnetic reversal events by means of a double-disk dynamo system.

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 ^{[Rössler1976]}. 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.

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).

sakarya(u0= [-2.8976045, 3.8877978, 3.07465];
    a = 1,
    b = 1,
    m = 1
)

\[\begin{aligned} \dot{x} &= ax + y + yz\\ \dot{y} &= - xz + yz \\ \dot{z} &= - z - mxy + b \end{aligned}\]

A system presenting robust chaos that varies from single wing to double wings to four wings. Its attractor arises due to merging of two disjoint bistable attractors ^[Li2015].

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 takes 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).

swinging_atwood(u0=[0.113296,1.5707963267948966,0.10992,0.17747]; m1=1.0, m2=4.5)

\[\begin{aligned} \dot{r} &= \frac{p_r}{M+m}\\ \dot{p}_r &= -Mg + mg\cos(\theta)\\ \dot{\theta} &= \frac{p_{\theta}}{mr^2}\\ \dot{p}_{\theta} &= -mgr\sin(\theta) \end{aligned}\]

A mechanical system consisting of two swinging weights connected by ropes and pulleys. This is only chaotic when m2 is sufficiently larger than m1, and there are nonzero initial momenta ^{[Tufillaro1984]}.

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 ^[Ruelle1980] 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. ^[Sprott2010]

ulam(N = 100, u0 = cos.(1:N); ε = 0.6)

A discrete system of N unidirectionally coupled maps on a circle, with equations

\[x^{(m)}_{n+1} = f(\varepsilon x_n^{(m-1)} + (1-\varepsilon)x_n^{(m)});\quad f(x) = 2 - x^2\]

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(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.