Skip to content

Reference

LorenzPy package:

Python package to simulate and Measure chaotic time series.

Modules exported by this package:

  • simulations: Simulate various discrete and continuous chaotic dynamical systems.
  • measures: Measures for the chaotic dynamical systems.

simulations module:

Simulate various continuous and discrete chaotic dynamical system.

Every dynamical system is represented as a class.

The available classes are: - Lorenz63 - MackeyGlass

The system's parameters are introduced in the class's constructor.

For example when creating a system object of the Lorenz63, the Lorenz parameters, sigma, rho, beta, and the timestep dt are parsed as:

sys_obj = Lorenz63(sigma=10, rho=10, beta=5, dt=1)

Each sys_obj contains a "simulate" function. To simulate 1000 time-steps of the Lorenz63 system call:

sys_obj.simulate(1000).

The general syntax to create a trajectory of a System is given as:

trajectory = (=). simulate(time_steps, starting_point=)

Examples:

>>> import lorenzpy.simulations as sims
>>> data = sims.Lorenz63().simulate(1000)
>>> data.shape
(1000, 3)

Chen

Bases: _BaseSimFlow

Simulation class for the Chen system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
class Chen(_BaseSimFlow):
    """Simulation class for the Chen system."""

    def __init__(
        self,
        a: float = 35.0,
        b: float = 3.0,
        c: float = 28.0,
        dt: float = 0.02,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the Chen simulation object.

        Args:
            a: a parameter of Chen equation.
            b: b parameter of Chen equation.
            c: c parameter of Chen equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.a = a
        self.b = b
        self.c = c

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of Chen equation."""
        return np.array(
            [
                self.a * (x[1] - x[0]),
                (self.c - self.a) * x[0] - x[0] * x[2] + self.c * x[1],
                x[0] * x[1] - self.b * x[2],
            ]
        )

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of Chen system."""
        return np.array([-10.0, 0.0, 37.0])

__init__(a=35.0, b=3.0, c=28.0, dt=0.02, solver='rk4')

Initialize the Chen simulation object.

Parameters:

Name Type Description Default
a float

a parameter of Chen equation.

 35.0
b float

b parameter of Chen equation.

 3.0
c float

c parameter of Chen equation.

 28.0
dt float

Time step to simulate.

 0.02
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
def __init__(
    self,
    a: float = 35.0,
    b: float = 3.0,
    c: float = 28.0,
    dt: float = 0.02,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the Chen simulation object.

    Args:
        a: a parameter of Chen equation.
        b: b parameter of Chen equation.
        c: c parameter of Chen equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.a = a
    self.b = b
    self.c = c

flow(x)

Return the flow of Chen equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
147
148
149
150
151
152
153
154
155
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of Chen equation."""
    return np.array(
        [
            self.a * (x[1] - x[0]),
            (self.c - self.a) * x[0] - x[0] * x[2] + self.c * x[1],
            x[0] * x[1] - self.b * x[2],
        ]
    )

get_default_starting_pnt()

Return default starting point of Chen system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
157
158
159
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of Chen system."""
    return np.array([-10.0, 0.0, 37.0])

ChuaCircuit

Bases: _BaseSimFlow

Simulation class for the ChuaCircuit system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
class ChuaCircuit(_BaseSimFlow):
    """Simulation class for the ChuaCircuit system."""

    def __init__(
        self,
        alpha: float = 9.0,
        beta: float = 100 / 7,
        a: float = 8 / 7,
        b: float = 5 / 7,
        dt: float = 0.1,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the ChuaCircuit simulation object.

        Args:
            alpha: alpha parameter of ChuaCircuit equation.
            beta: beta parameter of ChuaCircuit equation.
            a: a parameter of ChuaCircuit equation.
            b: b parameter of ChuaCircuit equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.alpha = alpha
        self.beta = beta
        self.a = a
        self.b = b

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of ChuaCircuit equation."""
        return np.array(
            [
                self.alpha
                * (
                    x[1]
                    - x[0]
                    + self.b * x[0]
                    + 0.5 * (self.a - self.b) * (np.abs(x[0] + 1) - np.abs(x[0] - 1))
                ),
                x[0] - x[1] + x[2],
                -self.beta * x[1],
            ]
        )

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of ChuaCircuit system."""
        return np.array([0.0, 0.0, 0.6])

__init__(alpha=9.0, beta=100 / 7, a=8 / 7, b=5 / 7, dt=0.1, solver='rk4')

Initialize the ChuaCircuit simulation object.

Parameters:

Name Type Description Default
alpha float

alpha parameter of ChuaCircuit equation.

 9.0
beta float

beta parameter of ChuaCircuit equation.

 100 / 7
a float

a parameter of ChuaCircuit equation.

 8 / 7
b float

b parameter of ChuaCircuit equation.

 5 / 7
dt float

Time step to simulate.

 0.1
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
def __init__(
    self,
    alpha: float = 9.0,
    beta: float = 100 / 7,
    a: float = 8 / 7,
    b: float = 5 / 7,
    dt: float = 0.1,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the ChuaCircuit simulation object.

    Args:
        alpha: alpha parameter of ChuaCircuit equation.
        beta: beta parameter of ChuaCircuit equation.
        a: a parameter of ChuaCircuit equation.
        b: b parameter of ChuaCircuit equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.alpha = alpha
    self.beta = beta
    self.a = a
    self.b = b

flow(x)

Return the flow of ChuaCircuit equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of ChuaCircuit equation."""
    return np.array(
        [
            self.alpha
            * (
                x[1]
                - x[0]
                + self.b * x[0]
                + 0.5 * (self.a - self.b) * (np.abs(x[0] + 1) - np.abs(x[0] - 1))
            ),
            x[0] - x[1] + x[2],
            -self.beta * x[1],
        ]
    )

get_default_starting_pnt()

Return default starting point of ChuaCircuit system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
206
207
208
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of ChuaCircuit system."""
    return np.array([0.0, 0.0, 0.6])

ComplexButterfly

Bases: _BaseSimFlow

Simulation class for the ComplexButterfly system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
class ComplexButterfly(_BaseSimFlow):
    """Simulation class for the ComplexButterfly system."""

    def __init__(
        self,
        a: float = 0.55,
        dt: float = 0.1,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the ComplexButterfly simulation object.

        Args:
            a: a parameter of ComplexButterfly equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.a = a

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of ComplexButterfly equation."""
        return np.array(
            [self.a * (x[1] - x[0]), -x[2] * np.sign(x[0]), np.abs(x[0]) - 1]
        )

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of ComplexButterfly system."""
        return np.array([0.2, 0.0, 0.0])

__init__(a=0.55, dt=0.1, solver='rk4')

Initialize the ComplexButterfly simulation object.

Parameters:

Name Type Description Default
a float

a parameter of ComplexButterfly equation.

 0.55
dt float

Time step to simulate.

 0.1
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
def __init__(
    self,
    a: float = 0.55,
    dt: float = 0.1,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the ComplexButterfly simulation object.

    Args:
        a: a parameter of ComplexButterfly equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.a = a

flow(x)

Return the flow of ComplexButterfly equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
111
112
113
114
115
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of ComplexButterfly equation."""
    return np.array(
        [self.a * (x[1] - x[0]), -x[2] * np.sign(x[0]), np.abs(x[0]) - 1]
    )

get_default_starting_pnt()

Return default starting point of ComplexButterfly system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
117
118
119
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of ComplexButterfly system."""
    return np.array([0.2, 0.0, 0.0])

DoublePendulum

Bases: _BaseSimFlow

Simulation class for the dimensionless double pendulum with m1 = m2 and l1=l2.

The state space is given by [angle1, angle2, angular_vel, angular_vel2].

Source code in src\lorenzpy\simulations\autonomous_flows.py 
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
class DoublePendulum(_BaseSimFlow):
    """Simulation class for the dimensionless double pendulum with m1 = m2 and l1=l2.

    The state space is given by [angle1, angle2, angular_vel, angular_vel2].
    """

    def __init__(
        self,
        dt: float = 0.1,
        solver: str
        | str
        | Callable[[Callable, float, np.ndarray], np.ndarray] = create_scipy_ivp_solver(
            "DOP853"
        ),
    ):
        """Initialize the Doueble Pendulum simulation object.

        Args:
            dt: Time step to simulate.
            solver: The solver. Default is DOP853 scipy solver here.
        """
        super().__init__(dt, solver)

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of double pendulum."""
        angle1, angle2, angle1_dot, angle2_dot = x[0], x[1], x[2], x[3]

        delta_angle = angle1 - angle2

        angle1_dotdot = (
            9 * np.cos(delta_angle) * np.sin(delta_angle) * angle1_dot**2
            + 6 * np.sin(delta_angle) * angle2_dot**2
            + 18 * np.sin(angle1)
            - 9 * np.cos(delta_angle) * np.sin(angle2)
        ) / (9 * np.cos(delta_angle) ** 2 - 16)

        angle2_dotdot = (
            24 * np.sin(delta_angle) * angle1_dot**2
            + 9 * np.cos(delta_angle) * np.sin(delta_angle) * angle2_dot**2
            + 27 * np.cos(delta_angle) * np.sin(angle1)
            - 24 * np.sin(angle2)
        ) / (16 - 9 * np.cos(delta_angle) ** 2)

        return np.array([angle1_dot, angle2_dot, angle1_dotdot, angle2_dotdot])

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of Double Pendulum."""
        return np.array([0.6, 2.04, 0, 0])

__init__(dt=0.1, solver=create_scipy_ivp_solver('DOP853'))

Initialize the Doueble Pendulum simulation object.

Parameters:

Name Type Description Default
dt float

Time step to simulate.

 0.1
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver. Default is DOP853 scipy solver here.

 create_scipy_ivp_solver('DOP853')
Source code in src\lorenzpy\simulations\autonomous_flows.py 
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
def __init__(
    self,
    dt: float = 0.1,
    solver: str
    | str
    | Callable[[Callable, float, np.ndarray], np.ndarray] = create_scipy_ivp_solver(
        "DOP853"
    ),
):
    """Initialize the Doueble Pendulum simulation object.

    Args:
        dt: Time step to simulate.
        solver: The solver. Default is DOP853 scipy solver here.
    """
    super().__init__(dt, solver)

flow(x)

Return the flow of double pendulum.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of double pendulum."""
    angle1, angle2, angle1_dot, angle2_dot = x[0], x[1], x[2], x[3]

    delta_angle = angle1 - angle2

    angle1_dotdot = (
        9 * np.cos(delta_angle) * np.sin(delta_angle) * angle1_dot**2
        + 6 * np.sin(delta_angle) * angle2_dot**2
        + 18 * np.sin(angle1)
        - 9 * np.cos(delta_angle) * np.sin(angle2)
    ) / (9 * np.cos(delta_angle) ** 2 - 16)

    angle2_dotdot = (
        24 * np.sin(delta_angle) * angle1_dot**2
        + 9 * np.cos(delta_angle) * np.sin(delta_angle) * angle2_dot**2
        + 27 * np.cos(delta_angle) * np.sin(angle1)
        - 24 * np.sin(angle2)
    ) / (16 - 9 * np.cos(delta_angle) ** 2)

    return np.array([angle1_dot, angle2_dot, angle1_dotdot, angle2_dotdot])

get_default_starting_pnt()

Return default starting point of Double Pendulum.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
420
421
422
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of Double Pendulum."""
    return np.array([0.6, 2.04, 0, 0])

DoubleScroll

Bases: _BaseSimFlow

Simulation class for the DoubleScroll system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
class DoubleScroll(_BaseSimFlow):
    """Simulation class for the DoubleScroll system."""

    def __init__(
        self,
        a: float = 0.8,
        dt: float = 0.3,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the DoubleScroll simulation object.

        Args:
            a: a parameter of DoubleScroll equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.a = a

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of DoubleScroll equation."""
        return np.array([x[1], x[2], -self.a * (x[2] + x[1] + x[0] - np.sign(x[0]))])

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of DoubleScroll system."""
        return np.array([0.01, 0.01, 0.0])

__init__(a=0.8, dt=0.3, solver='rk4')

Initialize the DoubleScroll simulation object.

Parameters:

Name Type Description Default
a float

a parameter of DoubleScroll equation.

 0.8
dt float

Time step to simulate.

 0.3
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
def __init__(
    self,
    a: float = 0.8,
    dt: float = 0.3,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the DoubleScroll simulation object.

    Args:
        a: a parameter of DoubleScroll equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.a = a

flow(x)

Return the flow of DoubleScroll equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
366
367
368
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of DoubleScroll equation."""
    return np.array([x[1], x[2], -self.a * (x[2] + x[1] + x[0] - np.sign(x[0]))])

get_default_starting_pnt()

Return default starting point of DoubleScroll system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
370
371
372
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of DoubleScroll system."""
    return np.array([0.01, 0.01, 0.0])

Halvorsen

Bases: _BaseSimFlow

Simulation class for the Halvorsen system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
class Halvorsen(_BaseSimFlow):
    """Simulation class for the Halvorsen system."""

    def __init__(
        self,
        a: float = 1.27,
        dt: float = 0.05,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the Halvorsen simulation object.

        Args:
            a: a parameter of Halvorsen equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.a = a

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of Halvorsen equation."""
        return np.array(
            [
                -self.a * x[0] - 4 * x[1] - 4 * x[2] - x[1] ** 2,
                -self.a * x[1] - 4 * x[2] - 4 * x[0] - x[2] ** 2,
                -self.a * x[2] - 4 * x[0] - 4 * x[1] - x[0] ** 2,
            ]
        )

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of Halvorsen system."""
        return np.array([-5.0, 0.0, 0.0])

__init__(a=1.27, dt=0.05, solver='rk4')

Initialize the Halvorsen simulation object.

Parameters:

Name Type Description Default
a float

a parameter of Halvorsen equation.

 1.27
dt float

Time step to simulate.

 0.05
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
def __init__(
    self,
    a: float = 1.27,
    dt: float = 0.05,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the Halvorsen simulation object.

    Args:
        a: a parameter of Halvorsen equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.a = a

flow(x)

Return the flow of Halvorsen equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
332
333
334
335
336
337
338
339
340
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of Halvorsen equation."""
    return np.array(
        [
            -self.a * x[0] - 4 * x[1] - 4 * x[2] - x[1] ** 2,
            -self.a * x[1] - 4 * x[2] - 4 * x[0] - x[2] ** 2,
            -self.a * x[2] - 4 * x[0] - 4 * x[1] - x[0] ** 2,
        ]
    )

get_default_starting_pnt()

Return default starting point of Halvorsen system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
342
343
344
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of Halvorsen system."""
    return np.array([-5.0, 0.0, 0.0])

Henon

Bases: _BaseSimIterate

Simulate the 2-dimensional dissipative map: Henon map.

Source code in src\lorenzpy\simulations\discrete_maps.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
class Henon(_BaseSimIterate):
    """Simulate the 2-dimensional dissipative map: Henon map."""

    def __init__(self, a: float = 1.4, b: float = 0.3) -> None:
        """Initialize the Logistic Map simulation object.

        Args:
            a: a parameter of the Henon map.
            b: b parameter of the Henon map.
        """
        self.a = a
        self.b = b

    def iterate(self, x: np.ndarray) -> np.ndarray:
        """Iterate the Henon map one step."""
        return np.array([1 - self.a * x[0] ** 2 + self.b * x[1], x[0]])

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of the Henon map."""
        return np.array([0.0, 0.9])

__init__(a=1.4, b=0.3)

Initialize the Logistic Map simulation object.

Parameters:

Name Type Description Default
a float

a parameter of the Henon map.

 1.4
b float

b parameter of the Henon map.

 0.3
Source code in src\lorenzpy\simulations\discrete_maps.py 
36
37
38
39
40
41
42
43
44
def __init__(self, a: float = 1.4, b: float = 0.3) -> None:
    """Initialize the Logistic Map simulation object.

    Args:
        a: a parameter of the Henon map.
        b: b parameter of the Henon map.
    """
    self.a = a
    self.b = b

get_default_starting_pnt()

Return default starting point of the Henon map.

Source code in src\lorenzpy\simulations\discrete_maps.py 
50
51
52
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of the Henon map."""
    return np.array([0.0, 0.9])

iterate(x)

Iterate the Henon map one step.

Source code in src\lorenzpy\simulations\discrete_maps.py 
46
47
48
def iterate(self, x: np.ndarray) -> np.ndarray:
    """Iterate the Henon map one step."""
    return np.array([1 - self.a * x[0] ** 2 + self.b * x[1], x[0]])

KuramotoSivashinsky

Bases: _BaseSimIterate

Simulate the n-dimensional Kuramoto-Sivashinsky PDE.

Note: dimension must be an even number.

PDE: y_t = -yy_x - (1+eps)y_xx - y_xxxx.

Reference for the numerical integration: "fourth order time stepping for stiff pde-kassam trefethen 2005" at https://people.maths.ox.ac.uk/trefethen/publication/PDF/2005_111.pdf

Python implementation at: https://github.com/E-Renshaw/kuramoto-sivashinsky

Literature values (doi:10.1017/S1446181119000105) for Lyapunov Exponents: - lyapunov exponents: (0.080, 0.056, 0.014, 0.003, -0.003 ...) They refer to: - Parameters: {"sys_length": 36.0, "eps": 0.0}

Source code in src\lorenzpy\simulations\others.py 
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
class KuramotoSivashinsky(_BaseSimIterate):
    """Simulate the n-dimensional Kuramoto-Sivashinsky PDE.

    Note: dimension must be an even number.

    PDE: y_t = -y*y_x - (1+eps)*y_xx - y_xxxx.

    Reference for the numerical integration:
    "fourth order time stepping for stiff pde-kassam trefethen 2005" at
    https://people.maths.ox.ac.uk/trefethen/publication/PDF/2005_111.pdf

    Python implementation at: https://github.com/E-Renshaw/kuramoto-sivashinsky

    Literature values (doi:10.1017/S1446181119000105) for Lyapunov Exponents:
    - lyapunov exponents: (0.080, 0.056, 0.014, 0.003, -0.003 ...)
    They refer to:
    - Parameters: {"sys_length": 36.0, "eps": 0.0}
    """

    def __init__(
        self,
        sys_dim: int = 50,
        sys_length: float = 36.0,
        eps: float = 0.0,
        dt: float = 0.1,
    ) -> None:
        """Initialize the Kuramoto-Sivashinsky simulation object.

        Args:
            sys_dim: The dimension of the system.
            sys_length: The physical length of the system.
            eps: A parameter in front of the y_xx term.
            dt: Time step to simulate.
        """
        self.sys_dim = sys_dim
        self.sys_length = sys_length
        self.eps = eps
        self.dt = dt

        self._prepare()

    def _prepare(self) -> None:
        """Calculate internal attributes.

        TODO: make auxiliary variables protected.
        """
        k = (
            np.transpose(
                np.conj(
                    np.concatenate(
                        (
                            np.arange(0, self.sys_dim / 2),
                            np.array([0]),
                            np.arange(-self.sys_dim / 2 + 1, 0),
                        )
                    )
                )
            )
            * 2
            * np.pi
            / self.sys_length
        )

        L = (1 + self.eps) * k**2 - k**4

        self.E = np.exp(self.dt * L)
        self.E_2 = np.exp(self.dt * L / 2)
        M = 64
        r = np.exp(1j * np.pi * (np.arange(1, M + 1) - 0.5) / M)
        LR = self.dt * np.transpose(np.repeat([L], M, axis=0)) + np.repeat(
            [r], self.sys_dim, axis=0
        )
        self.Q = self.dt * np.real(np.mean((np.exp(LR / 2) - 1) / LR, axis=1))
        self.f1 = self.dt * np.real(
            np.mean((-4 - LR + np.exp(LR) * (4 - 3 * LR + LR**2)) / LR**3, axis=1)
        )
        self.f2 = self.dt * np.real(
            np.mean((2 + LR + np.exp(LR) * (-2 + LR)) / LR**3, axis=1)
        )
        self.f3 = self.dt * np.real(
            np.mean((-4 - 3 * LR - LR**2 + np.exp(LR) * (4 - LR)) / LR**3, axis=1)
        )

        self.g = -0.5j * k

    def iterate(self, x: np.ndarray) -> np.ndarray:
        """Calculate next timestep x(t+1) with given x(t).

        Args:
            x: (x_0(i),x_1(i),..) coordinates. Needs to have shape (self.sys_dim,).

        Returns:
            : (x_0(i+1),x_1(i+1),..) corresponding to input x.

        """
        v = np.fft.fft(x)
        Nv = self.g * np.fft.fft(np.real(np.fft.ifft(v)) ** 2)
        a = self.E_2 * v + self.Q * Nv
        Na = self.g * np.fft.fft(np.real(np.fft.ifft(a)) ** 2)
        b = self.E_2 * v + self.Q * Na
        Nb = self.g * np.fft.fft(np.real(np.fft.ifft(b)) ** 2)
        c = self.E_2 * a + self.Q * (2 * Nb - Nv)
        Nc = self.g * np.fft.fft(np.real(np.fft.ifft(c)) ** 2)
        v = self.E * v + Nv * self.f1 + 2 * (Na + Nb) * self.f2 + Nc * self.f3
        return np.real(np.fft.ifft(v))

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of KS system."""
        x = (
            self.sys_length
            * np.transpose(np.conj(np.arange(1, self.sys_dim + 1)))
            / self.sys_dim
        )
        return np.array(
            np.cos(2 * np.pi * x / self.sys_length)
            * (1 + np.sin(2 * np.pi * x / self.sys_length))
        )

__init__(sys_dim=50, sys_length=36.0, eps=0.0, dt=0.1)

Initialize the Kuramoto-Sivashinsky simulation object.

Parameters:

Name Type Description Default
sys_dim int

The dimension of the system.

 50
sys_length float

The physical length of the system.

 36.0
eps float

A parameter in front of the y_xx term.

 0.0
dt float

Time step to simulate.

 0.1
Source code in src\lorenzpy\simulations\others.py 
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
def __init__(
    self,
    sys_dim: int = 50,
    sys_length: float = 36.0,
    eps: float = 0.0,
    dt: float = 0.1,
) -> None:
    """Initialize the Kuramoto-Sivashinsky simulation object.

    Args:
        sys_dim: The dimension of the system.
        sys_length: The physical length of the system.
        eps: A parameter in front of the y_xx term.
        dt: Time step to simulate.
    """
    self.sys_dim = sys_dim
    self.sys_length = sys_length
    self.eps = eps
    self.dt = dt

    self._prepare()

get_default_starting_pnt()

Return default starting point of KS system.

Source code in src\lorenzpy\simulations\others.py 
118
119
120
121
122
123
124
125
126
127
128
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of KS system."""
    x = (
        self.sys_length
        * np.transpose(np.conj(np.arange(1, self.sys_dim + 1)))
        / self.sys_dim
    )
    return np.array(
        np.cos(2 * np.pi * x / self.sys_length)
        * (1 + np.sin(2 * np.pi * x / self.sys_length))
    )

iterate(x)

Calculate next timestep x(t+1) with given x(t).

Parameters:

Name Type Description Default
x ndarray

(x_0(i),x_1(i),..) coordinates. Needs to have shape (self.sys_dim,).

 required

Returns:

Type Description
ndarray

(x_0(i+1),x_1(i+1),..) corresponding to input x.
Source code in src\lorenzpy\simulations\others.py 
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
def iterate(self, x: np.ndarray) -> np.ndarray:
    """Calculate next timestep x(t+1) with given x(t).

    Args:
        x: (x_0(i),x_1(i),..) coordinates. Needs to have shape (self.sys_dim,).

    Returns:
        : (x_0(i+1),x_1(i+1),..) corresponding to input x.

    """
    v = np.fft.fft(x)
    Nv = self.g * np.fft.fft(np.real(np.fft.ifft(v)) ** 2)
    a = self.E_2 * v + self.Q * Nv
    Na = self.g * np.fft.fft(np.real(np.fft.ifft(a)) ** 2)
    b = self.E_2 * v + self.Q * Na
    Nb = self.g * np.fft.fft(np.real(np.fft.ifft(b)) ** 2)
    c = self.E_2 * a + self.Q * (2 * Nb - Nv)
    Nc = self.g * np.fft.fft(np.real(np.fft.ifft(c)) ** 2)
    v = self.E * v + Nv * self.f1 + 2 * (Na + Nb) * self.f2 + Nc * self.f3
    return np.real(np.fft.ifft(v))

Logistic

Bases: _BaseSimIterate

Simulation class for the Logistic map.

Source code in src\lorenzpy\simulations\discrete_maps.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
class Logistic(_BaseSimIterate):
    """Simulation class for the Logistic map."""

    def __init__(self, r: float = 4.0) -> None:
        """Initialize the Logistic Map simulation object.

        Args:
            r: r parameter of the logistic map.
        """
        self.r = r

    def iterate(self, x: np.ndarray) -> np.ndarray:
        """Iterate the logistic map one step."""
        return np.array(
            [
                self.r * x[0] * (1 - x[0]),
            ]
        )

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of the Logistic map."""
        return np.array([0.1])

__init__(r=4.0)

Initialize the Logistic Map simulation object.

Parameters:

Name Type Description Default
r float

r parameter of the logistic map.

 4.0
Source code in src\lorenzpy\simulations\discrete_maps.py 
12
13
14
15
16
17
18
def __init__(self, r: float = 4.0) -> None:
    """Initialize the Logistic Map simulation object.

    Args:
        r: r parameter of the logistic map.
    """
    self.r = r

get_default_starting_pnt()

Return default starting point of the Logistic map.

Source code in src\lorenzpy\simulations\discrete_maps.py 
28
29
30
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of the Logistic map."""
    return np.array([0.1])

iterate(x)

Iterate the logistic map one step.

Source code in src\lorenzpy\simulations\discrete_maps.py 
20
21
22
23
24
25
26
def iterate(self, x: np.ndarray) -> np.ndarray:
    """Iterate the logistic map one step."""
    return np.array(
        [
            self.r * x[0] * (1 - x[0]),
        ]
    )

Lorenz63

Bases: _BaseSimFlow

Simulation class for the Lorenz63 system.

This function is able to simulate the chaotic dynamical system originally introduced by Lorenz.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
class Lorenz63(_BaseSimFlow):
    """Simulation class for the Lorenz63 system.

    This function is able to simulate the chaotic dynamical system originally
    introduced by Lorenz.
    """

    def __init__(
        self,
        sigma: float = 10.0,
        rho: float = 28.0,
        beta: float = 8 / 3,
        dt: float = 0.03,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the Lorenz63 simulation object.

        Args:
            sigma: Sigma parameter of Lorenz63 equation.
            rho: Rho parameter of Lorenz63 equation.
            beta: beta parameter of Lorenz63 equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.sigma = sigma
        self.rho = rho
        self.beta = beta

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of Lorenz63 equation."""
        return np.array(
            [
                self.sigma * (x[1] - x[0]),
                x[0] * (self.rho - x[2]) - x[1],
                x[0] * x[1] - self.beta * x[2],
            ]
        )

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of Lorenz63 system."""
        return np.array([0.0, -0.01, 9.0])

__init__(sigma=10.0, rho=28.0, beta=8 / 3, dt=0.03, solver='rk4')

Initialize the Lorenz63 simulation object.

Parameters:

Name Type Description Default
sigma float

Sigma parameter of Lorenz63 equation.

 10.0
rho float

Rho parameter of Lorenz63 equation.

 28.0
beta float

beta parameter of Lorenz63 equation.

 8 / 3
dt float

Time step to simulate.

 0.03
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def __init__(
    self,
    sigma: float = 10.0,
    rho: float = 28.0,
    beta: float = 8 / 3,
    dt: float = 0.03,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the Lorenz63 simulation object.

    Args:
        sigma: Sigma parameter of Lorenz63 equation.
        rho: Rho parameter of Lorenz63 equation.
        beta: beta parameter of Lorenz63 equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.sigma = sigma
    self.rho = rho
    self.beta = beta

flow(x)

Return the flow of Lorenz63 equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
41
42
43
44
45
46
47
48
49
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of Lorenz63 equation."""
    return np.array(
        [
            self.sigma * (x[1] - x[0]),
            x[0] * (self.rho - x[2]) - x[1],
            x[0] * x[1] - self.beta * x[2],
        ]
    )

get_default_starting_pnt()

Return default starting point of Lorenz63 system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
51
52
53
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of Lorenz63 system."""
    return np.array([0.0, -0.01, 9.0])

Lorenz96

Bases: _BaseSimFlow

Simulate the n-dimensional Lorenz 96 model.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
class Lorenz96(_BaseSimFlow):
    """Simulate the n-dimensional Lorenz 96 model."""

    def __init__(
        self,
        sys_dim: int = 30,
        force: float = 8.0,
        dt: float = 0.05,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ) -> None:
        """Initialize the Lorenz96 simulation object.

        Args:
            sys_dim: The dimension of the Lorenz96 system.
            force: The force value.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt=dt, solver=solver)
        self.sys_dim = sys_dim
        self.force = force

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of Lorenz96 equation."""
        system_dimension = x.shape[0]
        derivative = np.zeros(system_dimension)
        # Periodic Boundary Conditions for the 3 edge cases i=1,2,system_dimension
        derivative[0] = (x[1] - x[system_dimension - 2]) * x[system_dimension - 1] - x[
            0
        ]
        derivative[1] = (x[2] - x[system_dimension - 1]) * x[0] - x[1]
        derivative[system_dimension - 1] = (x[0] - x[system_dimension - 3]) * x[
            system_dimension - 2
        ] - x[system_dimension - 1]

        # TODO: Rewrite using numpy vectorization to make faster
        for i in range(2, system_dimension - 1):
            derivative[i] = (x[i + 1] - x[i - 2]) * x[i - 1] - x[i]

        derivative = derivative + self.force
        return derivative

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of Lorenz96 system."""
        return np.sin(np.arange(self.sys_dim))

__init__(sys_dim=30, force=8.0, dt=0.05, solver='rk4')

Initialize the Lorenz96 simulation object.

Parameters:

Name Type Description Default
sys_dim int

The dimension of the Lorenz96 system.

 30
force float

The force value.

 8.0
dt float

Time step to simulate.

 0.05
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
def __init__(
    self,
    sys_dim: int = 30,
    force: float = 8.0,
    dt: float = 0.05,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
) -> None:
    """Initialize the Lorenz96 simulation object.

    Args:
        sys_dim: The dimension of the Lorenz96 system.
        force: The force value.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt=dt, solver=solver)
    self.sys_dim = sys_dim
    self.force = force

flow(x)

Return the flow of Lorenz96 equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of Lorenz96 equation."""
    system_dimension = x.shape[0]
    derivative = np.zeros(system_dimension)
    # Periodic Boundary Conditions for the 3 edge cases i=1,2,system_dimension
    derivative[0] = (x[1] - x[system_dimension - 2]) * x[system_dimension - 1] - x[
        0
    ]
    derivative[1] = (x[2] - x[system_dimension - 1]) * x[0] - x[1]
    derivative[system_dimension - 1] = (x[0] - x[system_dimension - 3]) * x[
        system_dimension - 2
    ] - x[system_dimension - 1]

    # TODO: Rewrite using numpy vectorization to make faster
    for i in range(2, system_dimension - 1):
        derivative[i] = (x[i + 1] - x[i - 2]) * x[i - 1] - x[i]

    derivative = derivative + self.force
    return derivative

get_default_starting_pnt()

Return default starting point of Lorenz96 system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
467
468
469
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of Lorenz96 system."""
    return np.sin(np.arange(self.sys_dim))

MackeyGlass

Bases: _BaseSim

Simulate the Mackey-Glass delay differential system.

TODO: Add literature values for Lyapunov etc. TODO: Hint the differences between this class and the other Sim classes (delay). TODO: Check if the structure is really good? TODO: Add Proper Tests. TODO: Decide whether to use the simple forward-euler or RK4-style update.

Note: As the Mackey-Glass system is a delay-differential equation, the class does not contain a simple iterate function.

Source code in src\lorenzpy\simulations\others.py 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
class MackeyGlass(_BaseSim):
    """Simulate the Mackey-Glass delay differential system.

    TODO: Add literature values for Lyapunov etc.
    TODO: Hint the differences between this class and the other Sim classes (delay).
    TODO: Check if the structure is really good?
    TODO: Add Proper Tests.
    TODO: Decide whether to use the simple forward-euler or RK4-style update.

    Note: As the Mackey-Glass system is a delay-differential equation, the class does
    not contain a simple iterate function.
    """

    def __init__(
        self,
        a: float = 0.2,
        b: float = 0.1,
        c: int = 10,
        tau: float = 23.0,
        dt: float = 0.1,
        solver: str | Callable = "rk4",
    ) -> None:
        """Initialize the Mackey-Glass simulation object."""
        self.a = a
        self.b = b
        self.c = c
        self.tau = tau
        self.dt = dt
        self.solver = solver

        # The number of time steps between t=0 and t=-tau:
        self.history_steps = int(self.tau / self.dt)

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of MG system."""
        return np.array([0.9])

    def flow_mg(self, x: np.ndarray, x_past: np.ndarray) -> np.ndarray:
        """Calculate the flow of the Mackey-Glass equation.

        Args:
            x: The immediate value of the system. Needs to have shape (1,).
            x_past: The delayed value of the system. Needs to have shape (1,).

        Returns:
            : The flow corresponding to x and x_past.
        """
        return np.array(
            [self.a * x_past[0] / (1 + x_past[0] ** self.c) - self.b * x[0]]
        )

    def iterate_mg(self, x: np.ndarray, x_past: np.ndarray) -> np.ndarray:
        """Calculate the next time step in the Mackey-Glass equation.

        Args:
            x: The immediate value of the system. Needs to have shape (1,).
            x_past: The delayed value of the system. Needs to have shape (1,).

        Returns:
            : The next value given the immediate and delayed values.
        """

        def flow_like(x_use: np.ndarray) -> np.ndarray:
            return self.flow_mg(x_use, x_past)

        if isinstance(self.solver, str):
            if self.solver == "rk4":
                x_next = runge_kutta_4(flow_like, self.dt, x)
            elif self.solver == "forward_euler":
                x_next = forward_euler(flow_like, self.dt, x)
            else:
                raise ValueError(f"Unsupported solver: {self.solver}")
        else:
            x_next = self.solver(flow_like, self.dt, x)

        return x_next

    def simulate(
        self,
        time_steps: int,
        starting_point: np.ndarray | None = None,
        transient: int = 0,
    ) -> np.ndarray:
        """Simulate the Mackey-Glass trajectory.

        Args:
            time_steps: Number of time steps t to simulate.
            starting_point: Starting point of the trajectory shape (sys_dim,).
                            If None, take the default starting point.
            transient: Washout before storing the trajectory.

        Returns:
            Trajectory of shape (t, sys_dim).
        """
        if starting_point is None:
            starting_point = self.get_default_starting_pnt()

        if starting_point.size == self.history_steps + 1:
            initial_history = starting_point
        elif starting_point.size == 1:
            initial_history = np.repeat(starting_point, self.history_steps + 1)
        else:
            raise ValueError("Wrong size of starting point.")

        traj_w_hist = np.zeros((self.history_steps + time_steps, 1))
        traj_w_hist[: self.history_steps + 1, :] = initial_history[:, np.newaxis]

        for t in range(1, time_steps + transient):
            t_shifted = t + self.history_steps
            traj_w_hist[t_shifted] = self.iterate_mg(
                traj_w_hist[t_shifted - 1],
                traj_w_hist[t_shifted - 1 - self.history_steps],
            )

        return traj_w_hist[self.history_steps + transient :, :]

__init__(a=0.2, b=0.1, c=10, tau=23.0, dt=0.1, solver='rk4')

Initialize the Mackey-Glass simulation object.

Source code in src\lorenzpy\simulations\others.py 
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
def __init__(
    self,
    a: float = 0.2,
    b: float = 0.1,
    c: int = 10,
    tau: float = 23.0,
    dt: float = 0.1,
    solver: str | Callable = "rk4",
) -> None:
    """Initialize the Mackey-Glass simulation object."""
    self.a = a
    self.b = b
    self.c = c
    self.tau = tau
    self.dt = dt
    self.solver = solver

    # The number of time steps between t=0 and t=-tau:
    self.history_steps = int(self.tau / self.dt)

flow_mg(x, x_past)

Calculate the flow of the Mackey-Glass equation.

Parameters:

Name Type Description Default
x ndarray

The immediate value of the system. Needs to have shape (1,).

 required
x_past ndarray

The delayed value of the system. Needs to have shape (1,).

 required

Returns:

Type Description
ndarray

The flow corresponding to x and x_past.
Source code in src\lorenzpy\simulations\others.py 
168
169
170
171
172
173
174
175
176
177
178
179
180
def flow_mg(self, x: np.ndarray, x_past: np.ndarray) -> np.ndarray:
    """Calculate the flow of the Mackey-Glass equation.

    Args:
        x: The immediate value of the system. Needs to have shape (1,).
        x_past: The delayed value of the system. Needs to have shape (1,).

    Returns:
        : The flow corresponding to x and x_past.
    """
    return np.array(
        [self.a * x_past[0] / (1 + x_past[0] ** self.c) - self.b * x[0]]
    )

get_default_starting_pnt()

Return default starting point of MG system.

Source code in src\lorenzpy\simulations\others.py 
164
165
166
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of MG system."""
    return np.array([0.9])

iterate_mg(x, x_past)

Calculate the next time step in the Mackey-Glass equation.

Parameters:

Name Type Description Default
x ndarray

The immediate value of the system. Needs to have shape (1,).

 required
x_past ndarray

The delayed value of the system. Needs to have shape (1,).

 required

Returns:

Type Description
ndarray

The next value given the immediate and delayed values.
Source code in src\lorenzpy\simulations\others.py 
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
def iterate_mg(self, x: np.ndarray, x_past: np.ndarray) -> np.ndarray:
    """Calculate the next time step in the Mackey-Glass equation.

    Args:
        x: The immediate value of the system. Needs to have shape (1,).
        x_past: The delayed value of the system. Needs to have shape (1,).

    Returns:
        : The next value given the immediate and delayed values.
    """

    def flow_like(x_use: np.ndarray) -> np.ndarray:
        return self.flow_mg(x_use, x_past)

    if isinstance(self.solver, str):
        if self.solver == "rk4":
            x_next = runge_kutta_4(flow_like, self.dt, x)
        elif self.solver == "forward_euler":
            x_next = forward_euler(flow_like, self.dt, x)
        else:
            raise ValueError(f"Unsupported solver: {self.solver}")
    else:
        x_next = self.solver(flow_like, self.dt, x)

    return x_next

simulate(time_steps, starting_point=None, transient=0)

Simulate the Mackey-Glass trajectory.

Parameters:

Name Type Description Default
time_steps int

Number of time steps t to simulate.

 required
starting_point ndarray | None

Starting point of the trajectory shape (sys_dim,). If None, take the default starting point.

 None
transient int

Washout before storing the trajectory.

 0

Returns:

Type Description
ndarray

Trajectory of shape (t, sys_dim).
Source code in src\lorenzpy\simulations\others.py 
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
def simulate(
    self,
    time_steps: int,
    starting_point: np.ndarray | None = None,
    transient: int = 0,
) -> np.ndarray:
    """Simulate the Mackey-Glass trajectory.

    Args:
        time_steps: Number of time steps t to simulate.
        starting_point: Starting point of the trajectory shape (sys_dim,).
                        If None, take the default starting point.
        transient: Washout before storing the trajectory.

    Returns:
        Trajectory of shape (t, sys_dim).
    """
    if starting_point is None:
        starting_point = self.get_default_starting_pnt()

    if starting_point.size == self.history_steps + 1:
        initial_history = starting_point
    elif starting_point.size == 1:
        initial_history = np.repeat(starting_point, self.history_steps + 1)
    else:
        raise ValueError("Wrong size of starting point.")

    traj_w_hist = np.zeros((self.history_steps + time_steps, 1))
    traj_w_hist[: self.history_steps + 1, :] = initial_history[:, np.newaxis]

    for t in range(1, time_steps + transient):
        t_shifted = t + self.history_steps
        traj_w_hist[t_shifted] = self.iterate_mg(
            traj_w_hist[t_shifted - 1],
            traj_w_hist[t_shifted - 1 - self.history_steps],
        )

    return traj_w_hist[self.history_steps + transient :, :]

Roessler

Bases: _BaseSimFlow

Simulation class for the Roessler system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
class Roessler(_BaseSimFlow):
    """Simulation class for the Roessler system."""

    def __init__(
        self,
        a: float = 0.2,
        b: float = 0.2,
        c: float = 5.7,
        dt: float = 0.1,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the Roessler simulation object.

        Args:
            a: a parameter of Roessler equation.
            b: b parameter of Roessler equation.
            c: c parameter of Roessler equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.a = a
        self.b = b
        self.c = c

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of Roessler equation."""
        return np.array(
            [-x[1] - x[2], x[0] + self.a * x[1], self.b + x[2] * (x[0] - self.c)]
        )

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of Roessler system."""
        return np.array([-9.0, 0.0, 0.0])

__init__(a=0.2, b=0.2, c=5.7, dt=0.1, solver='rk4')

Initialize the Roessler simulation object.

Parameters:

Name Type Description Default
a float

a parameter of Roessler equation.

 0.2
b float

b parameter of Roessler equation.

 0.2
c float

c parameter of Roessler equation.

 5.7
dt float

Time step to simulate.

 0.1
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
def __init__(
    self,
    a: float = 0.2,
    b: float = 0.2,
    c: float = 5.7,
    dt: float = 0.1,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the Roessler simulation object.

    Args:
        a: a parameter of Roessler equation.
        b: b parameter of Roessler equation.
        c: c parameter of Roessler equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.a = a
    self.b = b
    self.c = c

flow(x)

Return the flow of Roessler equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
81
82
83
84
85
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of Roessler equation."""
    return np.array(
        [-x[1] - x[2], x[0] + self.a * x[1], self.b + x[2] * (x[0] - self.c)]
    )

get_default_starting_pnt()

Return default starting point of Roessler system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
87
88
89
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of Roessler system."""
    return np.array([-9.0, 0.0, 0.0])

Rucklidge

Bases: _BaseSimFlow

Simulation class for the Rucklidge system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
class Rucklidge(_BaseSimFlow):
    """Simulation class for the Rucklidge system."""

    def __init__(
        self,
        kappa: float = 2.0,
        lam: float = 6.7,
        dt: float = 0.1,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the Rucklidge simulation object.

        Args:
            kappa: kappa parameter of Rucklidge equation.
            lam: lambda parameter of Rucklidge equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.kappa = kappa
        self.lam = lam

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of Rucklidge equation."""
        return np.array(
            [
                -self.kappa * x[0] + self.lam * x[1] - x[1] * x[2],
                x[0],
                -x[2] + x[1] ** 2,
            ]
        )

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of Rucklidge system."""
        return np.array([1.0, 0.0, 4.5])

__init__(kappa=2.0, lam=6.7, dt=0.1, solver='rk4')

Initialize the Rucklidge simulation object.

Parameters:

Name Type Description Default
kappa float

kappa parameter of Rucklidge equation.

 2.0
lam float

lambda parameter of Rucklidge equation.

 6.7
dt float

Time step to simulate.

 0.1
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
def __init__(
    self,
    kappa: float = 2.0,
    lam: float = 6.7,
    dt: float = 0.1,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the Rucklidge simulation object.

    Args:
        kappa: kappa parameter of Rucklidge equation.
        lam: lambda parameter of Rucklidge equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.kappa = kappa
    self.lam = lam

flow(x)

Return the flow of Rucklidge equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
298
299
300
301
302
303
304
305
306
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of Rucklidge equation."""
    return np.array(
        [
            -self.kappa * x[0] + self.lam * x[1] - x[1] * x[2],
            x[0],
            -x[2] + x[1] ** 2,
        ]
    )

get_default_starting_pnt()

Return default starting point of Rucklidge system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
308
309
310
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of Rucklidge system."""
    return np.array([1.0, 0.0, 4.5])

SimplestDrivenChaotic

Bases: _BaseSimFlowDriven

Simulate the Simplest Driven Chaotic system from Sprott.

Taken from (Sprott, Julien Clinton, and Julien C. Sprott. Chaos and time-series analysis. Vol. 69. Oxford: Oxford university press, 2003.)

Source code in src\lorenzpy\simulations\driven_systems.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
class SimplestDrivenChaotic(_BaseSimFlowDriven):
    """Simulate the Simplest Driven Chaotic system from Sprott.

    Taken from (Sprott, Julien Clinton, and Julien C. Sprott. Chaos and time-series
    analysis. Vol. 69. Oxford: Oxford university press, 2003.)
    """

    def __init__(self, omega: float = 1.88, dt: float = 0.1, solver="rk4") -> None:
        """Initialize the SimplestDrivenChaotic simulation object."""
        super().__init__(dt, solver)
        self.omega = omega

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow ."""
        return np.array([x[1], -(x[0] ** 3) + np.sin(self.omega * x[2]), 1])

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point."""
        return np.array([0.0, 0.0, 0.0])

__init__(omega=1.88, dt=0.1, solver='rk4')

Initialize the SimplestDrivenChaotic simulation object.

Source code in src\lorenzpy\simulations\driven_systems.py 
16
17
18
19
def __init__(self, omega: float = 1.88, dt: float = 0.1, solver="rk4") -> None:
    """Initialize the SimplestDrivenChaotic simulation object."""
    super().__init__(dt, solver)
    self.omega = omega

flow(x)

Return the flow .

Source code in src\lorenzpy\simulations\driven_systems.py 
21
22
23
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow ."""
    return np.array([x[1], -(x[0] ** 3) + np.sin(self.omega * x[2]), 1])

get_default_starting_pnt()

Return default starting point.

Source code in src\lorenzpy\simulations\driven_systems.py 
25
26
27
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point."""
    return np.array([0.0, 0.0, 0.0])

Thomas

Bases: _BaseSimFlow

Simulation class for the Thomas system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
class Thomas(_BaseSimFlow):
    """Simulation class for the Thomas system."""

    def __init__(
        self,
        b: float = 0.18,
        dt: float = 0.3,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the Thomas simulation object.

        Args:
            b: b parameter of Thomas equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.b = b

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of Thomas equation."""
        return np.array(
            [
                -self.b * x[0] + np.sin(x[1]),
                -self.b * x[1] + np.sin(x[2]),
                -self.b * x[2] + np.sin(x[0]),
            ]
        )

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of Thomas system."""
        return np.array([0.1, 0.0, 0.0])

__init__(b=0.18, dt=0.3, solver='rk4')

Initialize the Thomas simulation object.

Parameters:

Name Type Description Default
b float

b parameter of Thomas equation.

 0.18
dt float

Time step to simulate.

 0.3
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
def __init__(
    self,
    b: float = 0.18,
    dt: float = 0.3,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the Thomas simulation object.

    Args:
        b: b parameter of Thomas equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.b = b

flow(x)

Return the flow of Thomas equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
230
231
232
233
234
235
236
237
238
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of Thomas equation."""
    return np.array(
        [
            -self.b * x[0] + np.sin(x[1]),
            -self.b * x[1] + np.sin(x[2]),
            -self.b * x[2] + np.sin(x[0]),
        ]
    )

get_default_starting_pnt()

Return default starting point of Thomas system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
240
241
242
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of Thomas system."""
    return np.array([0.1, 0.0, 0.0])

WindmiAttractor

Bases: _BaseSimFlow

Simulation class for the WindmiAttractor system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
class WindmiAttractor(_BaseSimFlow):
    """Simulation class for the WindmiAttractor system."""

    def __init__(
        self,
        a: float = 0.7,
        b: float = 2.5,
        dt: float = 0.2,
        solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
    ):
        """Initialize the WindmiAttractor simulation object.

        Args:
            a: a parameter of WindmiAttractor equation.
            b: b parameter of WindmiAttractor equation.
            dt: Time step to simulate.
            solver: The solver.
        """
        super().__init__(dt, solver)
        self.a = a
        self.b = b

    def flow(self, x: np.ndarray) -> np.ndarray:
        """Return the flow of WindmiAttractor equation."""
        return np.array([x[1], x[2], -self.a * x[2] - x[1] + self.b - np.exp(x[0])])

    def get_default_starting_pnt(self) -> np.ndarray:
        """Return default starting point of WindmiAttractor system."""
        return np.array([0.0, 0.8, 0.0])

__init__(a=0.7, b=2.5, dt=0.2, solver='rk4')

Initialize the WindmiAttractor simulation object.

Parameters:

Name Type Description Default
a float

a parameter of WindmiAttractor equation.

 0.7
b float

b parameter of WindmiAttractor equation.

 2.5
dt float

Time step to simulate.

 0.2
solver str | str | Callable[[Callable, float, ndarray], ndarray]

The solver.

 'rk4'
Source code in src\lorenzpy\simulations\autonomous_flows.py 
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
def __init__(
    self,
    a: float = 0.7,
    b: float = 2.5,
    dt: float = 0.2,
    solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
):
    """Initialize the WindmiAttractor simulation object.

    Args:
        a: a parameter of WindmiAttractor equation.
        b: b parameter of WindmiAttractor equation.
        dt: Time step to simulate.
        solver: The solver.
    """
    super().__init__(dt, solver)
    self.a = a
    self.b = b

flow(x)

Return the flow of WindmiAttractor equation.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
267
268
269
def flow(self, x: np.ndarray) -> np.ndarray:
    """Return the flow of WindmiAttractor equation."""
    return np.array([x[1], x[2], -self.a * x[2] - x[1] + self.b - np.exp(x[0])])

get_default_starting_pnt()

Return default starting point of WindmiAttractor system.

Source code in src\lorenzpy\simulations\autonomous_flows.py 
271
272
273
def get_default_starting_pnt(self) -> np.ndarray:
    """Return default starting point of WindmiAttractor system."""
    return np.array([0.0, 0.8, 0.0])

The solvers used to solve the flow equation.

create_scipy_ivp_solver(method='RK45', **additional_solve_ivp_args)

Create a scipy solver for initializing flow systems.

This function creates a scipy solver that can be used to initialize flow simulation classes. It wraps the scipy.integrate.solve_ivp function.

The scipy solvers often internally integrate more than 1 time step in the

range 0 to dt.

Parameters:

Name Type Description Default
method str

The integration method to use, e.g., 'RK45', 'RK23', 'DOP853', 'Radau', 'BDF', or 'LSODA'. See the documentation for scipy.integrate.solve_ivp for more information.

 'RK45'
**additional_solve_ivp_args

Additional arguments passed to scipy.integrate.solve_ivp as solve_ivp(fun, t_span, y0, method=method, **additional_solve_ivp_args). For example you can pass the relative and absolute tolerances as rtol=x and atol=x.

 {}

Returns:

Type Description
Callable[[Callable, float, ndarray], ndarray]

Callable[[Callable[[np.ndarray], np.ndarray], float, np.ndarray], np.ndarray]:
Callable[[Callable, float, ndarray], ndarray]

A solver function that takes three arguments:
Callable[[Callable, float, ndarray], ndarray]
  1. f (Callable[[np.ndarray], np.ndarray]): The flow function.
Callable[[Callable, float, ndarray], ndarray]
  1. dt (float): The time step for the integration.
Callable[[Callable, float, ndarray], ndarray]
  1. x (np.ndarray): The initial state.
Callable[[Callable, float, ndarray], ndarray]

The solver returns the integrated state at the end of the time step.
Source code in src\lorenzpy\simulations\solvers.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
def create_scipy_ivp_solver(
    method: str = "RK45", **additional_solve_ivp_args
) -> Callable[[Callable, float, np.ndarray], np.ndarray]:
    """Create a scipy solver for initializing flow systems.

    This function creates a scipy solver that can be used to initialize flow simulation
    classes. It wraps the scipy.integrate.solve_ivp function.

    Note: The scipy solvers often internally integrate more than 1 time step in the
        range 0 to dt.

    Args:
        method (str): The integration method to use, e.g., 'RK45', 'RK23', 'DOP853',
            'Radau', 'BDF', or 'LSODA'. See the documentation for
            scipy.integrate.solve_ivp for more information.
        **additional_solve_ivp_args: Additional arguments passed to
            scipy.integrate.solve_ivp as `solve_ivp(fun, t_span, y0, method=method,
            **additional_solve_ivp_args)`. For example you can pass the relative and
            absolute tolerances as rtol=x and atol=x.

    Returns:
        Callable[[Callable[[np.ndarray], np.ndarray], float, np.ndarray], np.ndarray]:
        A solver function that takes three arguments:
        1. f (Callable[[np.ndarray], np.ndarray]): The flow function.
        2. dt (float): The time step for the integration.
        3. x (np.ndarray): The initial state.

        The solver returns the integrated state at the end of the time step.
    """

    def solver(
        f: Callable[[np.ndarray], np.ndarray], dt: float, x: np.ndarray
    ) -> np.ndarray:
        """Solves the flow function for a time step using scipy.integrate.solve_ivp.

        Args:
            f (Callable[[np.ndarray], np.ndarray]): The flow function.
            dt (float): The time step for the integration.
            x (np.ndarray): The initial state.

        Returns:
            np.ndarray: The integrated state at the end of the time step.
        """

        def scipy_func_form(t, y):
            """Ignores the time argument for the flow f."""
            return f(y)

        out = solve_ivp(
            scipy_func_form, (0, dt), x, method=method, **additional_solve_ivp_args
        )
        return out.y[:, -1]

    return solver

forward_euler(f, dt, x)

Simulate one step for ODEs of the form dx/dt = f(x(t)) using the forward euler.

Parameters:

Name Type Description Default
f Callable[[ndarray], ndarray]

function used to calculate the time derivative at point x.

 required
dt float

time step size.

 required
x ndarray

d-dim position at time t.

 required

Returns:

Type Description
ndarray

d-dim position at time t+dt.
Source code in src\lorenzpy\simulations\solvers.py 
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def forward_euler(
    f: Callable[[np.ndarray], np.ndarray], dt: float, x: np.ndarray
) -> np.ndarray:
    """Simulate one step for ODEs of the form dx/dt = f(x(t)) using the forward euler.

    Args:
        f: function used to calculate the time derivative at point x.
        dt: time step size.
        x: d-dim position at time t.

    Returns:
       d-dim position at time t+dt.

    """
    return x + dt * f(x)

runge_kutta_4(f, dt, x)

Simulate one step for ODEs of the form dx/dt = f(x(t)) using Runge-Kutta.

Parameters:

Name Type Description Default
f Callable[[ndarray], ndarray]

function used to calculate the time derivative at point x.

 required
dt float

time step size.

 required
x ndarray

d-dim position at time t.

 required

Returns:

Type Description
ndarray

d-dim position at time t+dt.
Source code in src\lorenzpy\simulations\solvers.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def runge_kutta_4(
    f: Callable[[np.ndarray], np.ndarray], dt: float, x: np.ndarray
) -> np.ndarray:
    """Simulate one step for ODEs of the form dx/dt = f(x(t)) using Runge-Kutta.

    Args:
        f: function used to calculate the time derivative at point x.
        dt: time step size.
        x: d-dim position at time t.

    Returns:
       d-dim position at time t+dt.

    """
    k1: np.ndarray = dt * f(x)
    k2: np.ndarray = dt * f(x + k1 / 2)
    k3: np.ndarray = dt * f(x + k2 / 2)
    k4: np.ndarray = dt * f(x + k3)
    next_step: np.ndarray = x + 1.0 / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
    return next_step

timestep_iterator(f, time_steps, starting_point)

Iterate an iterator-function f: x(i+1) = f(x(i)) multiple times.

Source code in src\lorenzpy\simulations\solvers.py 
105
106
107
108
109
110
111
112
113
114
115
def timestep_iterator(
    f: Callable[[np.ndarray], np.ndarray], time_steps: int, starting_point: np.ndarray
) -> np.ndarray:
    """Iterate an iterator-function f: x(i+1) = f(x(i)) multiple times."""
    starting_point = np.array(starting_point)
    traj_size = (time_steps, starting_point.shape[0])
    traj = np.zeros(traj_size)
    traj[0, :] = starting_point
    for t in range(1, traj_size[0]):
        traj[t] = f(traj[t - 1])
    return traj

measures module:

Measures for chaotic dynamical systems.

largest_lyapunov_exponent(iterator_func, starting_point, deviation_scale=1e-10, steps=int(1000.0), part_time_steps=15, steps_skip=50, dt=1.0, initial_pert_direction=None, return_convergence=False)

Numerically calculate the largest lyapunov exponent given an iterator function.

See: Sprott, Julien Clinton, and Julien C. Sprott. Chaos and time-series analysis. Vol. 69. Oxford: Oxford university press, 2003.

Parameters:

Name Type Description Default
iterator_func Callable[[ndarray], ndarray]

Function to iterate the system to the next time step: x(i+1) = F(x(i))

 required
starting_point ndarray

The starting_point of the main trajectory.

 required
deviation_scale float

The L2-norm of the initial perturbation.

 1e-10
steps int

Number of renormalization steps.

 int(1000.0)
part_time_steps int

Time steps between renormalization steps.

 15
steps_skip int

Number of renormalization steps to perform, before tracking the log divergence. Avoid transients by using steps_skip.

 50
dt float

Size of time step.

 1.0
initial_pert_direction ndarray | None
  • If np.ndarray: The direction of the initial perturbation.
  • If None: The direction of the initial perturbation is assumed to be np.ones(..).
 None
return_convergence bool

If True, return the convergence of the largest LE; a numpy array of the shape (N, ).

 False

Returns:

Type Description
float | ndarray

The largest Lyapunov Exponent. If return_convergence is True: The convergence
float | ndarray

(np.ndarray), else just the float value, which is the last value in the
float | ndarray

convergence.
Source code in src\lorenzpy\measures.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
def largest_lyapunov_exponent(
    iterator_func: Callable[[np.ndarray], np.ndarray],
    starting_point: np.ndarray,
    deviation_scale: float = 1e-10,
    steps: int = int(1e3),
    part_time_steps: int = 15,
    steps_skip: int = 50,
    dt: float = 1.0,
    initial_pert_direction: np.ndarray | None = None,
    return_convergence: bool = False,
) -> float | np.ndarray:
    """Numerically calculate the largest lyapunov exponent given an iterator function.

    See: Sprott, Julien Clinton, and Julien C. Sprott. Chaos and time-series analysis.
    Vol. 69. Oxford: Oxford university press, 2003.

    Args:
        iterator_func: Function to iterate the system to the next time step:
                       x(i+1) = F(x(i))
        starting_point: The starting_point of the main trajectory.
        deviation_scale: The L2-norm of the initial perturbation.
        steps: Number of renormalization steps.
        part_time_steps: Time steps between renormalization steps.
        steps_skip: Number of renormalization steps to perform, before tracking the log
                    divergence. Avoid transients by using steps_skip.
        dt: Size of time step.
        initial_pert_direction:
            - If np.ndarray: The direction of the initial perturbation.
            - If None: The direction of the initial perturbation is assumed to be
                np.ones(..).
        return_convergence: If True, return the convergence of the largest LE; a numpy
                            array of the shape (N, ).

    Returns:
        The largest Lyapunov Exponent. If return_convergence is True: The convergence
        (np.ndarray), else just the float value, which is the last value in the
        convergence.
    """
    x_dim = starting_point.size

    if initial_pert_direction is None:
        initial_pert_direction = np.ones(x_dim)

    initial_perturbation = initial_pert_direction * (
        deviation_scale / np.linalg.norm(initial_pert_direction)
    )

    log_divergence = np.zeros(steps)

    x = starting_point
    x_pert = starting_point + initial_perturbation

    for i_n in range(steps + steps_skip):
        for i_t in range(part_time_steps):
            x = iterator_func(x)
            x_pert = iterator_func(x_pert)
        dx = x_pert - x
        norm_dx = np.linalg.norm(dx)
        x_pert = x + dx * (deviation_scale / norm_dx)
        if i_n >= steps_skip:
            log_divergence[i_n - steps_skip] = np.log(norm_dx / deviation_scale)

    if return_convergence:
        return np.array(
            np.cumsum(log_divergence) / (np.arange(1, steps + 1) * dt * part_time_steps)
        )
    else:
        return float(np.average(log_divergence) / (dt * part_time_steps))

lyapunov_exponent_spectrum(iterator_func, starting_point, deviation_scale=1e-10, steps=int(1000.0), part_time_steps=15, steps_skip=50, dt=1.0, m=None, initial_pert_directions=None, return_convergence=False)

Calculate the spectrum of m largest lyapunov exponent given an iterator function.

A mixture of: - The algorithm for the largest lyapunov exponent: Sprott, Julien Clinton, and Julien C. Sprott. Chaos and time-series analysis. Vol. 69. Oxford: Oxford university press, 2003. - The algorithm for the spectrum given in 1902.09651 "LYAPUNOV EXPONENTS of the KURAMOTO-SIVASHINSKY PDE".

Parameters:

Name Type Description Default
iterator_func Callable[[ndarray], ndarray]

Function to iterate the system to the next time step: x(i+1) = F(x(i))

 required
starting_point ndarray

The starting_point of the main trajectory.

 required
deviation_scale float

The L2-norm of the initial perturbation.

 1e-10
steps int

Number of renormalization steps.

 int(1000.0)
part_time_steps int

Time steps between renormalization steps.

 15
steps_skip int

Number of renormalization steps to perform, before tracking the log divergence. Avoid transients by using steps_skip.

 50
dt float

Size of time step.

 1.0
m int | None

Number of Lyapunov exponents to compute. If None: take all (m = x_dim).

 None
initial_pert_directions ndarray | None
  • If np.ndarray: The direction of the initial perturbations of shape (m, x_dim). Will be qr decomposed and the q matrix will be used.
  • If None: The direction of the initial perturbation is assumed to be np.eye(x_dim)[:m, :].
 None
return_convergence bool

If True, return the convergence of the largest LE; a numpy array of the shape (N, m).

 False

Returns:

Name Type Description
ndarray

The Lyapunov exponent spectrum of largest m values. If return_convergence is
True ndarray

The convergence (2D N x m np.ndarray), else a (1D m-size np.ndarray),
ndarray

which holds the last values in the convergence.
Source code in src\lorenzpy\measures.py 
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
def lyapunov_exponent_spectrum(
    iterator_func: Callable[[np.ndarray], np.ndarray],
    starting_point: np.ndarray,
    deviation_scale: float = 1e-10,
    steps: int = int(1e3),
    part_time_steps: int = 15,
    steps_skip: int = 50,
    dt: float = 1.0,
    m: int | None = None,
    initial_pert_directions: np.ndarray | None = None,
    return_convergence: bool = False,
) -> np.ndarray:
    """Calculate the spectrum of m largest lyapunov exponent given an iterator function.

    A mixture of:
    - The algorithm for the largest lyapunov exponent: Sprott, Julien Clinton, and
        Julien C. Sprott. Chaos and time-series analysis. Vol. 69. Oxford: Oxford
        university press, 2003.
    - The algorithm for the spectrum given in 1902.09651 "LYAPUNOV EXPONENTS of the
        KURAMOTO-SIVASHINSKY PDE".

    Args:
        iterator_func: Function to iterate the system to the next time step:
                        x(i+1) = F(x(i))
        starting_point: The starting_point of the main trajectory.
        deviation_scale: The L2-norm of the initial perturbation.
        steps: Number of renormalization steps.
        part_time_steps: Time steps between renormalization steps.
        steps_skip: Number of renormalization steps to perform, before tracking the log
                    divergence.
                Avoid transients by using steps_skip.
        dt: Size of time step.
        m: Number of Lyapunov exponents to compute. If None: take all (m = x_dim).
        initial_pert_directions:
            - If np.ndarray: The direction of the initial perturbations of shape
                            (m, x_dim). Will be qr decomposed
                             and the q matrix will be used.
            - If None: The direction of the initial perturbation is assumed to be
                        np.eye(x_dim)[:m, :].
        return_convergence: If True, return the convergence of the largest LE; a
                            numpy array of the shape (N, m).

    Returns:
        The Lyapunov exponent spectrum of largest m values. If return_convergence is
        True: The convergence (2D N x m np.ndarray), else a (1D m-size np.ndarray),
        which holds the last values in the convergence.
    """
    x_dim = starting_point.size

    if m is None:
        m = x_dim

    if initial_pert_directions is None:
        initial_perturbations = (
            np.eye(x_dim)[:m, :] * deviation_scale
        )  # each vector is of length deviation_scale
    else:
        q, _ = np.linalg.qr(initial_pert_directions)
        initial_perturbations = q * deviation_scale

    log_divergence_spec = np.zeros((steps, m))

    x = starting_point
    x_pert = starting_point + initial_perturbations

    for i_n in range(steps + steps_skip):
        for i_t in range(part_time_steps):
            x = iterator_func(x)
            for i_m in range(m):
                x_pert[i_m, :] = iterator_func(x_pert[i_m, :])
        dx = x_pert - x
        q, r = np.linalg.qr(dx.T, mode="reduced")
        x_pert = x + q.T * deviation_scale

        if i_n >= steps_skip:
            log_divergence_spec[i_n - steps_skip, :] = np.log(
                np.abs(np.diag(r)) / deviation_scale
            )

    if return_convergence:
        return np.array(
            np.cumsum(log_divergence_spec, axis=0)
            / (
                np.repeat(np.arange(1, steps + 1)[:, np.newaxis], m, axis=1)
                * dt
                * part_time_steps
            )
        )
    else:
        return np.average(log_divergence_spec, axis=0) / (dt * part_time_steps)