Finite Time Lyapunov Exponents
Mathematical Background
We change focus now to a continuous dynamical system, e.g. an ODE $\dot{u} = g(t, u)$ with solution $\Phi^{t,t_0} (x)$. Since $\Phi^{t,t_0} (x)$ is continuously dependent on the initial condition $x$, there exists an $\tilde{x}$ near $x$ with $sup_{t \in [t_0 , t_0 + T]} \| \Phi^{t,t_0} (\tilde{x}) - \Phi^{t,t_0} (x) \| < \epsilon$ for any fixed $\epsilon > 0$ and $T$ small enough. We wish to characterize this expansion term. We write $\Phi = \Phi^{t_0+T,t_0}$ and $y = x + \delta x_0$ where $\delta x_0 \in \mathbb{R}^d$ is infinitesimal. Then if $g$ is $\mathcal{C}^1$ w.r.t. $x$,
\[\delta x (t_0 + T) := \Phi (y) - \Phi (x) = D_x \Phi (x) \cdot \delta x_0 + \mathcal{O}(\| \delta x_0 \|^2)\]
Hence we can write
\[\| \delta x (t_0 + T) \|_2 = \| D_x \Phi (x) \cdot \delta x_0 \|_2 \leq \| D_x \Phi (x) \|_2 \cdot \| \delta x_0 \|_2\]
or equivalently
\[\frac{ \| \delta x (t_0 + T) \|_2 }{ \| \delta x_0 \|_2 } \leq \| D_x \Phi (x) \|_2\]
where equality holds if $\delta x_0$ is the eigenvector corresponding to the largest eigenvalue of
\[\Delta = \left( D_x \Phi (x) \right)^T \left( D_x \Phi (x) \right) . \]
Hence if we define
\[\sigma (x) = \sigma^{t_0 + T, t_0} (x) := \frac{1}{T} \ln \left( \sqrt{\lambda_{\text{max}}} (\Delta) \right) = \frac{1}{T} \ln \left( \sup_{\delta x_0} \frac{ \| \delta x (t_0 + T) \|_2 }{ \| \delta x_0 \|_2 } \right)\]
then
\[\| \delta x (t_0 + T) \|_2 \leq e^{T \cdot \sigma (x)} \cdot \| \delta x_0 \|_2 . \]
From this we see why $\sigma (x)$ is called the maximal finite-time lyapunov exponent (FTLE).
The definition of $\sigma (x)$ leads to a natural ansatz for approximating the FTLE: compute $\frac{1}{T} \ln \left( \sup_{\delta x_0} \frac{ \| \delta x (t_0 + T) \|_2 }{ \| \delta x_0 \|_2 } \right)$ for each of a set of test points $\| \delta x_0 \|$ of fixed order $\epsilon > 0$ and set $\sigma (x)$ to be the maximum over this set of test points.
An extension of this technique can be made for ergodic systems, as shown in [15]:
when calculating the maximal Lyapunov exponent for a discrete dynamical system $x_{n+1} = f(x_k)$ defined as
\[\lambda (x, v) = \lim_{n \to \infty} \frac{1}{n} \log \| Df^n (x) \cdot v \|\]
a known technique is to use a QR iteration. Let $A = Q(A) R(A)$ be the unique QR-decomposition of a nonsingular matrix $A$ into an orthogonal matrix $Q(A)$ and an upper-triangular matrix $R(A)$. Then from $\| Av \| = \| R(A) v \|$ we have
\[\lim_{n \to \infty} \frac{1}{n} \log \| Df^n(x) v \| = \lim_{n \to \infty} \frac{1}{n} \log \| R(Df^n(x))v \|. \]
Further, the Lyapunov exponents of the system $\lambda_1, \ldots, \lambda_d$ (which are costant over the phase space for an ergodic system) can be found via
\[\lambda_j = \lim_{n \to \infty} \frac{1}{n} \log R_{jj}(Df^n(x))\]
the $j$-th diagonal element of $R$, for $j = 1, \ldots, d$.
Using an extension of the Birkhoff ergodic theorem it can be proven that this method is equivalent to computing
\[\lambda_j = \lim_{n \to \infty} \int \log R_{jj}(Df^n(x)) \, d\mu\]
where $\mu$ is a measure which is ergodic and invariant under $f$.
GAIO.finite_time_lyapunov_exponents — Functionfinite_time_lyapunov_exponents(F::SampledBoxMap, boxset::BoxSet) -> BoxMeasureCompute the Finite Time Lyapunov Exponent for every box in boxset, where F represents a time-T integration of some continuous dynamical system. It is assumed that all boxes in boxset have radii of some fixed order ϵ.
finite_time_lyapunov_exponents(f, Df, μ::BoxMeasure; n=8) -> σCompute the Lyapunov exponents using a spatial integration method [1] based on Birkhoff's ergodic theorem. Computes
\[\sigma_j = \frac{1}{n} \int \log R_{jj}( Df^n (x) ) \, dμ (x), \quad j = 1, \ldots, d\]
with respect to an ergodic invariant measure $\mu$.
[1] Beyn, WJ., Lust, A. A hybrid method for computing Lyapunov exponents. Numer. Math. 113, 357–375 (2009). https://doi.org/10.1007/s00211-009-0236-4
Example
We will continue using the periodically driven double-gyre introduced in the section on Almost Invariant (metastable) Sets. See that code block for the definition of the map.
t₀, τ, steps = 0, 0.1, 20
t₁ = t₀ + τ * steps
Tspan = t₁ - t₀
Φₜ₀ᵗ¹(z) = Φ(z, t₀, τ, steps)
domain = Box((1.0, 0.5), (1.0, 0.5))
P = BoxPartition(domain, (256, 128))
S = cover(P, :)
F = BoxMap(:grid, Φₜ₀ᵗ¹, domain, n_points=(6,6))
γ = finite_time_lyapunov_exponents(F, S; T=Tspan)BoxMeasure in 256 x 128 - element BoxPartition with 32768 boxes in its suportusing Plots
p = plot(γ, clims=(0,2), colormap=:jet);
Since this map is time-dependent, the FTLE field will change over time as well.
n_frames = 120
times = range(t₀, t₁, length=n_frames)
anim = @animate for t in times
Φₜ(z) = Φ(z, t, τ, steps)
F = BoxMap(:grid, Φₜ, domain, n_points=(6,6))
γ = finite_time_lyapunov_exponents(F, S; T=Tspan)
plot(γ, clims=(0,2), colormap=:jet)
end;
gif(anim, "ftle1.gif", fps=20)
Example 2: An Ergodic System
using GAIO
using Plots
using StaticArrays
# the Henon map
a, b = 1.4, 0.3
f((x,y)) = SA{Float64}[ 1 - a*x^2 + y, b*x ]
Df((x,y)) = SA{Float64}[-2*a*x 1.;
b 0.]
center, radius = (0, 0), (3, 3)
P = BoxPartition(Box(center, radius))
F = BoxMap(f, P)
S = cover(P, :)
A = relative_attractor(F, S, steps = 20)
T = TransferOperator(F, A, A)
(λ, ev) = eigs(T)
μ = real ∘ ev[1]
σ16 = finite_time_lyapunov_exponents(f, Df, μ, n=16)
σ8 = finite_time_lyapunov_exponents(f, Df, μ, n=8)
σ = 2*σ16 - σ8