:orphan:

Examples
========

Runnable demos of lyapax's Lyapunov-exponent engine, in
`sphinx-gallery <https://sphinx-gallery.github.io/>`_ format (each file is
plain, directly-runnable Python, and also renders into the docs gallery --
see ``docs/``).

Demos cover: inspecting a model's raw simulated time series as a sanity
check before computing exponents (single systems and coupled networks
alike), benchmark systems with a known answer (Tiers 0-2 of
:doc:`/background/validation`), coupled networks (Tier 3), the extensibility
of custom coupling functions, a speed/accuracy characterization of the
current implementation (including RK4 vs RK6 convergence-order slopes),
delayed (DDE) networks -- both a per-edge linear
delay sweep against a closed-form (Lambert W) reference and a delayed
Kuramoto network showing what transmission delay does to synchronization
(Tier 4-5) -- a matrix-free (jvp/vmap) tangent-propagation speedup on
a large network, a jax.vmap parameter-sweep helper (one batched call
reproducing an earlier Python-loop G-sweep, faster and bit-for-bit
identical), the public front door (``ode_problem``, ``Network``/
``network_problem``, ``dde_problem``/``network_dde_problem``) that gives
plain, coupled, and delayed systems the same dynamics/network/coupling/
integrator/problem construction recipe, with ``state0``/``dt`` given once,
and grid-snapped vs. Hermite-interpolated DDE history reads
(``interpolate=True``) compared against a closed-form (Lambert W)
reference, showing why the latter converges smoothly as ``dt`` shrinks
and the former does not -- and a CPU-vs-GPU wall-time comparison across
growing network size, showing that ``lyapax``'s GPU support only pays off
once the per-step arithmetic is large enough to amortize a GPU's fixed
kernel-launch/transfer overhead, with the crossover point measured rather
than asserted -- and an adaptive-step ODE integrator
(``lyapax.adaptive.diffrax_adaptive_step``, backed by diffrax) that decouples
the renormalization sampling interval from the integrator's own internal
step size, with a tolerance-convergence sweep, a cross-check against
fixed-step rk4, and a demonstration that differentiating a Lyapunov exponent
through it requires ``jax.jacfwd`` rather than ``jax.grad`` -- and a
convergence diagnostic (``lyapax.core.convergence_drift``) that summarizes
how much a run's running exponent estimate moved over the tail of the run,
paired with ``result.checkpoint``/``lyapunov_spectrum(..., resume=...)``,
which continues a fixed-``n_steps`` run from where it left off instead of
restarting, letting a caller run in inspectable chunks and stop once
``convergence_drift`` says the estimate has settled -- and the same
run-inspect-resume loop for a DDE (``lyapax.dde.lyapunov_spectrum_dde(...,
resume=...)`` / ``lyapax.dde.DDECheckpoint``) on the Mackey-Glass chaotic
benchmark, showing that a DDE checkpoint must additionally carry the delay
ring buffer's state (not just trajectory state and tangent basis) for a
resumed run to continue correctly -- and differentiating a Lyapunov
exponent w.r.t. a system parameter with ``jax.grad``/``jax.jacfwd``,
demonstrating gradient-based tuning of a parameter toward a target
exponent on a non-chaotic system, then measuring how the same gradient
becomes numerically meaningless (grows by many orders of magnitude with
trajectory length rather than converging) once the underlying system is
genuinely chaotic -- and the Kaplan-Yorke (Lyapunov) dimension
(``lyapax.core.kaplan_yorke_dimension``), a pure post-processing readout
of an existing Lyapunov spectrum, computed for Lorenz and Rossler and
matching published values, with the cumulative-sum-crossing-zero mechanics
behind the formula plotted directly and the partial-spectrum guard
(``d_total=``) shown raising instead of silently understating the
dimension.

New capability -> new demo: as engine features land (a new coupling kind, a
new delay structure, a performance change, ...), add a runnable example
that exercises it, not just unit-test coverage.


.. raw:: html

  <div id='sg-tag-list' class='sphx-glr-tag-list'></div>


.. raw:: html

    <div class="sphx-glr-thumbnails">

.. thumbnail-parent-div-open

.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Before treating a Lyapunov-exponent estimate as meaningful, it&#x27;s worth first checking that the model behind step_fn is doing what you expect -- settling onto the attractor you had in mind, not diverging, not stuck at a fixed point. step_fn (however it&#x27;s built -- ode_problem, network_problem, ...) is a plain state -&gt; new_state map under the hood, so the exact function handed to lyapunov_spectrum can also just be run forward on its own and inspected as an ordinary time series -- no separate simulation code path to maintain. lyapax.utils.simulate_trajectory( step_fn, state0, n_steps, dt=dt) does this: it iterates step_fn via jax.lax.scan and returns (t, traj), the matching time axis and the whole trajectory (shape (n_steps + 1, d), including the initial state). Both problem objects expose their built step_fn/``state0`` as plain attributes, so simulate_trajectory and lyapunov_spectrum can share the exact same problem without re-deriving anything.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_00_time_series_sanity_check_thumb.png
    :alt:

  :doc:`/auto_examples/00_time_series_sanity_check`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Sanity check: watch the time series before trusting the exponents</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="The simplest possible correctness check for the Benettin/QR engine: for a linear system x&#x27; = A x, the Lyapunov spectrum is exactly the real parts of the eigenvalues of A -- no chaos, no literature value to trust, just linear algebra. See Tier 0.1 &lt;validation-tier-0-1&gt; of the validation guide.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_01_linear_ode_thumb.png
    :alt:

  :doc:`/auto_examples/01_linear_ode`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Linear ODE: exact Lyapunov spectrum</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="rhs is just a plain state -&gt; jnp.ndarray callable -- not a fixed set of named systems dispatched through hardcoded if/elif branches internally. That means using your own equations never requires touching the library&#x27;s source or registering anything: any function with the right signature works with ode_problem, exactly like the built-in lorenz/``rossler``/``linear_system`` in lyapax.systems.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_01b_custom_system_thumb.png
    :alt:

  :doc:`/auto_examples/01b_custom_system`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Writing your own system: no registry needed</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Discrete maps avoid integration-scheme error entirely -- useful for isolating bugs in the QR/renormalization bookkeeping itself, independent of any ODE-integrator accuracy question. See the validation guide, Tier 0.2 &lt;validation-tier-0-2&gt; and Tier 0.3 &lt;validation-tier-0-3&gt;.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_02_chaotic_maps_thumb.png
    :alt:

  :doc:`/auto_examples/02_chaotic_maps`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Chaotic maps: exact analytic Lyapunov exponents</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Full-spectrum Lyapunov exponents for the two classic chaotic-ODE benchmarks, cross-checked against a structural invariant (constant or directly-computable phase-space divergence) as well as published values. See the validation guide, Tier 1 &lt;validation-tier-1&gt; and Tier 2 &lt;validation-tier-2&gt;.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_03_chaotic_flows_thumb.png
    :alt:

  :doc:`/auto_examples/03_chaotic_flows`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Chaotic flows: Lorenz and Rossler</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Isolates &quot;did we wire the coupling term into the Jacobian correctly&quot; from &quot;is chaos numerically well-resolved&quot;: for a linear coupled network, the full Lyapunov spectrum is exactly the eigenvalues of the constant network Jacobian gammaI + GW. See Tier 3 &lt;validation-tier-3&gt; of the validation guide.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_04_linear_network_thumb.png
    :alt:

  :doc:`/auto_examples/04_linear_network`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Coupled linear network: exact eigenvalues</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Sweeps the global coupling strength G and tracks the top two Lyapunov exponents of a small heterogeneous-frequency Kuramoto network. At G=0 both are exactly 0 (decoupled oscillators: dtheta/dt = omega is state-independent, so the step map&#x27;s Jacobian is exactly the identity -- see tests/test_network.py::test_kuramoto_zero_coupling_gives_exactly_zero_spectrum).">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_05_kuramoto_sync_thumb.png
    :alt:

  :doc:`/auto_examples/05_kuramoto_sync`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Kuramoto network: synchronization transition</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="lyapax coupling is a plain callable coupling_fn(cvar_state, weights, params) -&gt; coupling -- not a fixed set of named coupling &quot;kinds&quot; dispatched through hardcoded if/elif branches internally. That means adding a new coupling rule never requires touching the library&#x27;s source: any function with the right signature is a first-class coupling, exactly like the built-in linear_coupling/ sigmoidal_coupling/``kuramoto_coupling`` in lyapax.coupling. This demo writes one from scratch -- no import from lyapax.coupling at all -- and plugs it directly into network_problem, reproducing the exact linear-network result from sphx_glr_auto_examples_04_linear_network.py.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_06_custom_coupling_thumb.png
    :alt:

  :doc:`/auto_examples/06_custom_coupling`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Custom coupling functions: no registry needed</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Two practical questions before trusting lyapax on a new problem:">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_07_speed_and_accuracy_thumb.png
    :alt:

  :doc:`/auto_examples/07_speed_and_accuracy`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Speed and accuracy characterization</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="The simplest possible &quot;delayed coupling&quot; demo: two identical linear nodes coupled through each other&#x27;s delayed state, x1&#x27; = gammax1 + Gx2(t-tau), x2&#x27; = gammax2 + Gx1(t-tau). This sweeps the delay tau itself, to show what delayed coupling does to the spectrum -- how it differs from the zero-delay case in sphx_glr_auto_examples_04_linear_network.py.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_08_delayed_coupling_thumb.png
    :alt:

  :doc:`/auto_examples/08_delayed_coupling`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Delayed coupling: two-node linear DDE network vs delay</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Extends sphx_glr_auto_examples_05_kuramoto_sync.py&#x27;s zero-delay synchronization sweep to a delayed Kuramoto network -- the same 6-oscillator, heterogeneous- frequency, all-to-all system, but now each oscillator feels its neighbors&#x27; phases as they were tau time units ago rather than instantaneously, via dtheta_i/dt = omega_i + (G/N) sum_j sin(theta_j(t-tau) - theta_i(t)). Runs the same G sweep two ways -- zero delay (exactly sphx_glr_auto_examples_05_kuramoto_sync.py) and tau=0.3 (the delayed/DDE engine) -- and overlays lambda_2 from both, since lambda_2 is what tracks synchronization (see sphx_glr_auto_examples_05_kuramoto_sync.py&#x27;s docstring for why: lambda_1 is pinned to 0 by the model&#x27;s exact rotational symmetry and never signals anything).">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_09_kuramoto_delayed_network_thumb.png
    :alt:

  :doc:`/auto_examples/09_kuramoto_delayed_network`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Kuramoto network with delayed coupling: what transmission delay does</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Tracking only the top k Lyapunov exponents of a d-dimensional system (``k &lt; d``) should only cost O(k) work per step, not O(d) -- that&#x27;s the whole point of a partial spectrum. The naive way to linearize a step function, jax.jacfwd(step_fn)(state), doesn&#x27;t get this for free: it always computes the full d x d Jacobian and only afterwards multiplies by the k tracked tangent columns (``jac @ Y``), so the d - k untracked columns are wasted work. lyapax.core.lyapunov_spectrum instead computes exactly the k columns that are actually needed, via one jax.jvp (forward-mode directional derivative) per tracked column, batched together with jax.vmap -- the same &quot;matrix-free&quot; idea underlying lyapax.dde.lyapunov_spectrum_dde&#x27;s delayed-network engine, here applied to the plain (non-delayed) case. See tests/test_lyapunov_core.py::test_tangent_propagation_matches_dense_jacfwd for the correctness check (this script is purely about cost).">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_10_matrix_free_scaling_thumb.png
    :alt:

  :doc:`/auto_examples/10_matrix_free_scaling`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Matrix-free tangent propagation: dense jacfwd vs jvp/vmap</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="sphx_glr_auto_examples_05_kuramoto_sync.py swept the coupling strength G with a plain Python for loop, one full lyapunov_spectrum call per G value. lyapax.sweep.sweep_lyapunov_spectrum does the same sweep as a single jax.vmap-batched call instead -- one XLA program that computes every grid point together, rather than len(G_values) separate Python-level dispatches (each of which re-enters the JAX/XLA call machinery even when, as sphx_glr_auto_examples_07_speed_and_accuracy.py found, the compiled executable itself is cached across same-shape calls).">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_11_vmap_parameter_sweep_thumb.png
    :alt:

  :doc:`/auto_examples/11_vmap_parameter_sweep`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">vmap parameter sweeps: the same Kuramoto transition, one XLA call</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Every example so far reaches the Lyapunov engine through a different door: plain systems build a step with rk4_step and pass state0/``dt`` a second time into lyapunov_spectrum, networks assemble a step through make_network_step_fn with topology/coupling/integrator all passed as separate positional arguments (and, again, a second dt), and delayed systems hand lyapunov_spectrum_dde a raw (state0, buf0, params, dt) tuple that exposes the ring buffer directly. Passing dt twice is not just repetitive -- nothing checks that the two copies agree, so a typo in either one silently computes the spectrum for the wrong step size.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_12_public_api_overview_thumb.png
    :alt:

  :doc:`/auto_examples/12_public_api_overview`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">The lyapax front door: problem objects for ODEs, networks, and DDEs</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Every delayed-coupling example so far has read history off an exact ring- buffer grid point: the physical delay tau is rounded to the nearest whole number of dt steps before the run even starts (``resolve_tau_steps``), and every later lookup pulls that exact stored sample. That&#x27;s simple and fully differentiable, but it means the delay actually being simulated is tau_eff = tau_steps  dt, not the tau that was asked for -- and, because rounding to the nearest integer is not a smooth function of dt, the size* of that mismatch does not shrink smoothly as dt gets finer. It can even get worse at a smaller dt if the rounding happens to land less favorably there than at a coarser one.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_13_dde_history_interpolation_thumb.png
    :alt:

  :doc:`/auto_examples/13_dde_history_interpolation`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">DDE history reads: grid-snapped vs Hermite-interpolated</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="lyapax runs on GPU with no code changes -- JAX picks its backend (CPU, GPU, ...) from an environment variable at startup, and none of lyapax&#x27;s own functions know or care which one is active. But &quot;runs on GPU&quot; and &quot;faster on GPU&quot; are two different claims. A GPU wins on throughput: it does many floating-point operations in parallel, at the cost of a fixed per-call overhead (kernel launch, host&lt;-&gt;device data transfer) that a CPU doesn&#x27;t pay. For a small problem -- a state vector with only a handful of entries, say -- that fixed overhead dwarfs the actual arithmetic, and the CPU wins outright. Once the per-step arithmetic is large enough to amortize the overhead, the GPU pulls ahead, often by a large margin.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_14_gpu_acceleration_thumb.png
    :alt:

  :doc:`/auto_examples/14_gpu_acceleration`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">GPU acceleration: when does it actually pay off?</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="Every other ODE example uses a fixed-step integrator (``&quot;rk4&quot;``, &quot;rk6&quot;, ...): the same dt is both the sampling interval between QR renormalizations and the integrator&#x27;s own internal step size, so getting a smaller local truncation error means shrinking dt globally, which also means more (cheap) raw steps between renormalizations.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_15_adaptive_ode_thumb.png
    :alt:

  :doc:`/auto_examples/15_adaptive_ode`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Adaptive-step ODE integration (diffrax)</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="lyapunov_spectrum always runs a fixed n_steps -- there is no built-in stopping criterion, adaptive or otherwise (see sphx_glr_auto_examples_01_linear_ode.py and every other example: they all pick n_steps up front and never revisit it). Two pieces turn that into a &quot;run a chunk, look, continue if not converged&quot; workflow:">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_16_convergence_drift_thumb.png
    :alt:

  :doc:`/auto_examples/16_convergence_drift`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Convergence diagnostics: convergence_drift + resume</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="sphx_glr_auto_examples_16_convergence_drift.py showed the &quot;run a chunk, look, continue if not converged&quot; workflow for a plain ODE. The same workflow applies to a fixed-delay DDE via lyapunov_spectrum_dde(..., resume=...) -- with one difference under the hood: an ODE&#x27;s Markovian state is just state, so its checkpoint (``lyapax.core.LyapunovCheckpoint``) only needs to carry state and the tangent basis Y. A DDE&#x27;s Markovian state is (state, buf) together -- buf is the ring buffer holding the recent history the delayed term reads from -- so lyapax.dde.DDECheckpoint additionally carries buf, its own tangent basis Y_buf, and the ring-buffer step counter t. Dropping any of those would mean the next chunk starts from an incomplete history and silently gets a biased exponent, not a crash.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_17_dde_resume_thumb.png
    :alt:

  :doc:`/auto_examples/17_dde_resume`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Resuming a DDE run: convergence_drift + resume, ring-buffer edition</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="lyapunov_spectrum is built entirely from jax.lax.scan and jnp.linalg.qr, both of which support ordinary JAX autodiff -- so a Lyapunov exponent can be differentiated w.r.t. a system parameter with jax.grad or jax.jacfwd, no different from any other JAX function. The practical payoff: instead of sweeping a parameter and hunting for where an exponent crosses a target value, gradient descent finds it directly -- and scales to many parameters at once, where a sweep does not.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_18_differentiate_lyapunov_exponent_thumb.png
    :alt:

  :doc:`/auto_examples/18_differentiate_lyapunov_exponent`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Differentiating a Lyapunov exponent w.r.t. a system parameter</div>
    </div>


.. raw:: html

    <div class="sphx-glr-thumbcontainer" tooltip="The Kaplan-Yorke dimension (Kaplan &amp; Yorke, 1979) estimates a chaotic attractor&#x27;s fractal dimension directly from its Lyapunov spectrum, with no extra simulation: sort the exponents descending, walk the cumulative sum lambda_1, lambda_1 + lambda_2, ... until it would go negative, and interpolate the fractional part from where it crosses zero. It&#x27;s a standard companion metric to a Lyapunov spectrum -- pure post-processing of LyapunovResult.exponents, exposed here as lyapax.core.kaplan_yorke_dimension.">

.. only:: html

  .. image:: /auto_examples/images/thumb/sphx_glr_19_kaplan_yorke_dimension_thumb.png
    :alt:

  :doc:`/auto_examples/19_kaplan_yorke_dimension`

.. raw:: html

      <div class="sphx-glr-thumbnail-title">Kaplan-Yorke (Lyapunov) dimension</div>
    </div>


.. thumbnail-parent-div-close

.. raw:: html

    </div>


.. toctree::
   :hidden:

   /auto_examples/00_time_series_sanity_check
   /auto_examples/01_linear_ode
   /auto_examples/01b_custom_system
   /auto_examples/02_chaotic_maps
   /auto_examples/03_chaotic_flows
   /auto_examples/04_linear_network
   /auto_examples/05_kuramoto_sync
   /auto_examples/06_custom_coupling
   /auto_examples/07_speed_and_accuracy
   /auto_examples/08_delayed_coupling
   /auto_examples/09_kuramoto_delayed_network
   /auto_examples/10_matrix_free_scaling
   /auto_examples/11_vmap_parameter_sweep
   /auto_examples/12_public_api_overview
   /auto_examples/13_dde_history_interpolation
   /auto_examples/14_gpu_acceleration
   /auto_examples/15_adaptive_ode
   /auto_examples/16_convergence_drift
   /auto_examples/17_dde_resume
   /auto_examples/18_differentiate_lyapunov_exponent
   /auto_examples/19_kaplan_yorke_dimension



.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_