Lyapunov exponents

For continuous-time dynamical systems, the Lyapunov exponents are computed by integrating simultaneaously the equations of motion and the variational equations. Given a \(d\)-dimensional continuous-time dynamical system described by the following differential equation \(\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x})\), where \(\mathbf{x}\in\mathbb{R}^d\), let \(\mathbf{J}(\mathbf{x}, t) = \mathbf{Df}(\mathbf{x}, t)\) be the Jacobian matrix of the vector field \(\mathbf{f}\) evaluated at the point \(\mathbf{x}\) at the instant of time \(t\).

Let \(\mathbf{v}_0 \equiv \mathbf{v}(0)\) be a deviation vector representing a small displacement from the initial condition \(\mathbf{x}_0 \equiv \mathbf{x}(0)\). The time evolution of \(\mathbf{v}_0\) under the linearized dynamics is given by the variational equation

\[\frac{d\mathbf{v}(t)}{dt} = J(\mathbf{x}(t), t)\mathbf{v}(t).\]

Let \(A(t)\in\mathbb{R}^{d\times d}\) be a matrix whose columns are the deviation vectors. We follow the evolution of the deviation vectors along the trajectory and reorthonormalize them using a QR decomposition: \(A(t) = Q(t) R(t)\), where \(Q(t)\) is an orthogonal matrix and \(R(t)\) is an upper triangular matrix. The Lyapunov exponents are then computed from the averages of the logarithm of the diagonal elements of the matrix \(R(t)\), \(|r_{ii}(t)|\):

\[\lambda_i = \lim_{t\rightarrow\infty}\frac{1}{t}\sum_{\tau = 0}^t\log|r_{ii}(\tau)|.\]

The Lyapunov exponents can be calculated using the lyapunov method from the ContinuousDynamicalSystem class.

To illustrate the Lyapunov exponents calculation, let’s consider the Rössler system:

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

where \(a\), \(b\), and \(c\) are the parameters of the system. Let’s first visualize the Rössler attractor:

from pynamicalsys import ContinuousDynamicalSystem as cds
from pynamicalsys import PlotStyler
import matplotlib.pyplot as plt

ds = cds(model="rossler system")
ds.integrator("rk45", atol=1e-15, rtol=1e-15)

# Parameters of the system
a, b, c = 0.15, 0.20, 10
parameters = [a, b, c]
ds.set_parameters(parameters)

# Initial conditions
u = [0.1, 0.1, 0.1]

# Total and transient time
total_time = 2000
transient_time = 1000

# Calculate the trajectory
trajectory = ds.trajectory(u, total_time, parameters=parameters, transient_time=transient_time)

# Set the plot style
ps = PlotStyler(fontsize=18, linewidth=0.1)
ps.apply_style()

# Create the 3D figure and axis
fig, ax = plt.subplots(subplot_kw={'projection': '3d'}, figsize=(5, 4))

# Plot the trajectory
ax.plot(trajectory[:, 1], trajectory[:, 2], trajectory[:, 3], label='Trajectory', color="k")

# Set the labels and view angle
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_zlabel("$z$")
ax.view_init(elev=30, azim=-140)

plt.show()

The Rössler attractor for \(a = 0.15\), \(b = 0.20\), and \(c = 10\).

Now, for the Lyapunov exponents:

total_time = 10000
lyapunov_exponents = ds.lyapunov(
    u,
    total_time,
    transient_time=transient_time,
    log_base=2,
)
print(lyapunov_exponents)

[1.27122925e-01  3.05263459e-04 -1.41384990e+01]

It is also possible to return the whole history of all Lyapunov exponents:

# Using return_history=True
lyapunov_exponents = ds.lyapunov(
    u,
    total_time,
    parameters=parameters,
    transient_time=transient_time,
    return_history=True,
)

# Set the plot style
ps = PlotStyler(fontsize=18, linewidth=0.75)
ps.apply_style()

# Create the figure and axes
fig, ax = plt.subplots(1, 2, figsize=(10, 3), sharex=True)

# Plot each Lyapunov exponent with a different color
colors = ["green", "gold", "blue"]
for i in range(3):
    ax[0].plot(
        lyapunov_exponents[:, 0],
        lyapunov_exponents[:, i + 1],
        color=colors[i],
    )

    ax[1].plot(
        lyapunov_exponents[:, 0],
        lyapunov_exponents[:, i + 1],
        color=colors[i],
        label=rf"$\lambda_{i + 1}$",
    )

# Set the legend, limits, and labels
ax[1].legend(frameon=False, ncol=3)
ax[0].set_ylim(-11, 1)
ax[1].set_ylim(1e-5, 1e0)
ax[0].set_xlim(transient_time, total_time)
ax[1].set_yscale("log")

ax[0].set_ylabel("Lyapunov exponents")
ax[0].set_xlabel("Time")
ax[1].set_xlabel("Time")

plt.show()

The Lyapunov exponents history for the Rössler system.