¹

An Elementary Definition

Note this just means a sink is when we can shrink any open ball around it to an arbitrarily small open ball by taking sufficiently large steps. Similarily source is a point when we can expand arbitrary open ball so that everything except the point itself goes outside the ball by taking sufficiently large steps.

We can simulate the action of Hénon map by starting with a grid of points and apply them for many iterations. Writing a small julia function:

function compute_attractor(a, b; N=1000, T=1000, h=100, w=100) 
    states = zeros(N, N, 2)
    for i = 1:N
        for j = 1:N
            for t = 1:T
                if t == 1
                    states[i,j,1] = -h/2 + (i-1)*h/N
                    states[i,j,2] = -w/2 + (j-1)*w/N
                else 
                    x = states[i,j, 1]
                    y = states[i,j, 2]
                    if (abs(x) < 1e2) && (abs(y) < 1e2)
                        states[i,j, 1] = a - x^2 + b*y 
                        states[i,j, 2] = x
                    else 
                        states[i, j, 1] = Inf
                        states[i, j, 2] = Inf
                    end
                end
            end
        end
    end
    return states 
end

For \(a=1.4, b=0.3\), the image looks like:

This gets an approximation of the attractor, which is a generalization of fixed point.

Now let’s see some continous-time examples

Another example is the Lorenz system, which originates from climate sciences(metereology), which arise after applying spectral methods to solve Oberbeck-Boussinesq approximation of the Navier-Stokes equation.

import Pkg; Pkg.add("GLMakie")
using GLMakie

Base.@kwdef mutable struct Lorenz
    dt::Float64 = 0.01
    σ::Float64 = 10
    ρ::Float64 = 28
    β::Float64 = 8/3
    x::Float64 = 1
    y::Float64 = 1
    z::Float64 = 1
end

function step!(l::Lorenz)
    dx = l.σ * (l.y - l.x)
    dy = l.x * (l.ρ - l.z) - l.y
    dz = l.x * l.y - l.β * l.z
    l.x += l.dt * dx
    l.y += l.dt * dy
    l.z += l.dt * dz
    Point3f(l.x, l.y, l.z)
end

attractor = Lorenz()

points = Point3f[]
colors = Int[]


fig, ax, l = lines(points, color = colors,
    colormap = :inferno, transparency = true,
    axis = (; type = Axis3, protrusions = (0, 0, 0, 0),
              viewmode = :fit, limits = (-30, 30, -30, 30, 0, 50)))

record(fig, "assets/videos/lorenz.gif", 1:120) do frame
    for i in 1:50
        push!(points, step!(attractor))
        push!(colors, frame)
    end
    ax.azimuth[]=atan(frame/90)
    Makie.update!(l, arg1 = points, color = colors) 
    l.colorrange = (0, frame)
end

Lyapuov Exponent

Consider first a discrete dynamical system \((X,f)\). We are interested in how a perturbation \(\delta x\) in the initial state \(x_0\) would impact the trajectory. By linear approximation:

where \(x_n = f^n(x_0)\) and \(x_n'\) denote the unperturbed and perturbed states at time \(n\). Similarily:

Suppose we denote \(u_n=x_n'-x_n\), by induction, we have that:

Or if \(f\) is time-dependent, \(x_{k+1} = f(x_k, k)\), we have

We’ll denote the product of matrices by ²

so that \(u_n \approx H(x_0, n)\, \delta x\); for autonomous \(f\), \(H(x_0, n) = Df(x_{n-1}) \cdots Df(x_0) = Df^n(x_0)\) by the chain rule.

Continuous Time

For a continuous system \(\frac{dx}{dt} = F(x, t)\) with trajectory \(x(t)\) starting at \(x_0\), a perturbation \(u(t) = x'(t) - x(t)\) with \(u(0) = \delta x\) obeys, to first order, the linear equation

where \(J(t)\) is the Jacobian of \(F\) along the trajectory. Its solution defines the continuous analog of \(H\):

Time-ordered Exponential. \(\mathcal{T}\exp\) is not the matrix exponential of \(\int_0^t J(s)\, ds\). Solving the equation above by repeated substitution gives the series

Now consider the time-ordering \(\mathcal{T}\) permutes each product so that later times stand to the left:

Note for each \((s_1,s_2,...)\) the \(\mathcal{T}\) enforces a possibly different order, i.e., \(\mathcal T\) is point-wise.

Now we claim that:

Editing notes: Need a proof here

However, if the jacobian \(J_j(s_j)\) commute with each other, \(\mathcal T\) would be redundant.

Some examples:

• Linear systems \(F(x) = Ax\) have \(J(s) = A\) for all \(s\), so \(H(x_0, t) = e^{tA}\).
• One-dimensional systems have scalar Jacobians, which always commute, so \(H(x_0, t) = \exp\left(\int_0^t F'(x(s))\, ds\right)\).
• Determinants turn matrices into numbers, so \(\det H\) always behaves like the scalar case. This is Liouville’s formula: \(\det H(x_0, t) = \exp\left(\int_0^t \operatorname{tr} J(s)\, ds\right)\).
• The Lorenz system has Jacobians that do not commute at different points, so its \(H\) is a genuine time-ordered exponential.

The Oseledets Theorem

In words: fix a typical starting point \(x\) (typical for \(\mu\)) and perturb it in a direction \(v\). The perturbation grows or shrinks like \(e^{n\lambda}\), where the rate \(\lambda\) is one of finitely many values \(\lambda_1 > \dots > \lambda_s\) and depends only on the direction: it is \(\lambda_i\) exactly when \(v \in V_i(x) \setminus V_{i+1}(x)\) (item 1). A randomly chosen direction gets the largest rate \(\lambda_1\), because \(V_2(x)\) is a lower-dimensional subspace. The rates belong to the whole orbit, not just to the point \(x\) (item 2), and they are also the growth rates of the singular values of \(T^n(x)\) (item 3).

Computing the Spectrum Numerically

The theorem suggests an algorithm (Benettin et al., 1980): follow a frame of \(d\) perturbation vectors along the trajectory and measure how fast it stretches. Two problems arise if the frame is evolved naively: every column turns toward the fastest direction \(V_1 \setminus V_2\), and the entries overflow. The fix is to re-orthonormalize with a QR decomposition every few steps. The \(k\)-th diagonal entry of \(R\) is the growth of the \(k\)-th column after the first \(k-1\) directions have been removed, so its logarithm, averaged over time, converges to \(\lambda_k\).

input: vector field F, Jacobian DF, initial state x, step Δt, block length m, number of blocks N
U ← I                                  # orthonormal frame of d perturbations
S ← 0                                  # accumulated log-stretching, one entry per exponent
run x forward for a while so it lands on the attractor
repeat N times:
    repeat m times:                    # one Euler step of x and of the frame
        x, U ← x + Δt F(x),  U + Δt DF(x) U
    Q, R ← qr(U)                       # U = QR, Q orthonormal, R upper triangular
    S ← S + log |diag(R)|              # how much each direction stretched
    U ← Q                              # restart from an orthonormal frame
return S / (N m Δt)                    # λ_1, ..., λ_d

The frame update \(U \leftarrow U + \Delta t\, DF(x)\, U = (I + \Delta t\, J)\, U\) is one Euler step of \(\frac{d}{dt}H = JH\), so between re-orthonormalizations \(U\) tracks \(H(x_0, t)\) applied to the frame.

# A script to compute the lyapunov spectrum for the lorenz system 
# @ author : Haoyu Tang 

using LinearAlgebra

Base.@kwdef mutable struct State
    x::Float64 = 1
    y::Float64 = 1
    z::Float64 = 1
    U::Matrix{Float64} = Matrix{Float64}(I,3,3)
end

Base.@kwdef mutable struct Lorenz
    dt::Float64 = 0.001
    σ::Float64 = 10
    R::Float64 = 28
    b::Float64 = 8/3
    v::State = State()
end

function compute_attractor!(L::Lorenz; N=100)
    v, σ, R, b, dt = L.v, L.σ, L.R, L.b, L.dt
    for T = 1:N
        v.x += σ * (v.y - v.x) * dt
        v.y += (R * v.x - v.y - v.x * v.z) * dt
        v.z += (- b * v.z + v.x * v.y) * dt
    end
end

function rhs(L::Lorenz)
    v, σ, R, b = L.v, L.σ, L.R, L.b
    x, y, z, U = v.x, v.y, v.z, v.U
    J = [-σ  σ  0 ; (R - z)  -1 -x ; y x -b]
    dX = [σ * (y - x), R * x - y - x * z, - b * z + x * y]
    dU = J * U 
    return hcat(dX, dU)
end

function benettin!(L::Lorenz; block_size=1, total_blocks=1_000_000)
    # set up accumulator for lyapunov exponent 
    S = [0 0 0]
    
    # integrate to get an approximation of the trajectory
    compute_attractor!(L)
    
    dt = L.dt
    v = L.v

    for n = 1:total_blocks
        # euler step for block_size steps
        for t = 1:block_size
            dv = rhs(L)
            v.x += dv[1, 1] * dt
            v.y += dv[2, 1] * dt
            v.z += dv[3, 1] * dt
            v.U += dv[:, 2:4] * dt
        end

        # reorthonormalization 
        q, r = qr(v.U)
        
        # adjust sign
        sign_vector = Diagonal(sign.(diag(r)))
        q = q * sign_vector
        r = r * sign_vector
        
        v.U = q 
        S += log.(diag(r)') 
    end
    return S ./ (total_blocks * block_size * dt)
end 

L = Lorenz()
λ = benettin!(L)

expected_sum = -(L.σ + L.b + 1)
computed_sum = sum(λ)

diff = (abs(expected_sum - computed_sum))/abs(computed_sum) * 100
println()
print("Expecting $expected_sum, got $computed_sum, error = $diff %")
println()
show(stdout, "text/plain", λ)

With \(\Delta t = 0.001\) over \(1000\) time units this prints \((\lambda_1, \lambda_2, \lambda_3) \approx (0.906, 0.001, -14.608)\), close to the reference values \((0.906, 0, -14.572)\). Three checks tell us whether to trust the result:

• The sum is \(-(\sigma + 1 + \beta)\). The frame’s volume is \(\lvert\det U\rvert\), which grows like \(e^{(\lambda_1 + \lambda_2 + \lambda_3)t}\). By Liouville’s formula it also equals \(\exp\left(\int_0^t \operatorname{tr} J\, ds\right)\), and \(\operatorname{tr} J = -\sigma - 1 - \beta\) at every point. So the sum must be \(-(\sigma + 1 + \beta) \approx -13.667\). The code gives \(-13.70\): the gap comes from Euler’s method and shrinks with \(\Delta t\). This check only sees the diagonal of \(J\), so passing it does not mean the individual exponents are right.
• One exponent is \(0\). A perturbation along the flow direction \(F(x)\) is just a shift in time along the same trajectory (\(H(x, t)\, F(x) = F(f^t(x))\)), so it neither grows nor shrinks exponentially. Every flow whose attractor is not a fixed point has this zero exponent.
• A positive exponent means chaos. If \(\lambda_1 > 0\), two nearby initial states separate like \(\lvert\delta x\rvert\, e^{\lambda_1 t}\), so an error of size \(\delta\) grows to the size \(L\) of the attractor after about \(\frac{1}{\lambda_1}\log\frac{L}{\delta}\) time units. After that, prediction is impossible. The attractor is bounded, so this stretching must be accompanied by folding, which is what produces its fractal structure. This sensitive dependence on initial conditions on a bounded attractor is what we call chaos. For Lorenz, \(1/\lambda_1 \approx 1.1\) time units: each extra digit of precision in the initial state buys only about \(2.5\) more time units of prediction.