I have a numerical model that I would like to calculate the Lyapunov spectrum for. The number of dimensions is in the hundreds, and I don’t have an analytical expression for the Jacobian available. The system is chaotic (positive largest Lyapunov exponent) and the fractal dimension is not huge (3–5). The numerical integration is done with an explicit Euler method.
What would be the best/simplest way to calculate the full spectrum of Lyapunov exponents? Does anyone know of Matlab scripts I could use and/or adapt? My goal is to calculate the Kaplan–Yorke dimension, and determine if the system is hyperchaotic.
You almost certainly do not want the entire Lyapunov spectrum, but only the largest two, three, or maybe four Lyapunov exponents. Strictly speaking, two suffice to show hyperchaos, but you may want to go for the first zero one as a sanity check. Getting, say, the hundredth Lyapunov exponent requires a lot of computation time and your result will likely be very inaccurate. Most importantly, you do not gain anything from it. You do not need it for the Kaplan–Yorke dimension either.
That being said, if you do not want to work with the Jacobian, orbit separation is the way to go. More specifically, you integrate many slightly different versions of the system and observe how their separation (the Lyapunov vectors) evolves. For Lyapunov exponents beyond the first one, you have to take care to regularly remove any components in the direction of the previous Lyapunov vectors.
Technically, the procedure for the first two Lyapunov exponents is:
Select some parameters $ε$, representing the size of a small perturbation, and $τ$ representing the rescaling interval. For most systems, $ε=10^{-12}$ is a good choice. Use a handful of oscillations of your dynamics for $τ$, but more about this choice later.
Create three instances of your system, and call denote their states at time $t$ by $y_0(t)$, $y_1(t)$, and $y_2(t)$.
Set the initial condition $y_0(0)$ to something on the attractor. Set $y_1(0)=y_0(0) +ε·r_1$ and $y_2(0)=y_0(0) +ε·r_2$, where $r_1$ and $r_2$ are random vectors.
Integrate for $τ$ time units.
Compute $$ v_1(t) = \frac{y_1(t)-y_0(t)}{|y_1(t)-y_0(t)|}; \qquad λ_1(t) = \frac{1}{τ}\log \left( \frac{|y_1(t)-y_0(t)|}{ε} \right).$$
$v_1$ is your first Lyapunov vector. $λ_1$ is your first local Lyapunov exponent. If $|y_1(t)-y_0(t)|$ is not much smaller than your attractor’s diameter, your choice of $τ$ is too high.
Compute $$ v_2(t) = \frac{y_2(t)-y_0(t)-\langle y_2, v_1 \rangle v_1}{\left | y_2(t)-y_0(t)-\langle y_2, v_1 \rangle v_1 \right|}; \qquad λ_2(t) = \frac{1}{τ} \log \left(\frac{\left| y_2(t)-y_0(t)-\langle y_2, v_1 \rangle v_1 \right|}{ε} \right).$$
$v_2$ is your second Lyapunov vector. $λ_2$ is your second local Lyapunov exponent. Note how this is analogous to Step 5 except for removing the projection on $v_1$.
Set $y_1(t) = y_0(t) + εv_1(t)$ and $y_2(t) = y_0(t) + εv_2(t)$. This rescaling is to ensure that your separations stay small.
Go to step 4.
After this you can average the $λ_1$ and $λ_2$ (except the first few ones) to obtain the respective Lyapunov exponents.
You could also estimate the Jacobian from finite differences, but that’s just replacing one $ε$ with another.