Source code for brainpy.datasets.chaos.chaotic_systems

# -*- coding: utf-8 -*-

import jax.numpy as jnp

from brainpy import math as bm, dyn
from brainpy.integrators import odeint, JointEq, IntegratorRunner

__all__ = [
'henon_map_series',
'logistic_map_series',
'modified_lu_chen_series',
'mackey_glass_series',

'rabinovich_fabrikant_series',
'chen_chaotic_series',
'lu_chen_chaotic_series',
'chua_chaotic_series',
'modified_chua_series',
'lorenz_series',
'modified_Lorenz_series',
'double_scroll_series',

'PWL_duffing_series',
]

[docs]def henon_map_series(num_step, a=1.4, b=0.3, inits=None, numpy_mon=False):
r"""The Hénon map time series.

The Hénon map is a discrete-time dynamical system. It is one of the
most studied examples of dynamical systems that exhibit chaotic behavior.

.. math::

\begin{split}\begin{cases}x_{n+1} = 1-a x_n^2 + y_n\\y_{n+1} = b x_n.\end{cases}\end{split}

The map depends on two parameters, a and b, which for the classical
Hénon map have values of a = 1.4 and b = 0.3. For the classical values
the Hénon map is chaotic.

References
----------

.. [1] https://en.wikipedia.org/wiki/H%C3%A9non_map
"""
if inits is None:
inits = {'x': bm.zeros(1), 'y': bm.zeros(1)}
elif isinstance(inits, dict):
assert 'x' in inits
assert 'y' in inits
inits = {'x': bm.asarray(inits['x']), 'y': bm.asarray(inits['y'])}
assert inits['x'].shape == inits['y'].shape
else:
raise ValueError(f'Please provide dict, and do not support {type(inits)}: {inits}')
map = _HenonMap(inits['x'].size, a=a, b=b)
runner = dyn.DSRunner(map, monitors=['x', 'y'], dt=1, progress_bar=False,
numpy_mon_after_run=numpy_mon)
runner.run(num_step)
return {'ts': runner.mon.ts,
'x': runner.mon.x,
'y': runner.mon.y}

[docs]def logistic_map_series(num_step, mu=3., inits=None, numpy_mon=False):
r"""The logistic map time series.

The logistic map is defined by the following equation:

.. math::

x_{n+1}=\lambda x_{n}\left(1-x_{n}\right) \quad \text { with } \quad n=0,1,2,3 \ldots

References
----------
.. [4] https://en.wikipedia.org/wiki/Logistic_map

"""
if inits is None:
inits = bm.ones(1) * 0.2
else:
inits = bm.asarray(inits)
runner = dyn.DSRunner(_LogisticMap(inits.size, mu=mu),
monitors=['x'], dt=1, progress_bar=False,
numpy_mon_after_run=numpy_mon)
runner.run(num_step)
return {'ts': runner.mon.ts, 'x': runner.mon.x}

[docs]def modified_lu_chen_series(duration, dt=0.001, a=36, c=20, b=3, d1=1, d2=0., tau=.2,
method='rk4', inits=None, numpy_mon=False):
"""Modified Lu Chen attractor.

References
----------
"""
if inits is None:
inits = {'x': bm.ones(1), 'y': bm.ones(1), 'z': bm.ones(1) * 14}
elif isinstance(inits, dict):
assert 'x' in inits
assert 'y' in inits
assert 'z' in inits
inits = {'x': bm.asarray(inits['x']),
'y': bm.asarray(inits['y']),
'z': bm.asarray(inits['z'])}
assert inits['x'].shape == inits['y'].shape == inits['z'].shape
else:
raise ValueError
eq = _ModifiedLuChenSystem(num=inits['x'].size, a=a, b=b, c=c, d1=d1, d2=d2, tau=tau, dt=dt, method=method)
eq.x[:] = inits['x']
eq.y[:] = inits['y']
eq.z[:] = inits['z']
runner = dyn.DSRunner(eq,
monitors=['x', 'y', 'z'],
dt=dt, progress_bar=False,
numpy_mon_after_run=numpy_mon)
runner.run(duration)
return {'ts': runner.mon.ts,
'x': runner.mon['x'],
'y': runner.mon['y'],
'z': runner.mon['z']}

[docs]def mackey_glass_series(duration, dt=0.1, beta=2., gamma=1., tau=2., n=9.65,
inits=None, method='rk4', seed=None,
progress_bar=False, numpy_mon=False):
"""The Mackey-Glass time series.

Its dynamics is governed by

.. math::

\frac{dP(t)}{dt} = \frac{\beta P(t - \tau)}{1 + P(t - \tau)^n} - \gamma P(t)

where $\beta = 0.2$, $\gamma = 0.1$, $n = 10$, and the time delay $\tau = 17$. $\tau$
controls the chaotic behaviour of the equations (the higher it is, the more chaotic
the timeserie becomes.)

Parameters
----------
duration: int
dt: float, int, optional
beta: float, JaxArray
gamma: float, JaxArray
tau: float, JaxArray
n: float, JaxArray
inits: optional, float, JaxArray
method: str
seed: optional, int
progress_bar: bool

Returns
-------
result: dict
The time series data which contain

References
----------

"""
if inits is None:
inits = bm.ones(1) * 1.2
elif isinstance(inits, (float, int)):
inits = bm.asarray([inits], dtype=bm.dftype())
else:
assert isinstance(inits, (bm.ndarray, jnp.ndarray))

rng = bm.random.RandomState(seed)
xdelay = bm.TimeDelay(inits, tau, dt=dt, interp_method='round')
xdelay.data.value = inits + 0.2 * (rng.random((xdelay.num_delay_step,) + inits.shape) - 0.5)

@odeint(method=method, state_delays={'x': xdelay})
def mg_eq(x, t):
xtau = xdelay(t - tau)
return beta * xtau / (1 + xtau ** n) - gamma * x

runner = IntegratorRunner(mg_eq,
inits={'x': inits},
monitors=['x'],
fun_monitors={'x(tau)': lambda t, _: xdelay(t - tau)},
progress_bar=progress_bar, dt=dt,
numpy_mon_after_run=numpy_mon)
runner.run(duration)
return {'ts': runner.mon.ts,
'x': runner.mon['x'],
'x_tau': runner.mon['x(tau)']}

[docs]def lorenz_series(duration, dt=0.001, sigma=10, beta=8 / 3, rho=28, method='rk4', inits=None,
numpy_mon=False):
"""The Lorenz system.

The Lorenz system is a system of ordinary differential equations first
studied by mathematician and meteorologist Edward Lorenz.

Returns
-------
data: dict
A dict data with the keys of ts, x, y, and z,
where ts is the history time value, x, y, z are history
values of the variable in the Lorenz system.

References
----------

"""
dx = lambda x, t, y: sigma * (y - x)
dy = lambda y, t, x, z: x * (rho - z) - y
dz = lambda z, t, x, y: x * y - beta * z
integral = odeint(JointEq([dx, dy, dz]), method=method)

return _three_variable_model(integral,
default_inits={'x': 8, 'y': 1, 'z': 1},
duration=duration, dt=dt, inits=inits,
numpy_mon=numpy_mon)

[docs]def rabinovich_fabrikant_series(duration, dt=0.001, alpha=1.1, gamma=0.803,
method='rk4', inits=None, numpy_mon=False):
"""Rabinovich-Fabrikant equations.

The Rabinovich–Fabrikant equations are a set of three coupled ordinary
differential equations exhibiting chaotic behaviour for certain values
of the parameters. They are named after Mikhail Rabinovich and Anatoly
Fabrikant, who described them in 1979.

References
----------

"""
@odeint(method=method)
def rf_eqs(x, y, z, t):
dx = y * (z - 1 + x * x) + gamma * x
dy = x * (3 * z + 1 - x * x) + gamma * y
dz = -2 * z * (alpha + x * y)
return dx, dy, dz

return _three_variable_model(rf_eqs,
default_inits={'x': -1, 'y': 0, 'z': 0.5},
duration=duration, dt=dt, inits=inits,
numpy_mon=numpy_mon)

[docs]def chen_chaotic_series(duration, dt=0.001, a=40, b=3, c=28,
method='euler', inits=None, numpy_mon=False):
"""Chen attractor.

References
----------
"""

@odeint(method=method)
def chen_system(x, y, z, t):
dx = a * (y - x)
dy = (c - a) * x - x * z + c * y
dz = x * y - b * z
return dx, dy, dz

return _three_variable_model(chen_system,
default_inits=dict(x=-0.1, y=0.5, z=-0.6),
duration=duration, dt=dt, inits=inits,
numpy_mon=numpy_mon)

[docs]def lu_chen_chaotic_series(duration, dt=0.001, a=36, c=20, b=3, u=-15.15,
method='rk4', inits=None, numpy_mon=False):
"""Lu Chen attractor.

References
----------
"""
@odeint(method=method)
def lu_chen_system(x, y, z, t):
dx = a * (y - x)
dy = x - x * z + c * y + u
dz = x * y - b * z
return dx, dy, dz

return _three_variable_model(lu_chen_system,
default_inits=dict(x=0.1, y=0.3, z=-0.6),
duration=duration, dt=dt, inits=inits,
numpy_mon=numpy_mon)

[docs]def chua_chaotic_series(duration, dt=0.001, alpha=10, beta=14.514, gamma=0, a=-1.197, b=-0.646, method='rk4',
inits=None, numpy_mon=False):
"""Chua’s system.

References
----------
"""
@odeint(method=method)
def chua_equation(x, y, z, t):
fx = b * x + 0.5 * (a - b) * (bm.abs(x + 1) - bm.abs(x - 1))
dx = alpha * (y - x) - alpha * fx
dy = x - y + z
dz = -beta * y - gamma * z
return dx, dy, dz

return _three_variable_model(chua_equation,
default_inits=dict(x=0.001, y=0, z=0.),
duration=duration, dt=dt, inits=inits,
numpy_mon=numpy_mon)

[docs]def modified_chua_series(duration, dt=0.001, alpha=10.82, beta=14.286, a=1.3, b=.11, d=0,
method='rk4', inits=None, numpy_mon=False):
"""Modified Chua chaotic attractor.

References
----------

"""

@odeint(method=method)
def modified_chua_system(x, y, z, t):
dx = alpha * (y + b * bm.sin(bm.pi * x / 2 / a + d))
dy = x - y + z
dz = -beta * y
return dx, dy, dz

return _three_variable_model(modified_chua_system,
default_inits=dict(x=1, y=1, z=0.),
duration=duration, dt=dt, inits=inits,
numpy_mon=numpy_mon)

[docs]def modified_Lorenz_series(duration, dt=0.001, a=10, b=8 / 3, c=137 / 5,
method='rk4', inits=None, numpy_mon=False):
"""Modified Lorenz chaotic system.

References
----------
"""

@odeint(method=method)
def modified_Lorenz(x, y, z, t):
temp = 3 * bm.sqrt(x * x + y * y)
dx = (-(a + 1) * x + a - c + z * y) / 3 + ((1 - a) * (x * x - y * y) + (2 * (a + c - z)) * x * y) / temp
dy = ((c - a - z) * x - (a + 1) * y) / 3 + (2 * (a - 1) * x * y + (a + c - z) * (x * x - y * y)) / temp
dz = (3 * x * x * y - y * y * y) / 2 - b * z
return dx, dy, dz

return _three_variable_model(modified_Lorenz,
default_inits=dict(x=-8, y=4, z=10),
duration=duration, dt=dt, inits=inits,
numpy_mon=numpy_mon)

[docs]def PWL_duffing_series(duration, dt=0.001, e=0.25, m0=-0.0845, m1=0.66, omega=1, i=-14,
method='rk4', inits=None, numpy_mon=False):
"""PWL Duffing chaotic attractor.

References
----------
"""

gamma = 0.14 + i / 20

@odeint(method=method)
def PWL_duffing_eq(x, y, t):
dx = y
dy = -m1 * x - (0.5 * (m0 - m1)) * (abs(x + 1) - abs(x - 1)) - e * y + gamma * bm.cos(omega * t)
return dx, dy

return _two_variable_model(PWL_duffing_eq,
default_inits=dict(x=0, y=0.),
duration=duration, dt=dt, inits=inits,
numpy_mon=numpy_mon)

[docs]def double_scroll_series(duration, dt=0.01,
R1=1.2, R2=3.44, R4=0.193, beta=11.6, Ir=2 * 2.25e-5,
method='rk4', inits=None, numpy_mon=False):
r"""Double-scroll electronic circuit attractor.

Its behavior is governed by

.. math::

{\dot{V}}_{1} ={V}_{1}/{R}_{1}-\varDelta V/{R}_{2}\,-\,2{I}_{r}\,\sinh (\beta \varDelta V),\\
\dot{{V}_{2}} =\varDelta V/{R}_{2}+2{I}_{r}\,\sinh (\beta \varDelta V)-I,\\
\dot{I} ={V}_{2}-{R}_{4}I

in dimensionless form.

References
----------
.. [1] Chang, A., Bienfang, J. C., Hall, G. M., Gardner, J. R. &
Gauthier, D. J. Stabilizing unstable steady states using
extended time-delay autosynchronization. Chaos 8, 782–790 (1998).
"""

@odeint(method=method)
def double_scroll(x, y, z, t):
delta = x - y
dV1 = x / R1 - delta / R2 - 2 * Ir * bm.sinh(beta * delta)
dV2 = delta / R2 + 2 * Ir * bm.sinh(beta * delta) - z
dI = y - R4 * z
return dV1, dV2, dI

return _three_variable_model(double_scroll,
default_inits=dict(x=0.37926545, y=0.058339, z=-0.08167691),
duration=duration, dt=dt, inits=inits,
numpy_mon=numpy_mon)

def _two_variable_model(integrator, duration, default_inits, inits=None,
args=None, dyn_args=None, dt=0.001, numpy_mon=False):
if inits is None:
inits = default_inits
elif isinstance(inits, dict):
assert 'x' in inits
assert 'y' in inits
inits = {'x': bm.asarray(inits['x']).flatten(),
'y': bm.asarray(inits['y']).flatten()}
assert inits['x'].shape == inits['y'].shape
else:
raise ValueError

runner = IntegratorRunner(integrator, monitors=['x', 'y'], inits=inits,
args=args, dyn_args=dyn_args, dt=dt, progress_bar=False,
numpy_mon_after_run=numpy_mon)
runner.run(duration)
return {'ts': runner.mon.ts,
'x': runner.mon.x,
'y': runner.mon.y}

def _three_variable_model(integrator, duration, default_inits, inits=None, args=None,
dyn_args=None, dt=0.001, numpy_mon=False):
if inits is None:
inits = default_inits  # {'x': -1, 'y': 0, 'z': 0.5}
elif isinstance(inits, dict):
assert 'x' in inits
assert 'y' in inits
assert 'z' in inits
inits = {'x': bm.asarray(inits['x']).flatten(),
'y': bm.asarray(inits['y']).flatten(),
'z': bm.asarray(inits['z']).flatten()}
assert inits['x'].shape == inits['y'].shape == inits['z'].shape
else:
raise ValueError

runner = IntegratorRunner(integrator, monitors=['x', 'y', 'z'], inits=inits,
args=args, dyn_args=dyn_args, dt=dt, progress_bar=False,
numpy_mon_after_run=numpy_mon)
runner.run(duration)
return {'ts': runner.mon.ts,
'x': runner.mon.x,
'y': runner.mon.y,
'z': runner.mon.z}

class _HenonMap(dyn.DynamicalSystem):
"""Hénon map."""

def __init__(self, num, a=1.4, b=0.3):
super(_HenonMap, self).__init__()

# parameters
self.a = a
self.b = b
self.num = num

# variables
self.x = bm.Variable(bm.zeros(num))
self.y = bm.Variable(bm.zeros(num))

def update(self, t, dt):
x_new = 1 - self.a * self.x * self.x + self.y
self.y.value = self.b * self.x
self.x.value = x_new

class _LogisticMap(dyn.DynamicalSystem):
def __init__(self, num, mu=3.):
super(_LogisticMap, self).__init__()

self.mu = mu
self.x = bm.Variable(bm.ones(num) * 0.2)

def update(self, t, dt):
self.x.value = self.mu * self.x * (1 - self.x)

class _ModifiedLuChenSystem(dyn.DynamicalSystem):
def __init__(self, num, a=35, b=3, c=28, d0=1, d1=1, d2=0., tau=.2, dt=0.1, method='rk4'):
super(_ModifiedLuChenSystem, self).__init__()

# parameters
self.a = a
self.b = b
self.c = c
self.d0 = d0
self.d1 = d1
self.d2 = d2
self.tau = tau

# variables
self.z_delay = dyn.ConstantDelay(num, delay=tau, dt=dt)
self.z_delay.data[:] = 14
self.z = bm.Variable(self.z_delay.latest)
self.x = bm.Variable(bm.ones(num))
self.y = bm.Variable(bm.ones(num))

# functions
def derivative(x, y, z, t):
dx = self.a * (y - x)
f = self.d0 * z + self.d1 * self.z_delay.pull() - self.d2 * bm.sin(self.z_delay.pull())
dy = (self.c - self.a) * x - x * f + self.c * y
dz = x * y - self.b * z
return dx, dy, dz

self.integral = odeint(derivative, method=method)

def update(self, t, dt):
self.x.value, self.y.value, self.z.value = self.integral(
self.x, self.y, self.z, t, dt)
self.z_delay.push(self.z)
self.z_delay.update(t, dt)