# How to calculate the maximal Lyapunov exponent(s) of a multidimensional system?

I've been reading up about Lyapunov exponents for a university group project on chaos theory and I'm a little confused as to how they are calculated for a system of multiple dimensions.

From what I can tell, the maximal Lyapunov exponent $\lambda$ for some 1-d map $f(x_{n})=x_{n+1}$ is:

$\lambda \approx \frac{1}{n}\sum\limits_{i=0}^{n-1}ln|f'(x_{i})|$

Where, if I understand things correctly, $f'(x_{i})$ is the derivative of $f$ at the ith value of $x$.

What my question is, is how does this extend to the multi-dimensional case? Say I had a pair of coupled maps representing a 2d system, $f(x_{n}, y_{n})=x_{n+1}$ and $g(x_{n}, y_{n})=y_{n+1}$, how is the Lyapunov exponent found?

The extension to several dimensions is natural:

• instead of a scalar, the state variable is a vector, and
• the derivative is substituted by the system's Jacobian.

Some earlier papers for numerical methods for calculating the Lyapunov exponents are mentioned in the Wikipedia entry and a more complete entry is found in the Scholarpedia, but textbooks on nonlinear dynamics (e.g., Ott, 1st ed., Eq. 4.40) offer more easier-to-follow how-to's.

• Thanks, for the links, I've looked through the pages in Ott on the Lyapunov exponents (around page 130) and I'd like to verify a few things. For the aforementioned project we want to find the maximum Lyapunov exponent for different Algorithms/maps applied to the same chaotic differential equations and look at the difference in the exponents. Later in the section it says that the maximum exponent dominates for large n, so to calculate the maximum exponent, would you simply use an approximation of (4.31) for sufficiently high n? Commented Feb 26, 2018 at 17:25
• @R.McGuigan, yes, it's simpler if you only need the largest exponent, but I'd use the 4.40 (or rather the approximation that follows it). Commented Feb 26, 2018 at 18:06
• Yes, eq (4.40) does seem more suitable, it's exactly what we're looking for. Thanks a bunch for your answer! Commented Feb 26, 2018 at 18:09
• If one is familiar with german, the german wikipedia article about the Lyapunov exponent offers a lot more details about the spectra computation for multidimensional systems (though I think, there's a minor mistake for the finite time LE computation, the logarithm is missing there): de.wikipedia.org/wiki/Ljapunow-Exponent Commented Jan 6, 2021 at 21:46

There is only one maximal Lyapunov exponent (MLE). In my answer I will tell you how to calculate this, and not many Lyapunov exponents (because it is a different method).

The method I am presenting is due to Benettin [1], and simply evolves in parallel 2 trajectories while constantly rescaling one of the two.

## MLE calculation

Two neighboring trajectories with initial distance d0 are evolved in time. At time $t_i$ their distance $d(t_i)$ either exceeds either some upper_threshold, or is lower than some lower_threshold, which initializes a rescaling of the test trajectory back to having distance d0 from the given one, while the rescaling keeps the difference vector along the maximal expansion/contraction direction: $$u_2 = u_1 + (u_2 - u_1)/(d(t_i)/d_0)$$ where $u_2, u_1$ are the test and reference states respectively (both vectors).

The maximum Lyapunov exponent is the average of the time-local Lyapunov exponents $$\lambda = \frac{1}{t_{n}}\sum_{i=1}^{n} \ln\left( a_i \right),\quad a_i = \frac{d(t_{i})}{d_0}.$$

## Software that does it

DynamicalSystems.jl is a open-sourced software library that offers functions that can calculate:

1. maximum lyapunov exponent.
2. lyapunov spectrum (all Lyapunov exponents). The library also describes how this method works in detail, in case your question was about computing many Lyapunov exponents instead of only the maximum one. The method you describe about how to find the MLE of a 1D map can be expanded into the method described in the link.
3. MLE from a numerical dataset.

[1] : G. Benettin et al., Phys. Rev. A 14, pp 2338 (1976)

• The method you are describing is not Benettin’s – it’s the classical orbit-separation method. Benettin’s method works in tangent space with infinitesimal distances; the one you are describing uses finite distances. Commented Feb 26, 2018 at 15:22
• Also, this is for continuous-time systems, while the question is about maps. Commented Feb 26, 2018 at 15:51
• To the best of my knowledge Bennetin was the first one to introduce the method in a paper, which also happens to be the one I am citing here. Maybe take a look at a the paper before making a hasty comment? Commented Feb 27, 2018 at 16:21
• Your comment "this is for continuous-time systems, while the question is about maps" is just wrong. The exact same method works for both maps and flows. In fact the numerical implementation I cite has identical source code for maps and flows. Commented Feb 27, 2018 at 16:22
• Your comment "this is for continuous-time systems, while the question is about maps" is just wrong. – Strictly speaking, you are correct and you can apply this method to maps. However, for maps, there is the method described in Stafusa’s answer which is far superior since it is simpler, you do have to bother with rescaling, and it is more accurate since you do not suffer from finite-size effects. Only when you get to continuous-time systems, it makes sense to think about orbit-separation methods in the first place. Commented Feb 27, 2018 at 20:37