I am trying to determine the evolution of a non-linear numerical model through time given certain perturbations in the initial conditions (I have the strong suspicion the system is dissipative, but would like to quantify this). As I understand, this can be done via Lyapunov exponents, approximated through he following formula for finite time series (the actual definition seems to be in the limit of infinity):
$$ \lambda^{(n)} = \left(\frac{1}{n}\right)ln\left[\frac{d(n)}{d(0)}\right]$$
Where $\lambda^{(n)}$ is the lyapunov exponent after the $n$th iteration and $d$ is the distance between two points in phase space. This definition seems intuitive to me: if the distance between the points decreases with time (i.e., on average: $\frac{d(n)}{d(0)}$ < 1), then the system is likely dissipative ('forgetting' its initial conditions), and the lyapunov exponent will be correctly negative. If the ratio is larger than one, the exponent will be positive, the points 'drift apart', and the system behaves chaotically.
So far, this makes sense to me. Now, most publications I have found so far suggest the following procedure:
- Define an initial condition
- Perturb it slightly in all directions of state space
- Propagate the points through time, receive new position of the 'main' point and its equivalent perturbed positions
- Use these points to determine the Jacobian matrix
- Send different perturbed points on the orthonormal axes of a unit hypersphere through the Jacobian, receive points on a hyper-ellipsis
- Determine equivalent orthonormal axes, re-scale them to get another hypersphere, save the scaling factors for each axis as Laypunov factors
- Rinse and repeat
Now my question: Are these steps and the linearization implied through use of the Jacobian really necessary? Numerically calculating the Jacobian is rather expensive, as interrupting the numerical model after each time step inflates the computing time a lot.
Couldn't I just take the initial perturbed points, send them all $t$ time steps forward, and then investigate their Euclidian distances in phase space to the central point at each of these time points according to the formula above? Wouldn't that be much simpler? Or am I missing (or misunderstanding) something important?
In general, you can.
Your simple method is what is often named orbit separation. The method involving the Jacobian is (if I understand it correctly) what is sometimes referred to as a tangent-space method. The latter have some advantages:
While orbit separation deals with tiny distances, tangent-space methods deal with truly infinitesimally small perturbations and are thus immune against some confounding effects. This becomes particularly problematic if you want more than one Lyapunov exponent.
With orbit separation you have to rescale regularly to keep your distances small. This introduces two parameters: How often you rescale and how small you rescale. With tangent-space methods you only need to do this to avoid numerical over- or underflow. As long as you rescale often enough, it doesn’t affect the result.
Now, if you avoid all the pitfalls and your system is not pathological, you can get proper estimates using orbit separation. (In particular if you already implemented a tangent-space method, you can compare the two.)