<commit_before>.. _nonlinear:

******************
Nonlinear dynamics
******************

Introduction
============

This module contains utilities for nonlinear dynamics. Currently, the only
implemented features enable you to compute estimates of the maximum
Lyapunov exponent for an orbit. In future releases, there will be features
for creating surface of sections and computing the full Lyapunov spectrum.

Some imports needed for the code below::

    import astropy.units as u
    import numpy as np
    import gary.potential as gp
    import gary.dynamics as gd
    from gary.units import galactic

Computing Lyapunov exponents
============================

There are two ways to compute Lyapunov exponents implemented in `gary.dynamics`.
In most cases, you'll want to use the `~gary.dynamics.fast_lyapunov_max` function
because the integration is implemented in C and is quite fast. This function only
works if the potential you are working with is implemented in C (e.g., it is a
`~gary.potential.CPotentialBase` subclass). With a potential object and
a set of initial conditions::

    >>> pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, r_h=0.1*u.kpc,
    ...                               q1=1., q2=0.8, q3=0.6, units=galactic)
    >>> w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,5.5]*u.kpc,
    ...                                     vel=[0.,100.,0]*u.km/u.s)
    >>> lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., nsteps=100000) # doctest: +SKIP

This returns two objects: an `~astropy.units.Quantity` object that
contains the maximum Lyapunov exponent estimate for each offset orbit,
(we can control the number of offset orbits with the ``noffset_orbits``
argument) and an `~gary.dynamics.CartesianOrbit` object that contains
the parent orbit and each offset orbit. Let's plot the parent orbit::

    >>> fig = orbit[:,0].plot(marker=',', alpha=0.1, linestyle='none') # doctest: +SKIP

.. plot::
    :align: center

    import astropy.units as u
    import matplotlib.pyplot as pl
    import gary.potential as gp
    import gary.dynamics as gd
    from gary.units import galactic

    pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, r_h=0.1*u.kpc,
                                  q1=1., q2=0.8, q3=0.6, units=galactic)
    w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,5.5]*u.kpc,
                                        vel=[0.,100.,0]*u.km/u.s)
    lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., nsteps=100000)
    fig = orbit[:,0].plot(marker=',', linestyle='none', alpha=0.1)

Visually, this looks like a chaotic orbit. This means the Lyapunov exponent
should saturate to some value. We'll now plot the estimate of the Lyapunov
exponent as a function of time -- because the algorithm re-normalizes every
several time-steps (controllable with the ``nsteps_per_pullback`` argument),
we have to down-sample the time array to align it with the Lyapunov exponent
array. This plots one line per offset orbit::

    >>> pl.figure() # doctest: +SKIP
    >>> pl.loglog(orbit.t[11::10], lyap, marker=None) # doctest: +SKIP
    >>> pl.xlabel("Time [{}]".format(orbit.t.unit)) # doctest: +SKIP
    >>> pl.ylabel(r"$\lambda_{{\rm max}}$ [{}]".format(lyap.unit)) # doctest: +SKIP
    >>> pl.tight_layout() # doctest: +SKIP

.. plot::
    :align: center

    import astropy.units as u
    import matplotlib.pyplot as pl
    import gary.potential as gp
    import gary.dynamics as gd
    from gary.units import galactic

    pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, r_h=0.1*u.kpc,
                                  q1=1., q2=0.8, q3=0.6, units=galactic)
    w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,5.5]*u.kpc,
                                        vel=[0.,100.,0]*u.km/u.s)
    lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., nsteps=100000)

    pl.figure()
    pl.loglog(orbit.t[11::10], lyap, marker=None)
    pl.xlabel("Time [{}]".format(orbit.t.unit))
    pl.ylabel(r"$\lambda_{{\rm max}}$ [{}]".format(lyap.unit))
    pl.tight_layout()

The estimate is clearly starting to diverge from a simple power law decay.
To compare, we will compute the estimate for a regular orbit as well::

    >>> w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,0.]*u.kpc,
    ...                                     vel=[0.,140.,25]*u.km/u.s)
    >>> lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., nsteps=100000) # doctest: +SKIP
    >>> fig = orbit[:,0].plot(marker=',', alpha=0.1, linestyle='none') # doctest: +SKIP

.. plot::
    :align: center

    import astropy.units as u
    import matplotlib.pyplot as pl
    import gary.potential as gp
    import gary.dynamics as gd
    from gary.units import galactic

    pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, r_h=0.1*u.kpc,
                                  q1=1., q2=0.8, q3=0.6, units=galactic)
    w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,0.]*u.kpc,
                                        vel=[0.,140.,25]*u.km/u.s)<commit_msg>subsections<commit_after>.. _nonlinear:

******************
Nonlinear dynamics
******************

Introduction
============

This module contains utilities for nonlinear dynamics. Currently, the only
implemented features enable you to compute estimates of the maximum
Lyapunov exponent for an orbit. In future releases, there will be features
for creating surface of sections and computing the full Lyapunov spectrum.

Some imports needed for the code below::

    import astropy.units as u
    import numpy as np
    import gary.potential as gp
    import gary.dynamics as gd
    from gary.units import galactic

Computing Lyapunov exponents
============================

Chaotic orbit
-------------

There are two ways to compute Lyapunov exponents implemented in `gary.dynamics`.
In most cases, you'll want to use the `~gary.dynamics.fast_lyapunov_max` function
because the integration is implemented in C and is quite fast. This function only
works if the potential you are working with is implemented in C (e.g., it is a
`~gary.potential.CPotentialBase` subclass). With a potential object and
a set of initial conditions::

    >>> pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, r_h=0.1*u.kpc,
    ...                               q1=1., q2=0.8, q3=0.6, units=galactic)
    >>> w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,5.5]*u.kpc,
    ...                                     vel=[0.,100.,0]*u.km/u.s)
    >>> lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., nsteps=100000) # doctest: +SKIP

This returns two objects: an `~astropy.units.Quantity` object that
contains the maximum Lyapunov exponent estimate for each offset orbit,
(we can control the number of offset orbits with the ``noffset_orbits``
argument) and an `~gary.dynamics.CartesianOrbit` object that contains
the parent orbit and each offset orbit. Let's plot the parent orbit::

    >>> fig = orbit[:,0].plot(marker=',', alpha=0.1, linestyle='none') # doctest: +SKIP

.. plot::
    :align: center

    import astropy.units as u
    import matplotlib.pyplot as pl
    import gary.potential as gp
    import gary.dynamics as gd
    from gary.units import galactic

    pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, r_h=0.1*u.kpc,
                                  q1=1., q2=0.8, q3=0.6, units=galactic)
    w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,5.5]*u.kpc,
                                        vel=[0.,100.,0]*u.km/u.s)
    lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., nsteps=100000)
    fig = orbit[:,0].plot(marker=',', linestyle='none', alpha=0.1)

Visually, this looks like a chaotic orbit. This means the Lyapunov exponent
should saturate to some value. We'll now plot the estimate of the Lyapunov
exponent as a function of time -- because the algorithm re-normalizes every
several time-steps (controllable with the ``nsteps_per_pullback`` argument),
we have to down-sample the time array to align it with the Lyapunov exponent
array. This plots one line per offset orbit::

    >>> pl.figure() # doctest: +SKIP
    >>> pl.loglog(orbit.t[11::10], lyap, marker=None) # doctest: +SKIP
    >>> pl.xlabel("Time [{}]".format(orbit.t.unit)) # doctest: +SKIP
    >>> pl.ylabel(r"$\lambda_{{\rm max}}$ [{}]".format(lyap.unit)) # doctest: +SKIP
    >>> pl.tight_layout() # doctest: +SKIP

.. plot::
    :align: center

    import astropy.units as u
    import matplotlib.pyplot as pl
    import gary.potential as gp
    import gary.dynamics as gd
    from gary.units import galactic

    pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, r_h=0.1*u.kpc,
                                  q1=1., q2=0.8, q3=0.6, units=galactic)
    w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,5.5]*u.kpc,
                                        vel=[0.,100.,0]*u.km/u.s)
    lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., nsteps=100000)

    pl.figure()
    pl.loglog(orbit.t[11::10], lyap, marker=None)
    pl.xlabel("Time [{}]".format(orbit.t.unit))
    pl.ylabel(r"$\lambda_{{\rm max}}$ [{}]".format(lyap.unit))
    pl.tight_layout()

The estimate is clearly starting to diverge from a simple power law decay.

Regular orbit
-------------

To compare, we will compute the estimate for a regular orbit as well::

    >>> w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,0.]*u.kpc,
    ...                                     vel=[0.,140.,25]*u.km/u.s)
    >>> lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=2., nsteps=100000) # doctest: +SKIP
    >>> fig = orbit[:,0].plot(marker=',', alpha=0.1, linestyle='none') # doctest: +SKIP

.. plot::
    :align: center

    import astropy.units as u
    import matplotlib.pyplot as pl
    import gary.potential as gp
    import gary.dynamics as gd
    from gary.units import galactic

    pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, r_h=0.1*u.kpc,
                                  q1=1., q2=0.8, q3=0.6, units=galactic)
    w0 = gd.CartesianPhaseSpacePosition(pos=[5.5,0.,0.]*u.kpc,
                                        vel=[0.,140.,25]*u.km/u.s)
{'fraction_non_alphanumeric': 0.14474881229151926, 'fraction_numerical': 0.03153745072273324, 'mean_word_length': 2.8891509433962264, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 75, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '745243', 'n_tokens_mistral': 3536, 'n_tokens_neox': 3246, 'n_words': 927}