Quickstart

This guide walks you through the basics of using pynamicalsys.

Discrete-time dynamical system

Creating a discrete-time dynamical system object

To get started, you need to create a discrete dynamical system object. This is done using the DiscreteDynamicalSystem class. For this example, we will use the logistic map, defined as:

\[\begin{equation*} x_{n+1} = r x_n (1 - x_n). \end{equation*}\]

This map is a discrete dynamical system that exhibits a wide range of behaviors depending on the parameter \(r\). It is often used as a classic example in chaos theory. To create the discrete-time dynamical system object, we need to instanciate the DiscreteDynamicalSystem class using the model parameter, since the logistic map is built-in within this class:

from pynamicalsys import DiscreteDynamicalSystem as dds  # Import the discrete-time system class
ds = dds(model="logistic map")  # Create the logistic map discrete system object

Generating a trajectory

We are going to generate a trajectory for this system using four different parameters values. Each one of these values produces a different dynamical behavior.

x0 = 0.2  # Initial condition for x
r = [2.6, 3.1, 3.5, 3.8]  # List of parameter values

# Generate trajectories for each parameter value (100 iterations each)
trajectories = [
    ds.trajectory(x0, 100, parameters=r[i])
    for i in range(len(r))
]

Visualizing the trajectory

To visualize the trajectory, we can use the PlotStyler class to customize our plots.

from pynamicalsys import PlotStyler  # For consistent plot styling
import seaborn as sns  # For color palettes
import matplotlib.pyplot as plt  # For plotting

# Apply the plot style
ps = PlotStyler()
ps.apply_style()

# Create the figure and axes
fig, ax = plt.subplots(figsize=(10, 4))

# Define colors for each trajectory
colors = sns.color_palette("hls", n_colors=len(r))

# Plot each trajectory with a different color and label
for i, traj in enumerate(trajectories):
    ax.plot(traj, "-o", color=colors[i], label=f"$r = {r[i]}$")

# Customize the plot labels and limits
plt.xlabel("$n$")  # Iteration index
plt.ylabel("$x$")  # State variable
plt.legend(
    loc="upper center",
    frameon=False,
    ncol=4,
    bbox_to_anchor=(0.5, 1.15)
)
plt.ylim(0.15, 1)
plt.xlim(-1, 100)

plt.show()  # Display the plot

_images/logistic_map_trajectories.png — Logistic map trajectories for different parameter values.

Continuous-time dynamical system

Creating a continuous-time dynamical system object

The continuous-time analysis is similar to the discrete-time analysis. To get started, you need to create a continuous-time dynamical system object. This is done using the ContinuousDynamicalSystem class. For this example, we will use the Lorenz system, defined as:

\[\begin{split}\begin{align*} \dot{x} &= \sigma(y - x),\\ \dot{y} &= x(\rho - z) - y,\\ \dot{z} &= xy - \beta z. \end{align*}\end{split}\]

For this example, we are going to use the same parameters Lorenz used in his original paper in 1963: \(\sigma = 10\), \(\sigma = 28\), and \(\beta = 8/3\). The system exhibits chaotic behavior for this set of parameters.

To create the continuous-time dynamical system object, we need to instanciate the ContinuousDynamicalSystem class using the model parameter, since the Lorenz system is built-in within this class:

from pynamicalsys import ContinuousDynamicalSystem as cds  # Import the continuous-time system class
ds = cds(model="lorenz system")  # Create the Lorenz system object

Generating a trajectory

We are going to generate a trajectory for this system using the mentioned parameters. The order in which the parameters must be given for the built-in system can be verified using the info property.

# Initial condition for (x, y, z)
u = [0.1, 0.1, 0.1]

# Parameters of the Lorenz system (σ, ρ, β)
sigma, rho, beta = 10, 28, 8/3
parameters = [sigma, rho, beta]
ds.set_parameters(parameters)

# Total integration time
total_time = 200

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

Visualizing the trajectory

The trajectory method returns the time samples and the coordinates of the system at the respective samples. If we don’t specify the integrator, it uses the 4th order Runge-Kutta method with a fixed time step of 0.01. We can then visualize the evolution of each coordiate:

from pynamicalsys import PlotStyler  # For consistent plot styling

# Apply the plot style
ps = PlotStyler(fontsize=18)
ps.apply_style()

# Create the figure and axes for x(t), y(t), z(t)
fig, ax = plt.subplots(3, 1, sharex=True, figsize=(10, 7))

# Plot x(t), y(t), z(t) separately
for i in range(3):
    ax[i].plot(trajectory[:, 0], trajectory[:, i + 1], "k")  # time vs coordinate

# Add axis labels and limits
ax[0].set_ylabel("$x(t)$")
ax[1].set_ylabel("$y(t)$")
ax[2].set_ylabel("$z(t)$")
ax[-1].set_xlabel("$t$")
ax[0].set_xlim(0, total_time)

plt.show()  # Display the plot

_images/lorenz_time_series.png — A chaotic trajectory of the Lorenz system.

We can also visualize the attractor (a projection onto the \(xz\) plane):

ps = PlotStyler(fontsize=18, linewidth=0.3)  # Style for attractor plot
ps.apply_style()

# Plot the Lorenz attractor projection on the x-z plane
plt.plot(trajectory[:, 1], trajectory[:, 3], "k-")

plt.xlabel("$x$")
plt.ylabel("$z$")

plt.show()  # Display the attractor plot

_images/lorenz_attractor.png — The Lorenz attractor.

Hamiltonian systems

Creating a Hamiltonian system object

To get started, you need to create a Hamiltonian system object. This is done using the HamiltonianSystem class. For this example, we will use the two degrees of freedom Hénon-Hailes system, defined by the Hamiltonian function:

\[\begin{align*} H(x, y, p_x, p_y) = \frac{1}{2}(p_x^2 + p_y^2) + \frac{1}{2}(x^2 + y^2) + x^2y - \frac{y^3}{3}. \end{align*}\]

with equations of motion:

\[\begin{split}\begin{align*} \dot{x} &= \frac{\partial H}{\partial p_x} = p_x,\\ \dot{y} &= \frac{\partial H}{\partial p_y} = p_y,\\ \dot{p}_x &= -\frac{\partial H}{\partial x} = x (2y - 1),\\ \dot{p}_y &= -\frac{\partial H}{\partial y} = y^2 - y - x^2.\\ \end{align*}\end{split}\]

This system is a paradigmatic example of a Hamiltonian system that exhibits both regular and chaotic solutions. The create the Hamiltonian system object, we need to instanciante the HamiltonianSystem class using the model parameter, since the Hénon-Heiles system is built in within this class:

from pynamicalsys import HamiltonianSystem  # Import the Hamiltonian system class
hs = HamiltonianSystem(model="henon heiles")  # Create the Hénon-Heiles Hamiltonian system object

Generating Poincaré section

To visualize the different behaviors of this system, we are going to generate the Poincaré section for an ensemble of randomly chosen initial conditions:

num_ic = 100  # Number of initial conditions
dof = 2  # Degrees of freedom of the system

# Allocate arrays for initial positions q and momenta p
q = np.zeros((num_ic, dof))
p = np.zeros((num_ic, dof))

E_ref = 1 / 8  # Total energy of the system
x = 0  # Fixed x = 0 for the Poincaré section

# Ranges for y and py sampling
y_range = (-0.5, 0.5)
py_range = (-0.5, 0.5)

# Randomly generate initial conditions consistent with energy
np.random.seed(13)
for i in range(num_ic):
    while True:
        py = np.random.uniform(*py_range)
        y = np.random.uniform(*y_range)
        # Energy constraint to solve for px
        px_squared = 2 * (E_ref - x**2 * y + y**3 / 3) - x**2 - y**2 - py**2
        if px_squared > 0:  # Ensure positive px_squared
            q[i] = [x, y]
            p[i] = [np.sqrt(px_squared), py]
            break

num_intersections = 10000  # Number of Poincaré section intersections to compute

We choose as our section the \(x = 0\) plane with \(\dot{x} > 0\) and integrate the system using three symplectic integrators: the second-order velocity-Verlet integrator (VV2), the fourth-order Yoshida method (SVY4), and the implicit midpoint method (IMP). The first two (VV2 and SVY4) are specific to separable Hamiltonians, i.e., \(H(\mathbf{q}, \mathbf{p}) = T(\mathbf{p}) + V(\mathbf{q})\), whereas the last one (IMP) is aplicable to general Hamiltonians.

time_step = 0.01  # Time step for the integrators

# Compute Poincaré section using the VV2 integrator
hs.integrator("vv2", time_step=time_step)
PS_vv2 = hs.poincare_section(q.copy(), p.copy(), num_intersections)

# Compute Poincaré section using the SVY4 integrator
hs.integrator("svy4", time_step=time_step)
PS_svy4 = hs.poincare_section(q.copy(), p.copy(), num_intersections)

# Compute Poincaré section using the IMP integrator
hs.integrator("imp", time_step=time_step)
PS_imp = hs.poincare_section(q.copy(), p.copy(), num_intersections)

Visualizing the Poincaré section

To visualize the trajectory, we can use the PlotStyler class to customize our plots. We plot the Poincaré section of each trajectory in different colors to highlight the different behaviors present in the system:

from pynamicalsys import PlotStyler  # For consistent plot styling
import seaborn as sns  # For color palettes
import matplotlib.pyplot as plt  # For plotting

# Apply plot style for scatter plots
fontsize = 18
ps = PlotStyler(fontsize=fontsize, markersize=0.25, markeredgewidth=0)
ps.apply_style()

# Define colors for each initial condition
colors = sns.color_palette("tab10", num_ic)

# Create side-by-side figures for VV2, SVY4, and IMP
fig, ax = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12, 4))

# Plot each Poincaré section point
for i in range(num_ic):
    ax[0].plot(PS_vv2[i, :, 2], PS_vv2[i, :, 4], "o", c=colors[i])  # VV2 plot
    ax[1].plot(PS_svy4[i, :, 2], PS_svy4[i, :, 4], "o", c=colors[i])  # SVY4 plot
    ax[2].plot(PS_imp[i, :, 2], PS_imp[i, :, 4], "o", c=colors[i])  # IMP plot

# Add axis labels and titles
ax[0].set_xlabel("$y$")
ax[0].set_ylabel("$p_y$")
ax[0].set_title("VV2", fontsize=fontsize)
ax[1].set_xlabel("$y$")
ax[1].set_title("SVY4", fontsize=fontsize)
ax[2].set_xlabel("$y$")
ax[2].set_title("IMP", fontsize=fontsize)

plt.tight_layout(pad=0.2)  # Improve layout spacing
plt.show()  # Display the figure

_images/henon_heiles_poincare_section.png — The Poincaré section for the Hénon-Heiles system using the VV2 integrator, the SVY4 integrator, and the IMP integrator. Each color corresponds to a different initial condition.

Quickstart

Discrete-time dynamical system

Creating a discrete-time dynamical system object

Generating a trajectory

Visualizing the trajectory

Continuous-time dynamical system

Creating a continuous-time dynamical system object

Generating a trajectory

Visualizing the trajectory

Hamiltonian systems

Creating a Hamiltonian system object

Generating Poincaré section

Visualizing the Poincaré section

Further reading