galpy.orbit.Orbit.lyapunov

Orbit.lyapunov(ts, pot=None, method='dop853_c', renorm_every=10, dxdv0=None, spectrum=False, progressbar=False, dt=None, numcores=2, force_map=False, rtol=None, atol=None, **kwargs)[source]

Estimate the largest Lyapunov exponent (or the full Lyapunov spectrum) using the variational equations.

The deviation vector is propagated with the linearized (variational) equations of motion (see integrate_dxdv) and renormalized to its initial norm every renorm_every output times (Benettin et al. 1980 method). The running estimate

\[\lambda(t) = \frac{1}{t-t_0}\,\sum \ln \frac{|\delta x|}{|\delta x_0|}\]

is returned at all output times, such that its convergence can be assessed; the final value is the best estimate of the largest Lyapunov exponent. For regular orbits, \(\lambda(t)\) decays as \(\ln(t)/t\), while for chaotic orbits it converges to a positive value.

When spectrum=True, the full Lyapunov spectrum (all phasedim exponents) is computed instead, using the classic generalization of the renormalization method (Benettin et al. 1980, part 2; Shimada & Nagashima 1979): an orthonormal set of phasedim deviation vectors is propagated simultaneously and re-orthonormalized with a QR decomposition every renorm_every output times; the running estimate of the \(i\)-th exponent is the accumulated \(\sum \ln |R_{ii}|\) of the QR growth factors divided by the elapsed time.

Parameters:

  • ts (list, numpy.ndarray or Quantity) – List of equispaced times at which to compute the running Lyapunov estimate. The initial condition is ts[0]. (note that for method=’odeint’, method=’dop853’, and method=’dop853_c’, the time array can be non-equispaced).

  • pot (Potential, DissipativeForce, or a combined force/potential formed using addition (pot1+pot2+force3+…), optional) – Gravitational field to integrate the orbit in. Default is the gravitational field used in the last orbit integration.

  • method (str, optional) – Integration method; see integrate_dxdv for the possible (non-symplectic) methods. Default is ‘dop853_c’.

  • renorm_every (int, optional) – Renormalize the deviation vector to its initial norm every renorm_every output time intervals. Default is 10.

  • dxdv0 (numpy.ndarray, optional) – Initial phase-space deviation vector in rectangular coordinates [dx,dy,dvx,dvy] (planar orbits) or [dx,dy,dz,dvx,dvy,dvz] (3D orbits); shape (phasedim,) [same deviation for all orbits] or (*input_shape, phasedim). Because the variational equations are linear, the overall normalization is irrelevant. Default is the unit vector with equal components along all phase-space directions. When spectrum=True, the (normalized) dxdv0 is used as the first deviation vector and completed to an orthonormal basis; the default initial basis is the identity.

  • spectrum (bool, optional) – If True, compute the full Lyapunov spectrum (all phasedim exponents) instead of only the largest exponent; see the Returns section for the output format. Default is False.

  • progressbar (bool, optional) – If True, display a tqdm progress bar when integrating multiple orbits (requires tqdm to be installed!); shown for each renormalization segment. Default is False.

  • dt (float, optional) – If set, force the integrator to use this basic stepsize; must be an integer divisor of output stepsize (only works for the C integrators that use a fixed stepsize) (can be Quantity).

  • numcores (int, optional) – Number of cores to use for Python-based multiprocessing (pure Python or using force_map=True); default = OMP_NUM_THREADS.

  • force_map (bool, optional) – If True, force use of Python-based multiprocessing (not recommended). Default is False.

  • rtol (float, optional) – Relative tolerance. Default is None.

  • atol (float, optional) – Absolute tolerance. Default is None.

  • ro (float or Quantity, optional) – Physical scale in kpc for distances to use to convert. Default is object-wide default.

  • vo (float or Quantity, optional) – Physical scale for velocities in km/s to use to convert. Default is object-wide default.

  • use_physical (bool, optional) – Use to override object-wide default for using a physical scale for output.

  • quantity (bool, optional) – If True, return an Astropy Quantity object. Default from configuration file.

Returns:

For spectrum=False (default): running estimate \(\lambda(t)\) of the largest Lyapunov exponent at the output times. For spectrum=True: running estimates \(\lambda_i(t)\) of all phasedim Lyapunov exponents, such that out[...,i,:] is the running estimate of the \(i\)-th exponent; the exponents are ordered from largest to smallest by construction of the algorithm (asymptotically; finite-time estimates of nearly-degenerate exponents may transiently swap). In both cases, the entry for the initial time ts[0] is NaN (no elapsed time yet).

Return type:

numpy.ndarray or Quantity [*input_shape,nt] or [*input_shape,phasedim,nt]

Notes

  • Only implemented for 4D (planar) and 6D (3D) orbits.

  • The previously integrated orbit stored in this instance (if any) is not affected.

  • For spectrum=True, the absolute values \(|R_{ii}|\) of the diagonal of the QR decomposition are used as the growth factors, so the sign convention of the QR routine is irrelevant. For the Hamiltonian flows integrated here, the exponents obey \(\sum_i \lambda_i = 0\) (phase-space volume conservation) and come in \(\pm\) pairs (\(\lambda_i = -\lambda_{\mathrm{phasedim}+1-i}\)); the deviation of the computed spectrum from these relations is a useful convergence/accuracy check.

  • References: Benettin et al. (1980, Meccanica 15, 9 and Meccanica 15, 21); Shimada & Nagashima (1979, Prog. Theor. Phys. 61, 1605); see also Skokos (2010, Lect. Notes Phys. 790, 63).

  • 2026-06-09 - Written - Bovy (UofT)

  • 2026-06-10 - Added spectrum= option - Bovy (UofT)