"""Generate the original visual atlas for the Dynamical Systems course. Every result is deterministic and is derived from the equations stated in the course. The script intentionally keeps the numerical method explicit so each figure can be reproduced or adapted without a proprietary plotting pipeline. Inside this repository the assets are written to their canonical public path. A downloaded copy instead writes to ``learning-figures/dinamica-caos`` below the current working directory, so it remains executable as a standalone file. """ from __future__ import annotations from pathlib import Path from typing import Callable import matplotlib as mpl mpl.use("Agg") import matplotlib.pyplot as plt import numpy as np from matplotlib.animation import FuncAnimation, PillowWriter from matplotlib.collections import LineCollection from scipy.signal import welch from scipy.special import gamma from scipy.spatial.distance import pdist, squareform SCRIPT_ROOT = Path(__file__).resolve().parents[1] IN_REPOSITORY = (SCRIPT_ROOT / "src" / "content" / "learning").is_dir() BASE_OUT = ( SCRIPT_ROOT / "public" / "images" / "learning" / "dinamica-caos" if IN_REPOSITORY else Path.cwd() / "learning-figures" / "dinamica-caos" ) OUT = BASE_OUT THEMES = { "dark": { "bg": "#07131f", "panel": "#0d2233", "grid": "#365065", "text": "#edf6ff", "muted": "#a8bfd0", "cyan": "#2dd4bf", "blue": "#60a5fa", "magenta": "#f472b6", "orange": "#fb923c", "yellow": "#facc15", "purple": "#a78bfa", }, "light": { "bg": "#ffffff", "panel": "#f8fafc", "grid": "#cbd5e1", "text": "#172033", "muted": "#475569", "cyan": "#0f766e", "blue": "#1d4ed8", "magenta": "#be185d", "orange": "#c2410c", "yellow": "#a16207", "purple": "#6d28d9", }, } BG = PANEL = GRID = TEXT = MUTED = CYAN = BLUE = MAGENTA = ORANGE = YELLOW = PURPLE = "" def select_theme(name: str) -> None: global OUT, BG, PANEL, GRID, TEXT, MUTED, CYAN, BLUE, MAGENTA, ORANGE, YELLOW, PURPLE palette = THEMES[name] OUT = BASE_OUT if name == "dark" else BASE_OUT / name BG = palette["bg"] PANEL = palette["panel"] GRID = palette["grid"] TEXT = palette["text"] MUTED = palette["muted"] CYAN = palette["cyan"] BLUE = palette["blue"] MAGENTA = palette["magenta"] ORANGE = palette["orange"] YELLOW = palette["yellow"] PURPLE = palette["purple"] def configure_style() -> None: mpl.rcParams.update( { "figure.facecolor": BG, "savefig.facecolor": BG, "axes.facecolor": PANEL, "axes.edgecolor": GRID, "axes.labelcolor": TEXT, "xtick.color": MUTED, "ytick.color": MUTED, "text.color": TEXT, "grid.color": GRID, "grid.alpha": 0.35, "font.family": "DejaVu Sans", "font.size": 11, } ) def style_3d_panes(ax) -> None: for axis in (ax.xaxis, ax.yaxis, ax.zaxis): axis.pane.set_facecolor(PANEL) axis.pane.set_edgecolor(GRID) axis.pane.set_alpha(1.0) ax.tick_params(colors=MUTED) def save(fig: mpl.figure.Figure, filename: str) -> None: fig.savefig(OUT / filename, dpi=190, bbox_inches="tight", pad_inches=0.18) plt.close(fig) def rk4( field: Callable[[np.ndarray], np.ndarray], x0: np.ndarray, dt: float, steps: int, ) -> np.ndarray: trajectory = np.empty((steps + 1, len(x0)), dtype=float) trajectory[0] = x0 for index in range(steps): x = trajectory[index] k1 = field(x) k2 = field(x + 0.5 * dt * k1) k3 = field(x + 0.5 * dt * k2) k4 = field(x + dt * k3) trajectory[index + 1] = x + dt * (k1 + 2 * k2 + 2 * k3 + k4) / 6 return trajectory def colored_line_2d(ax, x: np.ndarray, y: np.ndarray, cmap: str = "turbo") -> None: points = np.column_stack([x, y]).reshape(-1, 1, 2) segments = np.concatenate([points[:-1], points[1:]], axis=1) collection = LineCollection(segments, cmap=cmap, linewidth=2.5, alpha=0.94) collection.set_array(np.linspace(0, 1, len(segments))) ax.add_collection(collection) ax.autoscale() def lorenz_field(state: np.ndarray) -> np.ndarray: x, y, z = state return np.array([10.0 * (y - x), x * (28.0 - z) - y, x * y - (8.0 / 3.0) * z]) def sprott_b_field(state: np.ndarray) -> np.ndarray: x, y, z = state return np.array([y * z, x - y, 1.0 - x * y]) def rossler_field(state: np.ndarray) -> np.ndarray: x, y, z = state return np.array([-y - z, x + 0.2 * y, 0.2 + z * (x - 5.7)]) def hopf_field(mu: float) -> Callable[[np.ndarray], np.ndarray]: def field(state: np.ndarray) -> np.ndarray: x, y = state radius_squared = x * x + y * y return np.array([ mu * x - y - radius_squared * x, x + mu * y - radius_squared * y, ]) return field def double_well_field(state: np.ndarray) -> np.ndarray: x, velocity = state return np.array([velocity, x - x**3 - 0.45 * velocity]) def integrate_decay(method: str, step: float, end: float = 2.0) -> tuple[np.ndarray, np.ndarray]: count = int(round(end / step)) time = np.linspace(0.0, end, count + 1) values = np.empty(count + 1) values[0] = 1.0 for index in range(count): value = values[index] if method == "euler": values[index + 1] = value + step * (-2.0 * value) elif method == "heun": predictor = value + step * (-2.0 * value) values[index + 1] = value + 0.5 * step * (-2.0 * value - 2.0 * predictor) elif method == "rk4": k1 = -2.0 * value k2 = -2.0 * (value + 0.5 * step * k1) k3 = -2.0 * (value + 0.5 * step * k2) k4 = -2.0 * (value + step * k3) values[index + 1] = value + step * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0 else: raise ValueError(f"Unknown integration method: {method}") return time, values def classify_double_well_basins( x_values: np.ndarray, velocity_values: np.ndarray, steps: int = 650, dt: float = 0.032, ) -> np.ndarray: xx, vv = np.meshgrid(x_values, velocity_values) x = xx.copy() velocity = vv.copy() for _ in range(steps): k1x = velocity k1v = x - x**3 - 0.45 * velocity x2 = x + 0.5 * dt * k1x v2 = velocity + 0.5 * dt * k1v k2x = v2 k2v = x2 - x2**3 - 0.45 * v2 x3 = x + 0.5 * dt * k2x v3 = velocity + 0.5 * dt * k2v k3x = v3 k3v = x3 - x3**3 - 0.45 * v3 x4 = x + dt * k3x v4 = velocity + dt * k3v k4x = v4 k4v = x4 - x4**3 - 0.45 * v4 x += dt * (k1x + 2 * k2x + 2 * k3x + k4x) / 6.0 velocity += dt * (k1v + 2 * k2v + 2 * k3v + k4v) / 6.0 return (x >= 0.0).astype(float) def hopf_return_map(radius: np.ndarray, mu: float = 0.12, period: float = 2 * np.pi) -> np.ndarray: squared = radius**2 denominator = squared + (mu - squared) * np.exp(-2.0 * mu * period) return np.sqrt(np.maximum(mu * squared / np.maximum(denominator, 1e-14), 0.0)) def abm_pece_relaxation(q: float, step: float, end: float = 3.0) -> tuple[np.ndarray, np.ndarray]: count = int(round(end / step)) time = np.linspace(0.0, end, count + 1) values = np.empty(count + 1) force = np.empty(count + 1) values[0] = 1.0 force[0] = -values[0] predictor_scale = step**q / gamma(q + 1.0) corrector_scale = step**q / gamma(q + 2.0) for n in range(count): j = np.arange(n + 1) predictor_weights = (n + 1 - j) ** q - (n - j) ** q predicted = 1.0 + predictor_scale * np.dot(predictor_weights, force[: n + 1]) corrector_weights = np.empty(n + 1) corrector_weights[0] = n ** (q + 1.0) - (n - q) * (n + 1.0) ** q if n >= 1: interior = np.arange(1, n + 1) age = n - interior corrector_weights[1:] = ( (age + 2.0) ** (q + 1.0) + age ** (q + 1.0) - 2.0 * (age + 1.0) ** (q + 1.0) ) values[n + 1] = 1.0 + corrector_scale * ( -predicted + np.dot(corrector_weights, force[: n + 1]) ) force[n + 1] = -values[n + 1] return time, values def plot_flow_phase_space() -> None: x = np.linspace(-2.8 * np.pi, 2.8 * np.pi, 33) v = np.linspace(-3.2, 3.2, 27) xx, vv = np.meshgrid(x, v) dx = vv dv = -0.35 * vv - np.sin(xx) speed = np.hypot(dx, dv) fig, ax = plt.subplots(figsize=(10.4, 6.2)) stream = ax.streamplot( xx, vv, dx, dv, color=speed, cmap="viridis", density=1.25, linewidth=0.8, arrowsize=0.8, ) for theta0, v0, color in [(-2.4, 0.2, MAGENTA), (1.8, 2.2, ORANGE), (5.0, -1.1, CYAN)]: states = rk4( lambda state: np.array([state[1], -0.35 * state[1] - np.sin(state[0])]), np.array([theta0, v0]), 0.025, 950, ) ax.plot(states[:, 0], states[:, 1], color=color, lw=2.5, alpha=0.98) ax.scatter(states[0, 0], states[0, 1], color=color, s=34, zorder=5) ax.scatter([-2 * np.pi, 0, 2 * np.pi], [0, 0, 0], color=YELLOW, s=45, edgecolor=BG) ax.set(xlabel=r"Ángulo $\theta$", ylabel=r"Velocidad $\dot\theta$", xlim=(x.min(), x.max()), ylim=(v.min(), v.max())) fig.colorbar(stream.lines, ax=ax, label=r"Rapidez $\|f(\theta,\dot\theta)\|$") save(fig, "01_flujo_espacio_fases.png") def plot_saddle_manifolds() -> None: x = np.linspace(-3, 3, 31) y = np.linspace(-3, 3, 31) xx, yy = np.meshgrid(x, y) fig, ax = plt.subplots(figsize=(8.6, 6.4)) ax.streamplot(xx, yy, xx, -yy, color=np.hypot(xx, yy), cmap="magma", density=1.25, linewidth=0.8) ax.axhline(0, color=ORANGE, lw=3, label=r"$W^u(0)=\{y=0\}$") ax.axvline(0, color=CYAN, lw=3, label=r"$W^s(0)=\{x=0\}$") ax.scatter([0], [0], color=YELLOW, s=75, zorder=6, edgecolor=BG) ax.set(xlabel="$x$", ylabel="$y$", xlim=(-3, 3), ylim=(-3, 3)) ax.legend(loc="upper left", frameon=True, facecolor=BG, edgecolor=GRID) ax.grid(True) save(fig, "02_silla_variedades.png") def logistic_cobweb(ax, r: float, x0: float, steps: int) -> None: grid = np.linspace(0, 1, 600) curve = r * grid * (1 - grid) ax.plot(grid, curve, color=CYAN, lw=2.1, label=r"$f_r(x)=rx(1-x)$") ax.plot(grid, grid, color=MUTED, lw=1.25, ls="--", label=r"$y=x$") x = x0 px, py = [x], [0.0] for _ in range(steps): y = r * x * (1 - x) px.extend([x, y]) py.extend([y, y]) x = y ax.plot(px, py, color=MAGENTA, lw=1.0, alpha=0.86) ax.scatter([x0], [0], color=YELLOW, s=32, zorder=5) ax.set(xlim=(0, 1), ylim=(0, 1), xlabel="$x_n$", ylabel="$x_{n+1}$") ax.grid(True) def plot_logistic_cobweb() -> None: fig, axes = plt.subplots(1, 2, figsize=(11.2, 5.2), sharex=True, sharey=True) logistic_cobweb(axes[0], 3.2, 0.18, 36) logistic_cobweb(axes[1], 3.9, 0.18, 72) axes[0].legend(loc="upper left", frameon=True, facecolor=BG, edgecolor=GRID) save(fig, "03_logistica_telarana.png") def plot_logistic_bifurcation() -> None: r = np.linspace(2.5, 4.0, 3600) x = np.full_like(r, 0.217) lyap = np.zeros_like(r) samples: list[np.ndarray] = [] for index in range(1250): x = r * x * (1 - x) if index >= 700: derivative = np.maximum(np.abs(r * (1 - 2 * x)), 1e-14) lyap += np.log(derivative) if index % 3 == 0: samples.append(x.copy()) lyap /= 550 fig, (top, bottom) = plt.subplots(2, 1, figsize=(11.2, 7.7), sharex=True, gridspec_kw={"height_ratios": [2.25, 1]}) rr = np.tile(r, len(samples)) xx = np.concatenate(samples) top.scatter(rr, xx, s=0.045, color=CYAN, alpha=0.43, rasterized=True) top.set(ylabel=r"Estados asintóticos $x_n$", ylim=(0, 1)) bottom.plot(r, lyap, color=MAGENTA, lw=0.9) bottom.axhline(0, color=YELLOW, lw=1.1, ls="--") bottom.fill_between(r, 0, lyap, where=lyap > 0, color=MAGENTA, alpha=0.25) bottom.set(xlabel="Parámetro $r$", ylabel=r"$\lambda$") for ax in (top, bottom): ax.grid(True) save(fig, "04_logistica_bifurcacion.png") def lorenz_trajectory() -> tuple[np.ndarray, float]: dt = 0.005 trajectory = rk4(lorenz_field, np.array([1.0, 1.0, 1.0]), dt, 34000) return trajectory[4000:], dt def plot_lorenz(trajectory: np.ndarray) -> None: fig = plt.figure(figsize=(10, 7.1)) ax = fig.add_subplot(111, projection="3d") stride = 2 values = np.linspace(0, 1, len(trajectory[::stride])) ax.plot(trajectory[:, 0], trajectory[:, 1], trajectory[:, 2], color=CYAN, lw=0.48, alpha=0.72) ax.scatter( trajectory[::stride, 0], trajectory[::stride, 1], trajectory[::stride, 2], c=values, cmap="turbo", s=0.42, alpha=0.90, ) ax.set(xlabel="$x$", ylabel="$y$", zlabel="$z$") ax.view_init(24, -58) ax.grid(False) style_3d_panes(ax) save(fig, "05_lorenz_attractor.png") def plot_poincare(trajectory: np.ndarray) -> None: plane = 27.0 z0 = trajectory[:-1, 2] z1 = trajectory[1:, 2] mask = (z0 < plane) & (z1 >= plane) alpha = (plane - z0[mask]) / (z1[mask] - z0[mask]) crossings = trajectory[:-1][mask] + alpha[:, None] * (trajectory[1:][mask] - trajectory[:-1][mask]) fig = plt.figure(figsize=(11.2, 5.3)) ax3 = fig.add_subplot(121, projection="3d") sample = trajectory[::5] ax3.plot(sample[:, 0], sample[:, 1], sample[:, 2], color=BLUE, lw=0.62, alpha=0.82) gx, gy = np.meshgrid(np.linspace(-22, 22, 8), np.linspace(-30, 30, 8)) ax3.plot_surface(gx, gy, np.full_like(gx, plane), color=MAGENTA, alpha=0.22, edgecolor="none") ax3.scatter(crossings[:, 0], crossings[:, 1], crossings[:, 2], color=YELLOW, s=10) ax3.set(xlabel="$x$", ylabel="$y$", zlabel="$z$") ax3.view_init(24, -61) style_3d_panes(ax3) ax2 = fig.add_subplot(122) ax2.scatter(crossings[:, 0], crossings[:, 1], c=np.arange(len(crossings)), cmap="turbo", s=11, alpha=0.82) ax2.set(xlabel="$x$ al cruzar", ylabel="$y$ al cruzar") ax2.grid(True) save(fig, "06_poincare_section.png") def plot_lyapunov_divergence() -> None: dt = 0.0025 steps = 8000 first = rk4(lorenz_field, np.array([1.0, 1.0, 1.0]), dt, steps) second = rk4(lorenz_field, np.array([1.0 + 1e-8, 1.0, 1.0]), dt, steps) time = np.arange(steps + 1) * dt distance = np.linalg.norm(second - first, axis=1) distance = np.maximum(distance, 1e-16) fit_mask = (time >= 1.0) & (time <= 9.0) & (distance < 0.5) slope, intercept = np.polyfit(time[fit_mask], np.log(distance[fit_mask]), 1) reference = np.exp(intercept + slope * time) fig, ax = plt.subplots(figsize=(10.4, 6.2)) ax.semilogy(time, distance, color=CYAN, lw=1.45, label=r"$\|\delta x(t)\|$") ax.semilogy(time[fit_mask], reference[fit_mask], color=YELLOW, lw=2, ls="--", label=fr"ajuste inicial $\lambda\approx{slope:.2f}$") ax.axhspan(1, max(1.1, distance.max()), color=MAGENTA, alpha=0.08, label="saturación geométrica") ax.set(xlabel="Tiempo", ylabel="Separación entre trayectorias", ylim=(1e-9, 1e2)) ax.grid(True, which="both") ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) save(fig, "07_lyapunov_divergence.png") def plot_spectrum(trajectory: np.ndarray, dt: float) -> None: signal = trajectory[:, 0] - np.mean(trajectory[:, 0]) frequency, density = welch(signal, fs=1 / dt, window="hann", nperseg=4096, noverlap=2048) time = np.arange(len(signal)) * dt fig, (left, right) = plt.subplots(1, 2, figsize=(11.3, 5.1)) left.plot(time[:5000], signal[:5000], color=CYAN, lw=0.7) left.set(xlabel="Tiempo", ylabel=r"$x(t)-\bar x$") right.semilogy(frequency[1:500], density[1:500], color=ORANGE, lw=1.35) right.fill_between(frequency[1:500], density[1:500], color=ORANGE, alpha=0.18) right.set(xlabel="Frecuencia", ylabel="PSD") for ax in (left, right): ax.grid(True) save(fig, "08_fourier_spectrum.png") def plot_caputo_kernel() -> None: lag = np.geomspace(0.025, 12, 700) fig, ax = plt.subplots(figsize=(10.2, 6.1)) for q, color in [(0.3, MAGENTA), (0.55, ORANGE), (0.8, CYAN)]: kernel = lag ** (q - 1) / gamma(q) ax.loglog(lag, kernel, color=color, lw=2.4, label=fr"$q={q}$") ax.set(xlabel=r"Antigüedad $t-s$", ylabel=r"Peso $(t-s)^{q-1}/\Gamma(q)$") ax.grid(True, which="both") ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) save(fig, "09_caputo_memory_kernel.png") def plot_sprott_b() -> None: trajectory = rk4(sprott_b_field, np.array([0.1, 0.1, 0.1]), 0.006, 65000)[5000:] fig = plt.figure(figsize=(10, 7.1)) ax = fig.add_subplot(111, projection="3d") sample = trajectory[::2] ax.plot(sample[:, 0], sample[:, 1], sample[:, 2], color=ORANGE, lw=0.62, alpha=0.80) ax.scatter(sample[:, 0], sample[:, 1], sample[:, 2], c=np.linspace(0, 1, len(sample)), cmap="plasma", s=0.40, alpha=0.92) ax.set(xlabel="$x$", ylabel="$y$", zlabel="$z$") ax.view_init(22, -48) ax.grid(False) style_3d_panes(ax) save(fig, "10_sprott_b.png") def plot_pendulum_potential_phase() -> None: theta = np.linspace(-2.2 * np.pi, 2.2 * np.pi, 1200) potential = 1.0 - np.cos(theta) energies = [(0.65, CYAN), (1.55, ORANGE), (2.35, MAGENTA)] fig, (left, right) = plt.subplots(1, 2, figsize=(11.4, 5.2)) left.plot(theta, potential, color=BLUE, lw=3.0) for energy, color in energies: left.axhline(energy, color=color, lw=1.8, ls="--") left.set(xlabel=r"$\theta$", ylabel=r"$V(\theta)=1-\cos\theta$", xlim=(theta.min(), theta.max()), ylim=(-0.08, 2.55)) left.grid(True) for energy, color in energies: radicand = 2.0 * (energy - potential) valid = radicand >= 0.0 upper = np.where(valid, np.sqrt(np.maximum(radicand, 0.0)), np.nan) right.plot(theta, upper, color=color, lw=2.2) right.plot(theta, -upper, color=color, lw=2.2) damped = rk4( lambda state: np.array([state[1], -0.24 * state[1] - np.sin(state[0])]), np.array([1.85 * np.pi, 0.45]), 0.022, 1500, ) right.plot(damped[:, 0], damped[:, 1], color=YELLOW, lw=2.7, alpha=0.95) right.scatter([damped[0, 0]], [damped[0, 1]], color=YELLOW, s=38, edgecolor=BG, zorder=5) right.set(xlabel=r"$\theta$", ylabel=r"$\dot\theta$", xlim=(theta.min(), theta.max()), ylim=(-2.45, 2.45)) right.grid(True) save(fig, "11_pendulo_potencial_fase.png") def plot_pitchfork_bifurcation() -> None: negative = np.linspace(-1.25, 0.0, 360) positive = np.linspace(0.0, 1.45, 440) branch = np.sqrt(positive) fig, ax = plt.subplots(figsize=(9.6, 6.0)) ax.plot(negative, np.zeros_like(negative), color=CYAN, lw=3.2, label=r"$x^*=0$") ax.plot(positive, np.zeros_like(positive), color=MAGENTA, lw=2.5, ls="--") ax.plot(positive, branch, color=CYAN, lw=3.2, label=r"$x^*=\pm\sqrt{\mu}$") ax.plot(positive, -branch, color=CYAN, lw=3.2) ax.axvline(0.0, color=GRID, lw=1.2) ax.axhline(0.0, color=GRID, lw=1.2) ax.set(xlabel=r"$\mu$", ylabel=r"Equilibrio $x^*$", xlim=(-1.25, 1.45), ylim=(-1.35, 1.35)) ax.grid(True) ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) save(fig, "12_bifurcacion_horquilla.png") def plot_trace_determinant_plane() -> None: trace = np.linspace(-4.0, 4.0, 700) determinant = np.linspace(-2.0, 5.0, 620) tt, dd = np.meshgrid(trace, determinant) regions = np.zeros_like(tt, dtype=int) regions[dd < 0.0] = 0 node = (dd >= 0.0) & (dd <= tt**2 / 4.0) focus = dd > tt**2 / 4.0 regions[node & (tt < 0.0)] = 1 regions[focus & (tt < 0.0)] = 2 regions[node & (tt >= 0.0)] = 3 regions[focus & (tt >= 0.0)] = 4 colors = [MAGENTA, CYAN, BLUE, ORANGE, PURPLE] cmap = mpl.colors.ListedColormap(colors) fig, ax = plt.subplots(figsize=(9.4, 6.2)) ax.pcolormesh(trace, determinant, regions, shading="auto", cmap=cmap, alpha=0.36) ax.plot(trace, trace**2 / 4.0, color=YELLOW, lw=2.8, label=r"$\Delta=\tau^2/4$") ax.axhline(0.0, color=TEXT, lw=2.0, label=r"$\Delta=0$") ax.axvline(0.0, color=TEXT, lw=1.3, ls="--", label=r"$\tau=0$") ax.scatter([-2.8, -1.2, 2.3, 1.0, 0.0], [1.0, 2.2, 1.1, 2.4, 0.0], c=[CYAN, BLUE, ORANGE, PURPLE, YELLOW], s=55, edgecolor=BG, zorder=5) ax.set(xlabel=r"Traza $\tau$", ylabel=r"Determinante $\Delta$", xlim=(-4.0, 4.0), ylim=(-2.0, 5.0)) ax.grid(True) ax.legend(frameon=True, facecolor=BG, edgecolor=GRID, loc="upper center", ncol=3) save(fig, "13_plano_traza_determinante.png") def plot_nullclines_field() -> None: x = np.linspace(-1.9, 1.9, 35) velocity = np.linspace(-1.8, 1.8, 33) xx, vv = np.meshgrid(x, velocity) dx = vv dv = xx - xx**3 - 0.45 * vv speed = np.hypot(dx, dv) curve_x = np.linspace(-1.75, 1.75, 900) fig, ax = plt.subplots(figsize=(9.5, 6.2)) ax.streamplot(xx, vv, dx, dv, color=speed, cmap="viridis", density=1.3, linewidth=0.9, arrowsize=0.85) ax.axhline(0.0, color=CYAN, lw=3.0, label=r"$\dot x=0$") ax.plot(curve_x, (curve_x - curve_x**3) / 0.45, color=MAGENTA, lw=3.0, label=r"$\dot v=0$") ax.scatter([-1.0, 0.0, 1.0], [0.0, 0.0, 0.0], color=YELLOW, edgecolor=BG, s=65, zorder=6) ax.set(xlabel="$x$", ylabel="$v$", xlim=(-1.9, 1.9), ylim=(-1.8, 1.8)) ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) save(fig, "14_nullclines_campo.png") def plot_integrator_comparison() -> None: dense_time = np.linspace(0.0, 2.0, 900) exact = np.exp(-2.0 * dense_time) fig, (left, right) = plt.subplots(1, 2, figsize=(11.5, 5.2)) left.plot(dense_time, exact, color=TEXT, lw=3.0, label=r"$e^{-2t}$") methods = [("euler", MAGENTA, "o", "Euler"), ("heun", ORANGE, "s", "Heun"), ("rk4", CYAN, "^", "RK4")] for method, color, marker, label in methods: time, values = integrate_decay(method, 0.4) left.plot(time, values, color=color, marker=marker, ms=5, lw=2.2, label=label) error = np.abs(values - np.exp(-2.0 * time)) right.plot(time[1:], error[1:], color=color, marker=marker, ms=5, lw=2.2, label=label) left.set(xlabel="$t$", ylabel="$y(t)$", xlim=(0.0, 2.0)) right.set(xlabel="$t$", ylabel=r"$|y_h(t)-y(t)|$", xlim=(0.0, 2.0)) right.set_yscale("log") right.set_ylim(1e-7, 1.0) for ax in (left, right): ax.grid(True, which="both") ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) save(fig, "15_integradores_comparacion.png") def plot_step_convergence() -> None: steps = np.array([0.4, 0.2, 0.1, 0.05, 0.025]) methods = [("euler", MAGENTA, "o", "Euler"), ("heun", ORANGE, "s", "Heun"), ("rk4", CYAN, "^", "RK4")] fig, ax = plt.subplots(figsize=(9.4, 6.1)) for method, color, marker, label in methods: errors = [] for step in steps: time, values = integrate_decay(method, float(step)) errors.append(np.max(np.abs(values - np.exp(-2.0 * time)))) ax.loglog(steps, errors, color=color, marker=marker, ms=7, lw=2.7, label=label) reference = 0.18 * (steps / steps[0]) ax.loglog(steps, reference, color=MUTED, lw=1.4, ls="--", label=r"$O(h)$") ax.loglog(steps, reference * (steps / steps[0]), color=BLUE, lw=1.4, ls="--", label=r"$O(h^2)$") ax.loglog(steps, 0.015 * (steps / steps[0]) ** 4, color=YELLOW, lw=1.4, ls="--", label=r"$O(h^4)$") ax.set(xlabel="Paso $h$", ylabel="Error global máximo") ax.grid(True, which="both") ax.legend(frameon=True, facecolor=BG, edgecolor=GRID, ncol=2) ax.invert_xaxis() save(fig, "16_convergencia_paso.png") def plot_hopf_bifurcation() -> None: negative = np.linspace(-0.6, 0.0, 260) positive = np.linspace(0.0, 0.8, 360) fig, (left, right) = plt.subplots(1, 2, figsize=(11.4, 5.3)) left.plot(negative, np.zeros_like(negative), color=CYAN, lw=3.0) left.plot(positive, np.zeros_like(positive), color=MAGENTA, lw=2.3, ls="--") left.plot(positive, np.sqrt(positive), color=CYAN, lw=3.0) left.set(xlabel=r"$\mu$", ylabel=r"Amplitud $r$", xlim=(-0.6, 0.8), ylim=(-0.03, 0.98)) left.grid(True) mu = 0.28 circle = np.linspace(0.0, 2.0 * np.pi, 600) right.plot(np.sqrt(mu) * np.cos(circle), np.sqrt(mu) * np.sin(circle), color=YELLOW, lw=2.6, ls="--") for initial, color in [([0.12, 0.08], CYAN), ([1.05, 0.10], ORANGE), ([-0.75, 0.65], MAGENTA)]: trajectory = rk4(hopf_field(mu), np.array(initial, dtype=float), 0.025, 900) right.plot(trajectory[:, 0], trajectory[:, 1], color=color, lw=2.1) right.scatter([trajectory[0, 0]], [trajectory[0, 1]], color=color, s=28, edgecolor=BG) right.set(xlabel="$x$", ylabel="$y$", xlim=(-1.15, 1.15), ylim=(-1.15, 1.15), aspect="equal") right.grid(True) save(fig, "17_bifurcacion_hopf.png") def henon_step(points: np.ndarray) -> np.ndarray: x = points[:, 0] y = points[:, 1] return np.column_stack([1.0 - 1.4 * x**2 + y, 0.3 * x]) def plot_henon_folding() -> None: grid_x, grid_y = np.meshgrid(np.linspace(-1.1, 1.1, 64), np.linspace(-0.32, 0.32, 27)) points = np.column_stack([grid_x.ravel(), grid_y.ravel()]) colors = np.tile(np.linspace(0.0, 1.0, grid_x.shape[1]), grid_x.shape[0]) stages = [points] for _ in range(2): stages.append(henon_step(stages[-1])) fig, axes = plt.subplots(1, 3, figsize=(12.0, 4.2)) for ax, stage in zip(axes, stages): ax.scatter(stage[:, 0], stage[:, 1], c=colors, cmap="turbo", s=5.5, alpha=0.82) ax.set(xlabel="$x$", ylabel="$y$", xlim=(-1.8, 1.8), ylim=(-0.62, 0.62)) ax.grid(True) save(fig, "18_henon_plegamiento.png") def plot_bistable_basin() -> None: x_values = np.linspace(-2.1, 2.1, 270) velocity_values = np.linspace(-2.0, 2.0, 230) basins = classify_double_well_basins(x_values, velocity_values) cmap = mpl.colors.ListedColormap([BLUE, MAGENTA]) fig, ax = plt.subplots(figsize=(9.6, 6.2)) ax.imshow( basins, origin="lower", extent=(x_values.min(), x_values.max(), velocity_values.min(), velocity_values.max()), cmap=cmap, interpolation="nearest", alpha=0.82, aspect="auto", ) stable_eigenvalue = (-0.45 - np.sqrt(0.45**2 + 4.0)) / 2.0 for sign in (-1.0, 1.0): manifold = rk4(double_well_field, np.array([sign * 1e-4, sign * stable_eigenvalue * 1e-4]), -0.009, 900) mask = (np.abs(manifold[:, 0]) <= 2.1) & (np.abs(manifold[:, 1]) <= 2.0) ax.plot(manifold[mask, 0], manifold[mask, 1], color=YELLOW, lw=3.0) ax.scatter([-1.0, 0.0, 1.0], [0.0, 0.0, 0.0], color=[BLUE, YELLOW, MAGENTA], edgecolor=BG, s=70, zorder=5) ax.set(xlabel="$x_0$", ylabel="$v_0$", xlim=(-2.1, 2.1), ylim=(-2.0, 2.0)) save(fig, "19_cuenca_biestable.png") def plot_return_floquet() -> None: mu = 0.12 radius = np.linspace(0.0, 0.75, 700) returned = hopf_return_map(radius, mu) fixed = np.sqrt(mu) parameters = np.linspace(0.01, 0.55, 600) multiplier = np.exp(-4.0 * np.pi * parameters) fig, (left, right) = plt.subplots(1, 2, figsize=(11.3, 5.1)) left.plot(radius, returned, color=CYAN, lw=3.0, label=r"$P(r)$") left.plot(radius, radius, color=MUTED, lw=1.8, ls="--", label=r"$r$") left.scatter([fixed], [fixed], color=YELLOW, edgecolor=BG, s=75, zorder=5) left.set(xlabel="$r_n$", ylabel="$r_{n+1}$", xlim=(0.0, 0.75), ylim=(0.0, 0.75)) right.plot(parameters, multiplier, color=ORANGE, lw=3.0, label=r"$\rho=e^{-4\pi\mu}$") right.axhline(1.0, color=MUTED, lw=1.8, ls="--") right.set(xlabel=r"$\mu$", ylabel=r"Multiplicador $|\rho|$", xlim=(0.0, 0.55), ylim=(0.0, 1.06)) for ax in (left, right): ax.grid(True) ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) save(fig, "20_retorno_floquet.png") def plot_topological_degree() -> None: angle = np.linspace(0.0, 2.0 * np.pi, 1000) domain_x = np.cos(angle) domain_y = np.sin(angle) image_x = np.cos(2.0 * angle) + 0.35 * np.cos(angle) image_y = np.sin(2.0 * angle) + 0.35 * np.sin(angle) fig, (left, right) = plt.subplots(1, 2, figsize=(10.8, 5.1)) colored_line_2d(left, domain_x, domain_y) colored_line_2d(right, image_x, image_y) sample = np.arange(0, len(angle) - 1, 125) left.quiver(domain_x[sample], domain_y[sample], -domain_y[sample], domain_x[sample], color=YELLOW, angles="xy", scale_units="xy", scale=5.0, width=0.012) derivatives = np.column_stack([ -2.0 * np.sin(2.0 * angle) - 0.35 * np.sin(angle), 2.0 * np.cos(2.0 * angle) + 0.35 * np.cos(angle), ]) right.quiver(image_x[sample], image_y[sample], derivatives[sample, 0], derivatives[sample, 1], color=YELLOW, angles="xy", scale_units="xy", scale=9.0, width=0.012) for ax in (left, right): ax.scatter([0.0], [0.0], color=MAGENTA, s=55, edgecolor=BG, zorder=5) ax.set(xlabel="$x$", ylabel="$y$", xlim=(-1.55, 1.55), ylim=(-1.55, 1.55), aspect="equal") ax.grid(True) save(fig, "21_grado_topologico.png") def recurrence_data(trajectory: np.ndarray, count: int = 520) -> tuple[np.ndarray, np.ndarray, float]: sample = trajectory[:: max(1, len(trajectory) // count)][:count] normalized = (sample - np.mean(sample, axis=0)) / np.std(sample, axis=0) distances = squareform(pdist(normalized)) off_diagonal = distances[~np.eye(len(distances), dtype=bool)] threshold = float(np.quantile(off_diagonal, 0.085)) return sample, distances <= threshold, threshold def plot_recurrence(trajectory: np.ndarray, dt: float) -> None: sample, recurrence, _ = recurrence_data(trajectory) sample_stride = max(1, len(trajectory) // len(sample)) time = np.arange(len(sample)) * dt * sample_stride cmap = mpl.colors.ListedColormap([BG, CYAN]) fig, (left, right) = plt.subplots(1, 2, figsize=(11.3, 5.1)) left.plot(time, sample[:, 0], color=ORANGE, lw=1.25) left.set(xlabel="Tiempo", ylabel="$x(t)$") left.grid(True) right.imshow(recurrence, origin="lower", cmap=cmap, interpolation="nearest", aspect="equal") right.set(xlabel="$i$", ylabel="$j$") save(fig, "22_recurrencia.png") def plot_correlation_dimension(trajectory: np.ndarray) -> None: sample = trajectory[::20][:1500] normalized = (sample - np.mean(sample, axis=0)) / np.std(sample, axis=0) distances = pdist(normalized) positive = distances[distances > 0.0] epsilon = np.geomspace(np.quantile(positive, 0.003), np.quantile(positive, 0.34), 90) sorted_distances = np.sort(positive) correlation = np.searchsorted(sorted_distances, epsilon, side="right") / len(sorted_distances) log_epsilon = np.log(epsilon) log_correlation = np.log(np.maximum(correlation, 1e-14)) fit_mask = (correlation >= 0.012) & (correlation <= 0.16) slope, intercept = np.polyfit(log_epsilon[fit_mask], log_correlation[fit_mask], 1) local_slope = np.gradient(log_correlation, log_epsilon) fig, (left, right) = plt.subplots(1, 2, figsize=(11.2, 5.1)) left.loglog(epsilon, correlation, color=CYAN, lw=2.8, label=r"$C(\varepsilon)$") left.loglog(epsilon[fit_mask], np.exp(intercept + slope * log_epsilon[fit_mask]), color=YELLOW, lw=2.4, ls="--", label=fr"$D_2\approx{slope:.2f}$") left.set(xlabel=r"$\varepsilon$", ylabel=r"$C(\varepsilon)$") left.legend(frameon=True, facecolor=BG, edgecolor=GRID) right.semilogx(epsilon, local_slope, color=ORANGE, lw=2.6) right.axhline(slope, color=YELLOW, lw=1.8, ls="--") right.fill_between(epsilon, slope - 0.18, slope + 0.18, color=YELLOW, alpha=0.13) right.set(xlabel=r"$\varepsilon$", ylabel=r"$d\log C/d\log\varepsilon$", ylim=(0.0, 4.5)) for ax in (left, right): ax.grid(True, which="both") save(fig, "23_dimension_correlacion.png") def plot_rl_caputo_powers() -> None: time = np.geomspace(0.025, 4.0, 800) q = 0.65 rl_constant = time ** (-q) / gamma(1.0 - q) power_derivative = gamma(3.0) / gamma(3.0 - q) * time ** (2.0 - q) fig, (left, right) = plt.subplots(1, 2, figsize=(11.2, 5.1)) left.loglog(time, rl_constant, color=MAGENTA, lw=3.0, label=r"${}^{RL}D^q1$") left.plot(time, np.full_like(time, 1e-3), color=CYAN, lw=2.7, ls="--", label=r"${}^{C}D^q1=0$") left.set(xlabel="$t$", ylabel=r"$D^q1$", ylim=(7e-4, 6.0)) right.loglog(time, power_derivative, color=ORANGE, lw=3.0, label=r"$D^q t^2$") right.loglog(time, time**2, color=MUTED, lw=1.8, ls="--", label=r"$t^2$") right.set(xlabel="$t$", ylabel=r"$D^q t^2$") for ax in (left, right): ax.grid(True, which="both") ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) save(fig, "24_rl_caputo_potencias.png") def plot_abm_convergence() -> None: q = 0.72 dense_time = np.linspace(0.0, 3.0, 900) exact = mittag_leffler_relaxation(q, dense_time) steps = np.array([0.25, 0.125, 0.0625, 0.03125]) errors = [] fig, (left, right) = plt.subplots(1, 2, figsize=(11.4, 5.1)) left.plot(dense_time, exact, color=TEXT, lw=3.0, label=r"$E_q(-t^q)$") for step, color in zip(steps[:3], [MAGENTA, ORANGE, CYAN]): time, values = abm_pece_relaxation(q, float(step)) left.plot(time, values, color=color, marker="o", ms=3.5, lw=1.9, label=fr"$h={step:g}$") for step in steps: time, values = abm_pece_relaxation(q, float(step)) reference = mittag_leffler_relaxation(q, time) errors.append(np.max(np.abs(values - reference))) right.loglog(steps, errors, color=PURPLE, marker="o", ms=7, lw=2.8, label="ABM–PECE") fitted_slope, fitted_intercept = np.polyfit(np.log(steps), np.log(errors), 1) right.loglog(steps, np.exp(fitted_intercept) * steps**fitted_slope, color=YELLOW, lw=2.0, ls="--", label=fr"$p\approx{fitted_slope:.2f}$") left.set(xlabel="$t$", ylabel="$x(t)$") right.set(xlabel="Paso $h$", ylabel="Error máximo") right.invert_xaxis() right.set_xticks(steps) right.set_xticklabels(["0.25", "0.125", "0.0625", "0.03125"]) right.xaxis.set_minor_formatter(mpl.ticker.NullFormatter()) for ax in (left, right): ax.grid(True, which="both") ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) save(fig, "25_convergencia_abm.png") def plot_flow_family_comparison(lorenz: np.ndarray) -> None: rossler = rk4(rossler_field, np.array([0.1, 0.0, 0.0]), 0.018, 19000)[3500:] sprott = rk4(sprott_b_field, np.array([0.1, 0.1, 0.1]), 0.007, 42000)[4500:] samples = [lorenz[::6], rossler[::3], sprott[::5]] colors = [CYAN, MAGENTA, ORANGE] views = [(24, -58), (22, -54), (22, -48)] fig = plt.figure(figsize=(13.0, 4.6)) for index, (sample, color, view) in enumerate(zip(samples, colors, views), start=1): ax = fig.add_subplot(1, 3, index, projection="3d") ax.plot(sample[:, 0], sample[:, 1], sample[:, 2], color=color, lw=0.64, alpha=0.84) ax.set(xlabel="$x$", ylabel="$y$", zlabel="$z$") ax.view_init(*view) ax.grid(False) style_3d_panes(ax) save(fig, "26_comparacion_lorenz_rossler_sprott.png") def animate_logistic() -> None: fig, ax = plt.subplots(figsize=(7.4, 5.3)) values = np.linspace(2.8, 3.98, 38) def draw(frame: int): ax.clear() r = float(values[frame]) logistic_cobweb(ax, r, 0.183, 58) return [] animation = FuncAnimation(fig, draw, frames=len(values), interval=110, blit=False) animation.save(OUT / "logistic_cobweb_animation.gif", writer=PillowWriter(fps=9), dpi=92) plt.close(fig) def animate_lorenz_rotation(trajectory: np.ndarray) -> None: fig = plt.figure(figsize=(7.2, 5.4)) ax = fig.add_subplot(111, projection="3d") sample = trajectory[::4] ax.plot(sample[:, 0], sample[:, 1], sample[:, 2], color=CYAN, lw=0.45, alpha=0.72) ax.set(xlabel="$x$", ylabel="$y$", zlabel="$z$") ax.grid(False) style_3d_panes(ax) def draw(frame: int): ax.view_init(22, -70 + frame * 360 / 42) return [] animation = FuncAnimation(fig, draw, frames=42, interval=100, blit=False) animation.save(OUT / "lorenz_rotation.gif", writer=PillowWriter(fps=10), dpi=86) plt.close(fig) def animate_caputo_kernel() -> None: fig, ax = plt.subplots(figsize=(7.4, 5.2)) lag = np.geomspace(0.025, 12, 500) values = np.linspace(0.22, 0.96, 38) def draw(frame: int): ax.clear() q = float(values[frame]) kernel = lag ** (q - 1) / gamma(q) ax.loglog(lag, kernel, color=PURPLE, lw=2.8) ax.fill_between(lag, kernel, np.full_like(lag, kernel.min()), color=PURPLE, alpha=0.18) ax.set(xlabel=r"Antigüedad $t-s$", ylabel="Peso de memoria") ax.grid(True, which="both") return [] animation = FuncAnimation(fig, draw, frames=len(values), interval=105, blit=False) animation.save(OUT / "caputo_memory_animation.gif", writer=PillowWriter(fps=9), dpi=92) plt.close(fig) def animate_pendulum_phase() -> None: initial = [(-2.6, 0.2), (-1.1, 2.5), (1.9, -0.3), (4.4, -1.7)] trajectories = [ rk4( lambda state: np.array([state[1], -0.35 * state[1] - np.sin(state[0])]), np.array(state, dtype=float), 0.025, 1100, ) for state in initial ] colors = [MAGENTA, ORANGE, CYAN, YELLOW] fig, ax = plt.subplots(figsize=(7.6, 5.2)) x = np.linspace(-2.8 * np.pi, 2.8 * np.pi, 29) v = np.linspace(-3.2, 3.2, 23) xx, vv = np.meshgrid(x, v) ax.streamplot(xx, vv, vv, -0.35 * vv - np.sin(xx), color=GRID, density=1.0, linewidth=0.55, arrowsize=0.65) trails = [ax.plot([], [], color=color, lw=1.8)[0] for color in colors] dots = [ax.plot([], [], "o", color=color, ms=6)[0] for color in colors] ax.set(xlim=(x.min(), x.max()), ylim=(v.min(), v.max()), xlabel=r"$\theta$", ylabel=r"$\dot\theta$") def draw(frame: int): end = min(1 + frame * 13, len(trajectories[0])) start = max(0, end - 260) for trajectory, trail, dot in zip(trajectories, trails, dots): trail.set_data(trajectory[start:end, 0], trajectory[start:end, 1]) dot.set_data([trajectory[end - 1, 0]], [trajectory[end - 1, 1]]) return [*trails, *dots] animation = FuncAnimation(fig, draw, frames=80, interval=80, blit=True) animation.save(OUT / "pendulum_phase_animation.gif", writer=PillowWriter(fps=12), dpi=88) plt.close(fig) def animate_saddle_manifolds() -> None: angles = np.linspace(0, 2 * np.pi, 30, endpoint=False) points = 0.38 * np.column_stack([np.cos(angles), np.sin(angles)]) fig, ax = plt.subplots(figsize=(6.3, 5.3)) ax.axhline(0, color=ORANGE, lw=2.5, label=r"$W^u$") ax.axvline(0, color=CYAN, lw=2.5, label=r"$W^s$") scatter = ax.scatter(points[:, 0], points[:, 1], c=angles, cmap="turbo", s=32) ax.scatter([0], [0], color=YELLOW, s=65, edgecolor=BG, zorder=5) ax.set(xlim=(-3.2, 3.2), ylim=(-3.2, 3.2), xlabel="$x$", ylabel="$y$") ax.grid(True) ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) def draw(frame: int): time = frame * 0.045 evolved = np.column_stack([points[:, 0] * np.exp(time), points[:, 1] * np.exp(-time)]) scatter.set_offsets(evolved) return [scatter] animation = FuncAnimation(fig, draw, frames=46, interval=95, blit=True) animation.save(OUT / "saddle_manifolds_animation.gif", writer=PillowWriter(fps=10), dpi=92) plt.close(fig) def animate_lorenz_sensitivity() -> None: dt = 0.006 first = rk4(lorenz_field, np.array([1.0, 1.0, 1.0]), dt, 17000)[1800:] second = rk4(lorenz_field, np.array([1.00001, 1.0, 1.0]), dt, 17000)[1800:] fig = plt.figure(figsize=(7.4, 5.6)) ax = fig.add_subplot(111, projection="3d") ax.plot(first[::7, 0], first[::7, 1], first[::7, 2], color=GRID, lw=0.35, alpha=0.32) trail_a = ax.plot([], [], [], color=CYAN, lw=1.35, label="órbita A")[0] trail_b = ax.plot([], [], [], color=MAGENTA, lw=1.35, label="órbita B")[0] dot_a = ax.plot([], [], [], "o", color=CYAN, ms=5)[0] dot_b = ax.plot([], [], [], "o", color=MAGENTA, ms=5)[0] ax.set(xlim=(-22, 22), ylim=(-30, 30), zlim=(0, 52), xlabel="$x$", ylabel="$y$", zlabel="$z$") ax.view_init(23, -58) ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) style_3d_panes(ax) def draw(frame: int): end = min(1 + frame * 235, len(first)) start = max(0, end - 850) trail_a.set_data(first[start:end, 0], first[start:end, 1]) trail_a.set_3d_properties(first[start:end, 2]) trail_b.set_data(second[start:end, 0], second[start:end, 1]) trail_b.set_3d_properties(second[start:end, 2]) dot_a.set_data([first[end - 1, 0]], [first[end - 1, 1]]) dot_a.set_3d_properties([first[end - 1, 2]]) dot_b.set_data([second[end - 1, 0]], [second[end - 1, 1]]) dot_b.set_3d_properties([second[end - 1, 2]]) return [trail_a, trail_b, dot_a, dot_b] animation = FuncAnimation(fig, draw, frames=60, interval=95, blit=False) animation.save(OUT / "lorenz_sensitivity_animation.gif", writer=PillowWriter(fps=10), dpi=84) plt.close(fig) def poincare_crossings(trajectory: np.ndarray, plane: float = 27.0) -> np.ndarray: z0 = trajectory[:-1, 2] z1 = trajectory[1:, 2] mask = (z0 < plane) & (z1 >= plane) alpha = (plane - z0[mask]) / (z1[mask] - z0[mask]) return trajectory[:-1][mask] + alpha[:, None] * (trajectory[1:][mask] - trajectory[:-1][mask]) def animate_poincare_crossings(trajectory: np.ndarray) -> None: crossings = poincare_crossings(trajectory) fig, ax = plt.subplots(figsize=(6.8, 5.4)) scatter = ax.scatter([], [], s=15, c=[], cmap="turbo", vmin=0, vmax=max(1, len(crossings))) ax.set(xlim=(-18, 18), ylim=(-25, 25), xlabel="$x$ al cruzar", ylabel="$y$ al cruzar") ax.grid(True) def draw(frame: int): count = max(1, int((frame + 1) * len(crossings) / 55)) scatter.set_offsets(crossings[:count, :2]) scatter.set_array(np.arange(count)) return [scatter] animation = FuncAnimation(fig, draw, frames=55, interval=90, blit=True) animation.save(OUT / "poincare_crossings_animation.gif", writer=PillowWriter(fps=11), dpi=92) plt.close(fig) def animate_spectrum_window(trajectory: np.ndarray, dt: float) -> None: signal = trajectory[:, 0] - np.mean(trajectory[:, 0]) window_size = 2200 starts = np.linspace(0, len(signal) - window_size - 1, 44, dtype=int) time = np.arange(len(signal)) * dt fig, (top, bottom) = plt.subplots(2, 1, figsize=(7.6, 6.0)) top.plot(time, signal, color=GRID, lw=0.45) top.set(xlim=(time[0], time[-1]), ylim=(signal.min() * 1.08, signal.max() * 1.08), ylabel=r"$x(t)-\bar x$") span = top.axvspan(0, window_size * dt, color=ORANGE, alpha=0.25) spectral_line = bottom.semilogy([], [], color=ORANGE, lw=1.4)[0] bottom.set(xlim=(0.01, 8), ylim=(1e-5, 1e3), xlabel="Frecuencia", ylabel="PSD") bottom.grid(True, which="both") def draw(frame: int): start = starts[frame] segment = signal[start : start + window_size] frequency, density = welch(segment, fs=1 / dt, window="hann", nperseg=1024, noverlap=512) spectral_line.set_data(frequency[1:], density[1:]) left = time[start] right = time[start + window_size - 1] span.set_x(left) span.set_width(right - left) return [spectral_line, span] animation = FuncAnimation(fig, draw, frames=len(starts), interval=115, blit=False) animation.save(OUT / "spectrum_window_animation.gif", writer=PillowWriter(fps=9), dpi=88) plt.close(fig) def animate_sprott_orbit() -> None: trajectory = rk4(sprott_b_field, np.array([0.1, 0.1, 0.1]), 0.006, 65000)[5000:] fig = plt.figure(figsize=(7.2, 5.5)) ax = fig.add_subplot(111, projection="3d") line = ax.plot([], [], [], color=ORANGE, lw=1.05)[0] dot = ax.plot([], [], [], "o", color=YELLOW, ms=5)[0] ax.set(xlim=(trajectory[:, 0].min(), trajectory[:, 0].max()), ylim=(trajectory[:, 1].min(), trajectory[:, 1].max()), zlim=(trajectory[:, 2].min(), trajectory[:, 2].max())) ax.set(xlabel="$x$", ylabel="$y$", zlabel="$z$") ax.view_init(22, -48) ax.grid(False) style_3d_panes(ax) def draw(frame: int): end = min(1600 + frame * 880, len(trajectory)) sample = trajectory[:end:4] line.set_data(sample[:, 0], sample[:, 1]) line.set_3d_properties(sample[:, 2]) dot.set_data([trajectory[end - 1, 0]], [trajectory[end - 1, 1]]) dot.set_3d_properties([trajectory[end - 1, 2]]) return [line, dot] animation = FuncAnimation(fig, draw, frames=62, interval=85, blit=False) animation.save(OUT / "sprott_orbit_animation.gif", writer=PillowWriter(fps=11), dpi=86) plt.close(fig) def mittag_leffler_relaxation(q: float, time: np.ndarray) -> np.ndarray: result = np.ones_like(time, dtype=float) power = np.ones_like(time, dtype=float) argument = -(time ** q) for index in range(1, 180): power *= argument term = power / gamma(q * index + 1) result += term if np.max(np.abs(term)) < 2e-13: break return result def animate_fractional_relaxation() -> None: time = np.linspace(0, 3.0, 440) values = np.linspace(0.3, 1.0, 43) fig, ax = plt.subplots(figsize=(7.4, 5.2)) classical = np.exp(-time) ax.plot(time, classical, color=MUTED, lw=1.5, ls="--", label=r"$e^{-t}$") line = ax.plot([], [], color=CYAN, lw=2.7, label=r"$E_q(-t^q)$")[0] ax.set(xlim=(0, 3), ylim=(0, 1.04), xlabel="Tiempo", ylabel="Respuesta normalizada") ax.grid(True) ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) def draw(frame: int): q = float(values[frame]) line.set_data(time, mittag_leffler_relaxation(q, time)) return [line] animation = FuncAnimation(fig, draw, frames=len(values), interval=110, blit=True) animation.save(OUT / "fractional_relaxation_animation.gif", writer=PillowWriter(fps=9), dpi=92) plt.close(fig) def animate_abm_memory_weights() -> None: q = 0.65 history = 46 fig, ax = plt.subplots(figsize=(7.5, 5.2)) positions = np.arange(history) bars = ax.bar(positions, np.zeros(history), color=BLUE, width=0.85) ax.set(xlim=(-1, history), ylim=(0, 1.05), xlabel="Índice histórico $j$", ylabel="Peso normalizado") ax.grid(True, axis="y") current = ax.axvline(0, color=YELLOW, lw=2) def draw(frame: int): n = frame + 1 ages = np.arange(n, 0, -1, dtype=float) weights = ages**q - (ages - 1) ** q weights /= weights.max() heights = np.zeros(history) heights[:n] = weights for bar, height, age in zip(bars, heights, np.arange(history)): bar.set_height(height) bar.set_color(CYAN if age == n - 1 else BLUE) bar.set_alpha(0.92 if age < n else 0.1) current.set_xdata([n - 1, n - 1]) return [*bars, current] animation = FuncAnimation(fig, draw, frames=history, interval=105, blit=False) animation.save(OUT / "abm_memory_weights_animation.gif", writer=PillowWriter(fps=9), dpi=92) plt.close(fig) def animate_horseshoe_mechanism() -> None: gx, gy = np.meshgrid(np.linspace(-1, 1, 35), np.linspace(-0.42, 0.42, 16)) original = np.column_stack([gx.ravel(), gy.ravel()]) colors = np.repeat(np.linspace(0, 1, gx.shape[1]), gx.shape[0]).reshape(gx.shape[1], gx.shape[0]).T.ravel() fig, ax = plt.subplots(figsize=(6.8, 5.4)) scatter = ax.scatter(original[:, 0], original[:, 1], c=colors, cmap="turbo", s=10) ax.set(xlim=(-1.6, 1.6), ylim=(-1.35, 1.35), aspect="equal", xlabel="$x$", ylabel="$y$") ax.grid(True) def transformed(progress: float) -> np.ndarray: if progress < 1 / 3: u = progress * 3 target = np.column_stack([1.45 * original[:, 0], 0.45 * original[:, 1]]) return (1 - u) * original + u * target stretched = np.column_stack([1.45 * original[:, 0], 0.45 * original[:, 1]]) u = (progress - 1 / 3) * 1.5 angle = np.pi * (stretched[:, 0] / 2.9 + 0.5) folded = np.column_stack([0.95 * np.cos(angle), 1.05 * np.sin(angle) + stretched[:, 1]]) if progress < 2 / 3: return (1 - u) * stretched + u * folded u = (progress - 2 / 3) * 3 shifted = folded + np.column_stack([np.zeros(len(folded)), -0.18 * np.ones(len(folded))]) return (1 - u) * folded + u * shifted def draw(frame: int): points = transformed(frame / 59) scatter.set_offsets(points) return [scatter] animation = FuncAnimation(fig, draw, frames=60, interval=95, blit=True) animation.save(OUT / "horseshoe_mechanism_animation.gif", writer=PillowWriter(fps=10), dpi=92) plt.close(fig) def animate_pitchfork_parameter() -> None: negative = np.linspace(-1.1, 0.0, 260) positive = np.linspace(0.0, 1.2, 320) parameters = np.linspace(-0.85, 1.05, 44) x = np.linspace(-1.45, 1.45, 700) fig, (left, right) = plt.subplots(1, 2, figsize=(8.4, 4.3)) left.plot(negative, np.zeros_like(negative), color=CYAN, lw=2.8) left.plot(positive, np.zeros_like(positive), color=MAGENTA, lw=2.2, ls="--") left.plot(positive, np.sqrt(positive), color=CYAN, lw=2.8) left.plot(positive, -np.sqrt(positive), color=CYAN, lw=2.8) left.set(xlabel=r"$\mu$", ylabel=r"$x^*$", xlim=(-1.1, 1.2), ylim=(-1.25, 1.25)) left.grid(True) parameter_line = left.axvline(parameters[0], color=YELLOW, lw=2.4) equilibria = left.scatter([], [], color=YELLOW, edgecolor=BG, s=55, zorder=6) potential_line = right.plot([], [], color=PURPLE, lw=3.0)[0] critical_points = right.scatter([], [], color=YELLOW, edgecolor=BG, s=55, zorder=6) right.set(xlabel="$x$", ylabel=r"$V_\mu(x)$", xlim=(-1.45, 1.45), ylim=(-0.45, 1.45)) right.grid(True) def draw(frame: int): mu = float(parameters[frame]) parameter_line.set_xdata([mu, mu]) if mu > 0.0: roots = np.array([-np.sqrt(mu), 0.0, np.sqrt(mu)]) else: roots = np.array([0.0]) equilibria.set_offsets(np.column_stack([np.full(len(roots), mu), roots])) potential = 0.25 * x**4 - 0.5 * mu * x**2 potential_line.set_data(x, potential) critical_points.set_offsets(np.column_stack([roots, 0.25 * roots**4 - 0.5 * mu * roots**2])) return [parameter_line, equilibria, potential_line, critical_points] animation = FuncAnimation(fig, draw, frames=len(parameters), interval=105, blit=False) animation.save(OUT / "pitchfork_parameter_animation.gif", writer=PillowWriter(fps=9), dpi=88) plt.close(fig) def animate_nullclines_flow() -> None: initial = [(-1.65, 1.15), (-0.7, -1.25), (0.25, 1.45), (1.4, -1.25), (0.6, 0.25)] trajectories = [rk4(double_well_field, np.array(state, dtype=float), 0.024, 1500) for state in initial] colors = [MAGENTA, ORANGE, CYAN, YELLOW, PURPLE] x = np.linspace(-1.9, 1.9, 29) velocity = np.linspace(-1.8, 1.8, 27) xx, vv = np.meshgrid(x, velocity) fig, ax = plt.subplots(figsize=(7.2, 5.2)) ax.streamplot(xx, vv, vv, xx - xx**3 - 0.45 * vv, color=GRID, density=1.05, linewidth=0.65, arrowsize=0.72) curve_x = np.linspace(-1.75, 1.75, 700) ax.axhline(0.0, color=CYAN, lw=2.6) ax.plot(curve_x, (curve_x - curve_x**3) / 0.45, color=MAGENTA, lw=2.6) trails = [ax.plot([], [], color=color, lw=2.0)[0] for color in colors] dots = [ax.plot([], [], "o", color=color, ms=5.5)[0] for color in colors] ax.set(xlabel="$x$", ylabel="$v$", xlim=(-1.9, 1.9), ylim=(-1.8, 1.8)) def draw(frame: int): end = min(1 + frame * 24, len(trajectories[0])) start = max(0, end - 320) for trajectory, trail, dot in zip(trajectories, trails, dots): trail.set_data(trajectory[start:end, 0], trajectory[start:end, 1]) dot.set_data([trajectory[end - 1, 0]], [trajectory[end - 1, 1]]) return [*trails, *dots] animation = FuncAnimation(fig, draw, frames=62, interval=85, blit=True) animation.save(OUT / "nullclines_flow_animation.gif", writer=PillowWriter(fps=11), dpi=88) plt.close(fig) def animate_integrator_step_refinement() -> None: steps = [0.5, 0.4, 0.25, 0.2, 0.125, 0.1, 0.0625, 0.05] dense_time = np.linspace(0.0, 2.0, 900) fig, ax = plt.subplots(figsize=(7.2, 5.2)) ax.plot(dense_time, np.exp(-2.0 * dense_time), color=TEXT, lw=3.0, label=r"$e^{-2t}$") lines = { "euler": ax.plot([], [], color=MAGENTA, marker="o", ms=4, lw=2.0, label="Euler")[0], "heun": ax.plot([], [], color=ORANGE, marker="s", ms=4, lw=2.0, label="Heun")[0], "rk4": ax.plot([], [], color=CYAN, marker="^", ms=4, lw=2.0, label="RK4")[0], } ax.set(xlabel="$t$", ylabel="$y(t)$", xlim=(0.0, 2.0), ylim=(-0.08, 1.06)) ax.grid(True) ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) def draw(frame: int): step = steps[min(frame // 4, len(steps) - 1)] for method, line in lines.items(): time, values = integrate_decay(method, step) line.set_data(time, values) return list(lines.values()) animation = FuncAnimation(fig, draw, frames=len(steps) * 4, interval=125, blit=True) animation.save(OUT / "integrator_step_refinement_animation.gif", writer=PillowWriter(fps=8), dpi=90) plt.close(fig) def animate_hopf_cycle_birth() -> None: parameters = np.linspace(-0.34, 0.58, 44) negative = np.linspace(-0.5, 0.0, 220) positive = np.linspace(0.0, 0.7, 260) trajectories = [rk4(hopf_field(float(mu)), np.array([0.86, 0.13]), 0.025, 1100) for mu in parameters] fig, (left, right) = plt.subplots(1, 2, figsize=(8.4, 4.3)) left.plot(negative, np.zeros_like(negative), color=CYAN, lw=2.7) left.plot(positive, np.zeros_like(positive), color=MAGENTA, lw=2.1, ls="--") left.plot(positive, np.sqrt(positive), color=CYAN, lw=2.7) left.set(xlabel=r"$\mu$", ylabel="$r$", xlim=(-0.5, 0.7), ylim=(-0.03, 0.9)) left.grid(True) parameter_line = left.axvline(parameters[0], color=YELLOW, lw=2.3) branch_dot = left.scatter([], [], color=YELLOW, edgecolor=BG, s=52, zorder=5) trail = right.plot([], [], color=ORANGE, lw=2.2)[0] cycle = right.plot([], [], color=CYAN, lw=2.4, ls="--")[0] right.scatter([0.0], [0.0], color=YELLOW, edgecolor=BG, s=48, zorder=5) right.set(xlabel="$x$", ylabel="$y$", xlim=(-1.05, 1.05), ylim=(-1.05, 1.05), aspect="equal") right.grid(True) angle = np.linspace(0.0, 2.0 * np.pi, 500) def draw(frame: int): mu = float(parameters[frame]) trajectory = trajectories[frame] parameter_line.set_xdata([mu, mu]) branch_dot.set_offsets([[mu, np.sqrt(mu) if mu > 0.0 else 0.0]]) trail.set_data(trajectory[:, 0], trajectory[:, 1]) if mu > 0.0: radius = np.sqrt(mu) cycle.set_data(radius * np.cos(angle), radius * np.sin(angle)) else: cycle.set_data([], []) return [parameter_line, branch_dot, trail, cycle] animation = FuncAnimation(fig, draw, frames=len(parameters), interval=105, blit=False) animation.save(OUT / "hopf_cycle_birth_animation.gif", writer=PillowWriter(fps=9), dpi=88) plt.close(fig) def animate_henon_iteration() -> None: grid_x, grid_y = np.meshgrid(np.linspace(-1.05, 1.05, 52), np.linspace(-0.28, 0.28, 24)) points = np.column_stack([grid_x.ravel(), grid_y.ravel()]) colors = np.tile(np.linspace(0.0, 1.0, grid_x.shape[1]), grid_x.shape[0]) stages = [points] for _ in range(8): next_points = henon_step(stages[-1]) invalid = (np.abs(next_points[:, 0]) > 2.2) | (np.abs(next_points[:, 1]) > 0.75) next_points[invalid] = np.nan stages.append(next_points) fig, ax = plt.subplots(figsize=(7.0, 5.1)) scatter = ax.scatter(points[:, 0], points[:, 1], c=colors, cmap="turbo", s=8, alpha=0.84) ax.set(xlabel="$x_n$", ylabel="$y_n$", xlim=(-1.8, 1.8), ylim=(-0.62, 0.62)) ax.grid(True) def draw(frame: int): scatter.set_offsets(stages[min(frame // 4, len(stages) - 1)]) return [scatter] animation = FuncAnimation(fig, draw, frames=len(stages) * 4, interval=115, blit=True) animation.save(OUT / "henon_iteration_animation.gif", writer=PillowWriter(fps=9), dpi=90) plt.close(fig) def animate_basin_refinement() -> None: resolutions = [(24, 22), (34, 30), (48, 42), (66, 58), (90, 78), (122, 104), (158, 134)] basin_images = [] for nx, ny in resolutions: x_values = np.linspace(-2.1, 2.1, nx) velocity_values = np.linspace(-2.0, 2.0, ny) basin_images.append(classify_double_well_basins(x_values, velocity_values, steps=560)) cmap = mpl.colors.ListedColormap([BLUE, MAGENTA]) fig, ax = plt.subplots(figsize=(7.0, 5.2)) image = ax.imshow(basin_images[0], origin="lower", extent=(-2.1, 2.1, -2.0, 2.0), cmap=cmap, interpolation="nearest", aspect="auto") ax.scatter([-1.0, 0.0, 1.0], [0.0, 0.0, 0.0], color=[BLUE, YELLOW, MAGENTA], edgecolor=BG, s=55, zorder=5) ax.set(xlabel="$x_0$", ylabel="$v_0$") def draw(frame: int): image.set_data(basin_images[min(frame // 4, len(basin_images) - 1)]) return [image] animation = FuncAnimation(fig, draw, frames=len(basin_images) * 4, interval=135, blit=True) animation.save(OUT / "basin_refinement_animation.gif", writer=PillowWriter(fps=7), dpi=90) plt.close(fig) def animate_return_cobweb() -> None: radius = np.linspace(0.0, 0.75, 700) returned = hopf_return_map(radius) orbit = [0.69] for _ in range(16): orbit.append(float(hopf_return_map(np.array([orbit[-1]]))[0])) cobweb_x = [orbit[0]] cobweb_y = [0.0] for current, following in zip(orbit[:-1], orbit[1:]): cobweb_x.extend([current, following]) cobweb_y.extend([following, following]) fig, ax = plt.subplots(figsize=(6.6, 5.2)) ax.plot(radius, returned, color=CYAN, lw=2.8, label=r"$P(r)$") ax.plot(radius, radius, color=MUTED, lw=1.7, ls="--", label=r"$r$") line = ax.plot([], [], color=MAGENTA, lw=2.0)[0] dot = ax.plot([], [], "o", color=YELLOW, ms=6)[0] ax.set(xlabel="$r_n$", ylabel="$r_{n+1}$", xlim=(0.0, 0.75), ylim=(0.0, 0.75)) ax.grid(True) ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) def draw(frame: int): count = min(2 + frame, len(cobweb_x)) line.set_data(cobweb_x[:count], cobweb_y[:count]) dot.set_data([cobweb_x[count - 1]], [cobweb_y[count - 1]]) return [line, dot] animation = FuncAnimation(fig, draw, frames=len(cobweb_x) - 1, interval=125, blit=True) animation.save(OUT / "return_cobweb_animation.gif", writer=PillowWriter(fps=8), dpi=92) plt.close(fig) def animate_topological_degree() -> None: angle = np.linspace(0.0, 2.0 * np.pi, 64) domain = np.column_stack([np.cos(angle), np.sin(angle)]) image_points = np.column_stack([ np.cos(2.0 * angle) + 0.35 * np.cos(angle), np.sin(2.0 * angle) + 0.35 * np.sin(angle), ]) fig, (left, right) = plt.subplots(1, 2, figsize=(8.2, 4.2)) left.plot(domain[:, 0], domain[:, 1], color=GRID, lw=1.3) right.plot(image_points[:, 0], image_points[:, 1], color=GRID, lw=1.3) domain_trail = left.plot([], [], color=CYAN, lw=3.0)[0] image_trail = right.plot([], [], color=ORANGE, lw=3.0)[0] domain_dot = left.plot([], [], "o", color=YELLOW, ms=6)[0] image_dot = right.plot([], [], "o", color=YELLOW, ms=6)[0] for ax in (left, right): ax.scatter([0.0], [0.0], color=MAGENTA, s=44, edgecolor=BG, zorder=5) ax.set(xlabel="$x$", ylabel="$y$", xlim=(-1.55, 1.55), ylim=(-1.55, 1.55), aspect="equal") ax.grid(True) def draw(frame: int): count = frame + 1 domain_trail.set_data(domain[:count, 0], domain[:count, 1]) image_trail.set_data(image_points[:count, 0], image_points[:count, 1]) domain_dot.set_data([domain[count - 1, 0]], [domain[count - 1, 1]]) image_dot.set_data([image_points[count - 1, 0]], [image_points[count - 1, 1]]) return [domain_trail, image_trail, domain_dot, image_dot] animation = FuncAnimation(fig, draw, frames=len(angle), interval=95, blit=True) animation.save(OUT / "topological_degree_animation.gif", writer=PillowWriter(fps=10), dpi=88) plt.close(fig) def animate_recurrence_build(trajectory: np.ndarray) -> None: _, recurrence, _ = recurrence_data(trajectory, count=420) progressive = np.zeros_like(recurrence) cmap = mpl.colors.ListedColormap([BG, CYAN]) fig, ax = plt.subplots(figsize=(6.1, 5.4)) image = ax.imshow( progressive, origin="lower", cmap=cmap, vmin=0, vmax=1, interpolation="nearest", aspect="equal", ) ax.set(xlabel="$i$", ylabel="$j$") def draw(frame: int): count = max(1, int((frame + 1) * len(recurrence) / 45)) progressive[:count, :] = recurrence[:count, :] image.set_data(progressive.copy()) return [image] animation = FuncAnimation(fig, draw, frames=45, interval=95, blit=True) animation.save(OUT / "recurrence_build_animation.gif", writer=PillowWriter(fps=10), dpi=90) plt.close(fig) def animate_fractional_order_variation() -> None: time = np.geomspace(0.03, 4.0, 650) orders = np.linspace(0.18, 0.94, 42) fig, (left, right) = plt.subplots(1, 2, figsize=(8.4, 4.3)) rl_line = left.loglog([], [], color=MAGENTA, lw=2.8, label=r"${}^{RL}D^q1$")[0] left.plot(time, np.full_like(time, 1e-3), color=CYAN, lw=2.5, ls="--", label=r"${}^{C}D^q1=0$") power_line = right.loglog([], [], color=ORANGE, lw=2.8, label=r"$D^q t^2$")[0] right.loglog(time, time**2, color=MUTED, lw=1.5, ls="--", label=r"$t^2$") left.set(xlabel="$t$", ylabel=r"$D^q1$", xlim=(time.min(), time.max()), ylim=(7e-4, 7.0)) right.set(xlabel="$t$", ylabel=r"$D^q t^2$", xlim=(time.min(), time.max()), ylim=(0.002, 25.0)) for ax in (left, right): ax.grid(True, which="both") ax.legend(frameon=True, facecolor=BG, edgecolor=GRID) def draw(frame: int): q = float(orders[frame]) rl_line.set_data(time, time ** (-q) / gamma(1.0 - q)) power_line.set_data(time, gamma(3.0) / gamma(3.0 - q) * time ** (2.0 - q)) return [rl_line, power_line] animation = FuncAnimation(fig, draw, frames=len(orders), interval=105, blit=True) animation.save(OUT / "fractional_order_variation_animation.gif", writer=PillowWriter(fps=9), dpi=88) plt.close(fig) def animate_abm_convergence() -> None: q = 0.72 steps = [0.375, 0.3, 0.25, 0.1875, 0.15, 0.125, 0.1, 0.075, 0.0625, 0.05] dense_time = np.linspace(0.0, 3.0, 800) exact = mittag_leffler_relaxation(q, dense_time) solutions = [] errors = [] for step in steps: time, values = abm_pece_relaxation(q, step) solutions.append((time, values)) errors.append(np.max(np.abs(values - mittag_leffler_relaxation(q, time)))) fig, (left, right) = plt.subplots(1, 2, figsize=(8.5, 4.3)) left.plot(dense_time, exact, color=TEXT, lw=2.8, label=r"$E_q(-t^q)$") numerical = left.plot([], [], color=PURPLE, marker="o", ms=3, lw=2.0, label="ABM–PECE")[0] left.set(xlabel="$t$", ylabel="$x(t)$", xlim=(0.0, 3.0), ylim=(0.0, 1.04)) left.grid(True) left.legend(frameon=True, facecolor=BG, edgecolor=GRID) convergence = right.loglog([], [], color=ORANGE, marker="o", ms=6, lw=2.4)[0] right.set(xlabel="Paso $h$", ylabel="Error máximo", xlim=(0.42, 0.045), ylim=(min(errors) * 0.65, max(errors) * 1.5)) right.set_xticks([0.3, 0.1, 0.05]) right.set_xticklabels(["0.3", "0.1", "0.05"]) right.xaxis.set_minor_formatter(mpl.ticker.NullFormatter()) right.grid(True, which="both") def draw(frame: int): index = min(frame // 3, len(steps) - 1) time, values = solutions[index] numerical.set_data(time, values) convergence.set_data(steps[: index + 1], errors[: index + 1]) return [numerical, convergence] animation = FuncAnimation(fig, draw, frames=len(steps) * 3, interval=120, blit=True) animation.save(OUT / "abm_convergence_animation.gif", writer=PillowWriter(fps=8), dpi=88) plt.close(fig) def generate_current_theme(include_animations: bool) -> None: OUT.mkdir(parents=True, exist_ok=True) configure_style() plot_flow_phase_space() plot_saddle_manifolds() plot_logistic_cobweb() plot_logistic_bifurcation() lorenz, dt = lorenz_trajectory() plot_lorenz(lorenz) plot_poincare(lorenz) plot_lyapunov_divergence() plot_spectrum(lorenz, dt) plot_caputo_kernel() plot_sprott_b() plot_pendulum_potential_phase() plot_pitchfork_bifurcation() plot_trace_determinant_plane() plot_nullclines_field() plot_integrator_comparison() plot_step_convergence() plot_hopf_bifurcation() plot_henon_folding() plot_bistable_basin() plot_return_floquet() plot_topological_degree() plot_recurrence(lorenz, dt) plot_correlation_dimension(lorenz) plot_rl_caputo_powers() plot_abm_convergence() plot_flow_family_comparison(lorenz) if include_animations: animate_logistic() animate_lorenz_rotation(lorenz) animate_caputo_kernel() animate_pendulum_phase() animate_saddle_manifolds() animate_lorenz_sensitivity() animate_poincare_crossings(lorenz) animate_spectrum_window(lorenz, dt) animate_sprott_orbit() animate_fractional_relaxation() animate_abm_memory_weights() animate_horseshoe_mechanism() animate_pitchfork_parameter() animate_nullclines_flow() animate_integrator_step_refinement() animate_hopf_cycle_birth() animate_henon_iteration() animate_basin_refinement() animate_return_cobweb() animate_topological_degree() animate_recurrence_build(lorenz) animate_fractional_order_variation() animate_abm_convergence() print(f"Generated course assets in {OUT}") def main() -> None: select_theme("dark") generate_current_theme(include_animations=True) select_theme("light") generate_current_theme(include_animations=False) if __name__ == "__main__": main()