Added examples folder and updated project readme.

- also added new attractor_animation.gif.

README.md

@@ -58,6 +58,8 @@ lle = lpy.measures.largest_lyapunov_exponent(
 The calculated largest Lyapunov exponent of *0.9051...* is very close to the literature
 value of *0.9056*[^SprottChaos].
 
+For more examples see the [examples folder](examples/README.md).
+
 ## 💫 Supported systems
 
 

docs/reference.md

@@ -3,6 +3,7 @@
 
 ### simulations module:
 ::: lorenzpy.simulations
+::: lorenzpy.simulations.solvers
 
 ### measures module:
 ::: lorenzpy.measures

examples/README.md

@@ -0,0 +1,8 @@
+# Examples:
+
+### Contents
+- ``double_pendulum_lyapunov_spectrum.py``: Plot the lyapunov spectrum of the double pendulum.
+- ``lyapunov_spectra_of_3d_autonomous_flows.py``: Plot the Lyapunov spectrum of all 3D autonomous dissipative flow systems.
+
+⚠️ Not many examples here yet, and examples might be flawed.
+

examples/double_pendulum_lyapunov_spectrum.py

@@ -0,0 +1,68 @@
+"""Lyapunov Spectrum of single system."""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from lorenzpy import measures as meas
+from lorenzpy import simulations as sims
+
+sys_obj = sims.DoublePendulum(dt=0.1)
+dt = sys_obj.dt
+
+# Calculate exponents:
+m = 4
+deviation_scale = 1e-10
+steps = 1000
+part_time_steps = 15
+steps_skip = 0
+
+iterator_func = sys_obj.iterate
+starting_point = sys_obj.get_default_starting_pnt()
+
+le_spectrum = meas.lyapunov_exponent_spectrum(
+    iterator_func=iterator_func,
+    starting_point=starting_point,
+    deviation_scale=deviation_scale,
+    steps=steps,
+    part_time_steps=part_time_steps,
+    steps_skip=steps_skip,
+    dt=dt,
+    m=m,
+    initial_pert_directions=None,
+    return_convergence=True,
+)
+
+fig, ax = plt.subplots(
+    1, 1, figsize=(6, 6), layout="constrained", sharex=True, sharey=False
+)
+
+# x and y titles:
+fig.supxlabel("Number of renormalization steps")
+fig.supylabel("Lyapunov exponent convergence")
+
+x = np.arange(1, steps + 1)
+ax.plot(
+    x,
+    le_spectrum,
+    linewidth=1,
+)
+ax.grid(True)
+
+final_les = np.round(le_spectrum[-1, :], 4).tolist()
+final_les = [str(x) for x in final_les]
+le_string = "\n".join(final_les)
+le_string = "Final LEs: \n" + le_string
+x_position = 0.1  # X-coordinate of the upper-left corner for each subplot
+y_position = 0.5
+ax.text(
+    x_position,
+    y_position,
+    le_string,
+    fontsize=10,
+    bbox=dict(facecolor="white", edgecolor="black", boxstyle="round"),
+    verticalalignment="center",
+    horizontalalignment="left",
+    transform=ax.transAxes,
+)
+
+plt.show()

examples/lyapunov_spectra_of_3d_autonomous_flows.py

@@ -0,0 +1,91 @@
+"""Calculate and plot the Lyapunov Spectra of all three-dimensional chaotic flows."""
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from lorenzpy import measures as meas
+from lorenzpy import simulations as sims
+
+systems = [
+    "Lorenz63",
+    "Chen",
+    "ChuaCircuit",
+    "ComplexButterfly",
+    "DoubleScroll",
+    "Halvorsen",
+    "Roessler",
+    "Rucklidge",
+    "Thomas",
+    "WindmiAttractor",
+]
+
+# Calculate exponents:
+m = 3
+deviation_scale = 1e-10
+steps = 1000
+part_time_steps = 15
+steps_skip = 50
+
+solver = "rk4"
+# solver = sims.solvers.create_scipy_ivp_solver(method="RK45")
+
+lyap_dict = {}
+for i_sys, system in enumerate(systems):
+    print(system)
+    sys_obj = getattr(sims, system)(solver=solver)
+    iterator_func = sys_obj.iterate
+    starting_point = sys_obj.get_default_starting_pnt()
+    dt = sys_obj.dt
+
+    lyap_dict[system] = meas.lyapunov_exponent_spectrum(
+        iterator_func=iterator_func,
+        starting_point=starting_point,
+        deviation_scale=deviation_scale,
+        steps=steps,
+        part_time_steps=part_time_steps,
+        steps_skip=steps_skip,
+        dt=dt,
+        m=m,
+        initial_pert_directions=None,
+        return_convergence=True,
+    )
+
+fig, axs = plt.subplots(
+    2, 5, figsize=(15, 8), layout="constrained", sharex=True, sharey=False
+)
+
+# x and y titles:
+fig.supxlabel("Number of renormalization steps")
+fig.supylabel("Lyapunov exponent convergence")
+
+axs = axs.flatten()
+x = np.arange(1, steps + 1)
+for i_ax, ax in enumerate(axs):
+    system = systems[i_ax]
+    le_spectrum = lyap_dict[system]
+    ax.title.set_text(system)
+    ax.plot(
+        x,
+        le_spectrum,
+        linewidth=1,
+    )
+    ax.grid(True)
+
+    final_les = np.round(le_spectrum[-1, :], 4).tolist()
+    final_les = [str(x) for x in final_les]
+    le_string = "\n".join(final_les)
+    le_string = "Final LEs: \n" + le_string
+    x_position = 0.1  # X-coordinate of the upper-left corner for each subplot
+    y_position = 0.5
+    ax.text(
+        x_position,
+        y_position,
+        le_string,
+        fontsize=10,
+        bbox=dict(facecolor="white", edgecolor="black", boxstyle="round"),
+        verticalalignment="center",
+        horizontalalignment="left",
+        transform=ax.transAxes,
+    )
+
+plt.show()

src/lorenzpy/simulations/autonomous_flows.py

@@ -95,7 +95,7 @@ class ComplexButterfly(_BaseSimFlow):
     def __init__(
         self,
         a: float = 0.55,
-        dt: float = 0.05,
+        dt: float = 0.1,
         solver: str | str | Callable[[Callable, float, np.ndarray], np.ndarray] = "rk4",
     ):
         """Initialize the ComplexButterfly simulation object.

static/gen_attractor_animation.py

@@ -21,7 +21,7 @@
 ]
 
 # create data:
-N = 1000
+N = 1500
 data_dict = {}
 for i_sys, system in enumerate(systems):
     sys_class = sim_class = getattr(sims, system)
@@ -37,7 +37,7 @@
     system = systems[i_ax]
     data = data_dict[system]
     ax.title.set_text(system)
-    ax.plot(data[:, 0], data[:, 1], data[:, 2], linewidth=0.2, color="white")
+    ax.plot(data[:, 0], data[:, 1], data[:, 2], linewidth=0.12, color="white")
     ax.set_xticks([])  # Remove x-axis ticks
     ax.set_yticks([])  # Remove y-axis ticks
     ax.set_zticks([])  # Remove y-axis ticks
@@ -50,9 +50,9 @@
     ax.axis("off")
 
 
-def update_graph(num):
+def update_graph(angle):
     """Rotate camera angle."""
-    azim = num
+    azim = angle
     roll = 0
     elev = 0
     for i_ax, ax in enumerate(axs):
@@ -67,9 +67,12 @@ def update_graph(num):
 for prev_file in previous_files:
     os.remove(prev_file)
 
-for num in range(0, 360):
-    update_graph(num)
-    plt.savefig(f"{frames_dir}/frame_{str(num).zfill(3)}.png", transparent=True, dpi=50)
+for num in range(0, 90):
+    angle = 4 * num
+    update_graph(angle)
+    plt.savefig(
+        f"{frames_dir}/frame_{str(num).zfill(3)}.png", transparent=True, dpi=100
+    )
 
 # gif:
 image_filenames = [
@@ -86,7 +89,7 @@ def update_graph(num):
     format="GIF",
     save_all=True,
     append_images=images[1:],
-    duration=20,
+    duration=40,
     loop=0,
     transparency=0,
     disposal=2,