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?
 A: 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.
A: 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:


*

*maximum lyapunov exponent.

*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.

*MLE from a numerical dataset.



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