# 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? – R.McGuigan Feb 26 '18 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). – stafusa Feb 26 '18 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! – R.McGuigan Feb 26 '18 at 18:09

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. – Wrzlprmft Feb 26 '18 at 15:22
• Also, this is for continuous-time systems, while the question is about maps. – Wrzlprmft Feb 26 '18 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? – George Datseris Feb 27 '18 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. – George Datseris Feb 27 '18 at 16:22
• I looked into the literature again and it is indeed a confusing mess: The method usually referred to as Benettin’s is from Benettin et al (1980) and works in tangent space, which has several advantages over orbit separation. That paper states that “the method for computing the maximal [Lyapunov exponent] was already given in [Benettin et al (1976)]” which is not true, as the method given there is based on orbit separation and thus different. (I suppose they confused definite and indefinite articles.) – Wrzlprmft Feb 27 '18 at 20:28