The temporal evolution of recurrent chaotic systems is marked by self-similar patterns that reoccur in finite time. Due to the sensitivity of such systems to both initial conditions and the precision of the numerical scheme, simply integrating their governing equations over long periods of time adds little additional information. Under this circumstance, such a systems temporal evolution can be extrapolated from a given point on by simply replaying sequences from the systems history.

In this post, I will apply a procedure adapted from literature to the Lorentz systems. First, a dataset is generated from its governing equations, then a recurrence plot generated and finally extrapolated.

The Lorentz system

The famous Lorentz system is given by: \begin{align} \frac{dX}{dt} = \begin{bmatrix} \sigma (X_2-X_1) \\ X_1 (\rho-X_3) - X_2 \\ X_1*X_2 - \beta X_3 \end{bmatrix} \end{align} Where $\sigma=10$, $\beta=8/3$ and $\rho=28$ are parameter choices that give recurrent behavior.

In Python, this can be implemented in a few lines:

def lorentz_attractor(t, W, sigma, rho, beta):
    x, y, z = W
    return [
        sigma*(y-x),
        x*(rho-z)-y,
        x*y - beta*z,
    ]

import numpy as np
import pandas as pd
from scipy.integrate import solve_ivp

t_span = [0, 50]
integration_result = solve_ivp(
    lambda t, W: lorentz_attractor(t, W, sigma=10, beta=8 / 3, rho=28),
    t_span=t_span,
    y0=[1, 0, 0],
    t_eval=np.linspace(*t_span, 5000),
)
trajectory = pd.DataFrame(
    np.transpose([integration_result.t, *integration_result.y]),
    columns="t x y z".split(),
)

This trajectory can be plotted like this:

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(ncols=2, sharey=True)
trajectory.plot('t', ['y', 'z', 'x'], ax=ax1, label='$X_2$ $X_3$ $X_1$'.split())
trajectory.plot('x', ['y', 'z'], ax=ax2)
ax2.set(aspect='equal', xlabel='$X_1$')
ax2.get_legend().remove()

The dataframe trajectory now contains a sample of the state space of the Lorentz system. After a spin-up period that runs up to $t=18$, the $X_1$ and $X_2$ coordinates oscillate around either $+9$ or $-9$, with attraction point switch interval of $1–5$, and $X_3$ oscillating around $28$.

Recurrence plot generation

The recurrence plot $\mathcal{R}(t_i, t_j)$ quantifies the difference between two states of the system. Here, we use the $\mathcal{l}^2$-norm to compare the states (= coordinates $X$) of the system at different times: \begin{align} \mathcal{R}(t_i, t_j) = \frac{1}{3}\sqrt{\sum^3_{k=1}{\left(X_k(t_i)-X_k(t_j)\right)^2}} \end{align} The result of this is a symmetric matrix, meaning that only the half of the entries need to be calculated:

from numpy.linalg import norm # l2 - norm

def recurrence_matrix(trajectory, normed=False):
    rec_matrix = np.zeros((len(trajectory), len(trajectory)))
    np_traj_values = trajectory[["x", "y", "z"]].to_numpy()
    for i in range(1, rec_matrix.shape[0]):
        values = norm(
            np_traj_values[i, :] - np_traj_values[(i + 1) :, :],
            axis=-1,
        )
        rec_matrix[i, (i + 1) :] = values
        rec_matrix[(i + 1) :, i] = values
    if normed:
        rec_matrix /= rec_matrix.max()
    return trajectory["t"], rec_matrix

_, rec_matrix = recurrence_matrix(trajectory, normed=True)

t_range = (trajectory.t.min(), trajectory.t.max())
plt.imshow(
    rec_matrix,
    vmin=0, vmax=1,
    aspect="equal",
    origin="lower",
    extent=(*t_range, *t_range),
)
plt.colorbar(label='$\\mathcal{R}(t_i,t_j)$')
plt.xlabel('$t_i$')
plt.ylabel('$t_j$')