.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/01_linear_ode.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_01_linear_ode.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_01_linear_ode.py:


Linear ODE: exact Lyapunov spectrum
====================================

The simplest possible correctness check for the Benettin/QR engine: for a
linear system ``x' = 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 :ref:`Tier 0.1 <validation-tier-0-1>` of
the validation guide.

**The system.** ``A`` here is diagonal with entries -1, -2, -5, so the
trajectory decays to the origin along the coordinate axes, and the
eigenvalues (already real, since ``A`` is diagonal) are exactly -1, -2, -5.
Every solution contracts, so all three exponents are negative -- there is
no chaos to detect, which is the point: this case isolates the numerics of
the Lyapunov engine from any question about the dynamics themselves.

**The method (Benettin/QR).** Alongside the state ``x(t)``, lyapax evolves
a set of ``d`` tangent vectors under the linearized dynamics ``dY/dt = A
Y`` (for this linear system the "linearization" is just ``A`` itself).
Left alone, all tangent vectors collapse onto the single fastest-growing
(here, slowest-decaying) direction and numerically underflow. To prevent
that, every ``renorm_every`` steps the tangent matrix ``Y`` is
QR-decomposed, ``Y = Q R``; ``Q`` (an orthonormal basis) replaces ``Y`` for
the next stretch, and the log of the diagonal of ``R`` gives that
stretch's contribution to each exponent. Averaging those contributions
over time and dividing by elapsed time gives the running estimate plotted
below, which converges to ``Re(eigvals(A))`` as the average washes out
the initial transient.

.. GENERATED FROM PYTHON SOURCE LINES 32-47

.. code-block:: Python

    import os
    import time

    os.environ["JAX_PLATFORMS"] = "cpu"

    import jax
    import jax.numpy as jnp
    import matplotlib.pyplot as plt
    import numpy as np

    jax.config.update("jax_enable_x64", True)

    from lyapax import systems
    from lyapax.core import lyapunov_spectrum, ode_problem








.. GENERATED FROM PYTHON SOURCE LINES 48-53

Distinct real eigenvalues, no coupling, no chaos. ``ode_problem`` turns
the continuous-time ``rhs`` into a fixed-``dt`` update (``integrator=
"rk4"`` by default) and bundles ``state0``/``dt`` alongside it, so
``lyapunov_spectrum`` -- the Benettin/QR engine described above -- only
needs ``n_steps`` and the run controls, not a second copy of ``dt``.

.. GENERATED FROM PYTHON SOURCE LINES 53-71

.. code-block:: Python

    A = jnp.diag(jnp.array([-1.0, -2.0, -5.0]))
    rhs = systems.linear_system(A)
    dt = 1e-3
    problem = ode_problem(rhs, state0=jnp.array([0.3, -0.2, 0.5]), dt=dt)

    t0 = time.perf_counter()
    result = lyapunov_spectrum(
        problem, n_steps=20_000, renorm_every=10, t_transient=5.0,
    )
    elapsed = time.perf_counter() - t0

    expected = np.array([-1.0, -2.0, -5.0])
    estimate = np.array(result.exponents)
    print(f"exact eigenvalues:      {expected}")
    print(f"lyapax estimate:        {estimate}")
    print(f"max abs error:          {np.max(np.abs(estimate - expected)):.2e}")
    print(f"wall time (incl. JIT):  {elapsed:.2f}s")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    exact eigenvalues:      [-1. -2. -5.]
    lyapax estimate:        [-1.00000095 -1.99999905 -5.        ]
    max abs error:          9.45e-07
    wall time (incl. JIT):  0.67s




.. GENERATED FROM PYTHON SOURCE LINES 72-76

Convergence of the running estimate toward the exact eigenvalues. The
transient (first 5 time units, not shown) is where the randomly
initialized tangent vectors align to the eigendirections -- see the note
on this in src/lyapax/core.py.

.. GENERATED FROM PYTHON SOURCE LINES 76-88

.. code-block:: Python

    fig, ax = plt.subplots(figsize=(6, 4))
    history = np.array(result.history)
    times = np.array(result.times)
    for i, lam in enumerate(expected):
        ax.plot(times, history[:, i], label=rf"$\lambda_{i + 1}$ estimate", color=f"C{i}")
        ax.axhline(lam, color=f"C{i}", linestyle="--", alpha=0.5)
    ax.set_xlabel("time (post-transient)")
    ax.set_ylabel("running Lyapunov exponent estimate")
    ax.set_title("Linear ODE: convergence to exact eigenvalues")
    ax.legend()
    fig.tight_layout()
    plt.show()



.. image-sg:: /auto_examples/images/sphx_glr_01_linear_ode_001.png
   :alt: Linear ODE: convergence to exact eigenvalues
   :srcset: /auto_examples/images/sphx_glr_01_linear_ode_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 1.137 seconds)


.. _sphx_glr_download_auto_examples_01_linear_ode.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: 01_linear_ode.ipynb <01_linear_ode.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: 01_linear_ode.py <01_linear_ode.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: 01_linear_ode.zip <01_linear_ode.zip>`


.. only:: html

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

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