📂Dynamics

Systems Represented by Differential Equations and Their Lyapunov Spectrum with Numerical Calculation Methods

Algorithm

Assume that a vector field $X = \mathbb{R}^{n}$ for space and function $f : X \to X$ is given by a differential equation $$ \dot{x} = f(x) $$. Define a variational equation for the Jacobian matrix $J$ of $f$ as follows $$ \dot{Y} = J Y $$. Here, the initial condition for the matrix function $Y = Y(t) \in \mathbb{R}^{n \times n}$ is set as the identity matrix $Y(0) = I_{n}$. The Lyapunov spectrum $\lambda_{1} , \cdots , \lambda_{n}$ of the system $\dot{x} = f(x)$ is approximated by repeatedly solving the variational equation for the numerical solution (e.g., calculated by methods like RK4) $\left\{ x_{k} \right\}_{k=0}^{K}$ of $\dot{x} = f(x)$. $J(x)$ is dependent on the trajectory $x(t)$ of the original system.

• Solving the variational equation $\dot{Y} = J Y$ implies using a method like RK4, which was used to obtain the numerical solution of matrices or vectors, known as perturbations $Y$. If the sufficiently small time interval used to obtain $x_{k}$ is $h \approx 0$, the same time interval $h$ should be used for the variational equation.
• The Runge-Kutta method $RK4(J, U, h)$ is defined for matrices $J$ and $U$ as follows $$ \begin{align*} U(t + h) \approx RK4(J, U, h) =& {\frac{ h }{ 6 }} \left( V_{1} + 2 V_{2} + 2 V_{3} + V_{4} \right)\\ V_{1} =& J U \\ V_{2} =& J \left( U + {\frac{ h }{ 2 }} V_{1} \right) \\ V_{3} =& J \left( U + {\frac{ h }{ 2 }} V_{2} \right) \\ V_{4} =& J \left( U + h V_{3} \right) \end{align*} $$.

Now, introduce the traditional method using Gram-Schmidt orthogonalization and a method using QR decomposition based on geometric definitions.

Based on Gram-Schmidt Orthogonalization ¹

• Represent the $j$-th column vector of the matrix $U$ as $U_{:j}$.
• Represent the transpose matrix of the matrix $U$ as $U^{T}$.
• Represent the norm of the matrix $U$ as $\left\| U \right\|$.

Define a function $GS$ to obtain $U$ with normalized axis length and a matrix $V$ that is merely orthogonal without normalization, unlike the widely known Gram-Schmidt orthogonalization in linear algebra. Here, each column of $V$ in the variational equation $\dot{Y} = J Y$ represents how much the perturbation $Y$ stretches or compresses the axis length at specific points $x(t)$ of the original equation $\dot{x} = f(x)$ via the linear transformation $J \left( x(t) \right)$.

Based on QR Decomposition ²

Assume that a function $QR(A)$ for the matrix $A$ returns the results of the QR decomposition, $Q$ and $R$, of $A$.

Explanation

Calculating the Lyapunov spectrum for a dynamical system expressed by ODEs is much more challenging than calculating a bifurcation diagram. Since the concept of chaos itself is defined by the Lyapunov spectrum, it’s an extremely important skill for studying dynamics, and it’s also a part where one can feel a significant computational barrier. Naturally, there is a big difference between having a vague notion of Lyapunov exponents and actually being able to calculate them.

Although the two introduced methods may seem unrelated at first glance, a detailed look at the proof of the existence of QR decomposition reveals that they are fundamentally very similar. The differences can be summarized as follows:

• GS-based: Geometrically intuitive and easy to understand. Reading through the algorithm reveals that it’s essentially measuring the length of axes of an ellipsoid, transformed by $J$ from a matrix $Y$ having $n$ axes of a unit ball as column vectors.
• QR-based: Relatively easier to implement. After QR decomposition, just gather the diagonal elements of $R$, and since you’re likely using a package for QR decomposition from the start, the reliability of the program is high.

Comparison of Results

The results of calculating the Lyapunov spectrum of the Lorenz attractor at $\rho \in [0, 100]$ using two methods are as follows:

• Based on Gram-Schmidt:

• Based on QR decomposition:

• Reference ³

As seen, it’s hard to distinguish visually, obtaining similar results, and it’s also confirmed to be reliable when compared to the reference. Now, let’s compare as concrete data. Blue (left) is calculated based on Gram-Schmidt, while orange (right) is based on QR decomposition.

Performance Comparison

• Based on Gram-Schmidt (top): 2:39:31
• Based on QR decomposition (bottom): 9:50:20

However, the two techniques show considerable differences in speed. The result itself can be achieved by either method, but benchmarking to calculate these numerical values reveals a speed difference of about 3 to 4 times.

The following is a measure of the time taken to obtain the Lorenz attractor for a more accurate comparison. It took about 9 minutes, confirming that calculating the Lyapunov spectrum itself did not significantly affect it.

Code ^{4 5}

Full Code

Julia code calculating the Lyapunov spectrum from $\rho \in [0, 100]$ of the Lorenz attractor.

using CSV, DataFrames, LinearAlgebra, ProgressMeter, Plots, Dates; @info now()
using Base.Threads: @threads # Base.Threads.nthreads()

function RK4(f::Function, v::AbstractVector, h=1e-3)
    V1 = f(v)
    V2 = f(v + (h/2)*V1)
    V3 = f(v + (h/2)*V2)
    V4 = f(v + h*V3)
    return v + (h/6)*(V1 + 2V2 + 2V3 + V4), V1
end
function RK4(J::AbstractMatrix, U::AbstractMatrix, dt=1e-3)
    V1 = J*U
    V2 = J*(U + (dt/2)*V1)
    V3 = J*(U + (dt/2)*V2)
    V4 = J*(U + dt*V3)
    return U + (dt/6)*(V1 + 2V2 + 2V3 + V4)
end
function factory_lorenz(ρ::Number; ic = [10.,10.,10.], tspan = [0., 1000.], dt = 1e-3)
    σ = 10
    β = 8/3
    function sys(v::AbstractVector)
        x, y, z = v
        dx = σ*(y - x)
        dy = x*(ρ - z) - y
        dz = x*y - β*z
        return [dx, dy, dz]
    end
        
    t_ = first(tspan):dt:last(tspan)
    len_t_ = length(t_)
    
    t, tk = .0, 0
    v = ic; DIM = length(v)
    traj = zeros(len_t_+2, 2DIM)
    while tk ≤ len_t_
        v, dv = RK4(sys, v, dt)

        if t ≥ first(t_)
            tk += 1
            traj[tk+1,         1:DIM ] =  v
            traj[tk  , DIM .+ (1:DIM)] = dv
        end
    end
    return traj[2:(end-2), :]
end
factory_lorenz(T::Type, args...; kargs...) =
DataFrame(factory_lorenz(args...; kargs...), ["x", "y", "z", "dx", "dy", "dz"])

function gram_schmidt(J)
    N = size(J, 1)
    U, V = deepcopy(J), deepcopy(J)
    U[:,1] = V[:,1] / norm(V[:,1])
    for j in 2:N
        for jp in 1:(j-1)
            V[:,j] -= (J[:,j]'U[:,jp])*U[:,jp]
        end
        U[:,j] = V[:,j] / norm(V[:,j])
    end
    return U, V
end
J_lorenz(x,y,z,ρ) = [
     -10   10  0
      ρ-z -1  -x
        y  x  -8/3
]

const dt = 1e-3
ρ_range = 0:.1:100
cd("G:/DDM")

function lyapunov_lorenz_GS()
    schedules = DataFrame(idx = eachindex(ρ_range), ρ = ρ_range)
    schedules[!, :λ1] .= .0; schedules[!, :λ2] .= .0; schedules[!, :λ3] .= .0; 
    @showprogress @threads for dr = eachrow(schedules)
        data = factory_lorenz(DataFrame, dr.ρ; dt)
        
        λ = zeros(3);
        J = J_lorenz(collect(data[1, 1:3])..., dr.ρ)
        U = I(3)
        for i in 2:nrow(data)
            U, V = gram_schmidt(U)
            λ += V |> eachcol .|> norm .|> log          

            U = RK4(J, U, dt)
            J = J_lorenz(collect(data[i, 1:3])..., dr.ρ)
        end

        λ ./= dt*nrow(data)
        dr[[:λ1, :λ2, :λ3]] .= sort(λ, rev = true)
    end
    sort!(schedules, :ρ)
    CSV.write("lyapunov/lorenzGS.csv", schedules, bom = true)
    plot(xticks = 0:20:100, legend = :none, xlabel = "ρ", ylabel = "λ", size = [600, 300])
    plot!(schedules.ρ, schedules.λ1, lw = 2, color = 1)
    plot!(schedules.ρ, schedules.λ2, lw = 2, color = 2)
    plot!(schedules.ρ, schedules.λ3, lw = 2, color = 3)
    png("lyapunov/lorenzGS.png")
end

function lyapunov_lorenz_QR()
    schedules = DataFrame(idx = eachindex(ρ_range), ρ = ρ_range)
    schedules[!, :λ1] .= .0; schedules[!, :λ2] .= .0; schedules[!, :λ3] .= .0; 
    @showprogress @threads for dr = eachrow(schedules)
        data = factory_lorenz(DataFrame, dr.ρ; dt)
        
        λ = zeros(3);
        J = J_lorenz(collect(data[1, 1:3])..., dr.ρ)
        Q, _ = qr(J)
        for i in 2:nrow(data)
            Q, R = qr(Q)
            λ += R .|> abs |> diag .|> log   

            Q = RK4(J, Matrix(Q), dt)
            J = J_lorenz(collect(data[i, 1:3])..., dr.ρ)
        end

        λ ./= dt*nrow(data)
        dr[[:λ1, :λ2, :λ3]] .= sort(λ, rev = true)
    end
    sort!(schedules, :ρ)
    CSV.write("lyapunov/lorenzQR.csv", schedules, bom = true)
    plot(xticks = 0:20:100, legend = :none, xlabel = "ρ", ylabel = "λ", size = [600, 300])
    plot!(schedules.ρ, schedules.λ1, lw = 2, color = 1)
    plot!(schedules.ρ, schedules.λ2, lw = 2, color = 2)
    plot!(schedules.ρ, schedules.λ3, lw = 2, color = 3)
    png("lyapunov/lorenzQR.png")
end

lyapunov_lorenz_GS()
lyapunov_lorenz_QR()

See Also

• Lyapunov Exponents of Multidimensional Maps and Their Numerical Calculation