Calculating Lyapunov Exponents with ForwardDiff.jl and DifferentialEquations.jl

examples.ipynb

LyapunovExponentsWithForwardDiff.jl

module LyapunovExponentsWithForwardDiff

using DifferentialEquations
using ForwardDiff
using ParameterizedFunctions
using ProgressMeter
using RecipesBase

type LyapunovExponentsResult
    sol_p
    sol_pt
    values
end

Base.show(io::IO, res::LyapunovExponentsResult) =
    print(io, "Lyapunov Exponents: ", res.values[:, end])

type PhaseTangentParam
    phase_dynamics
    Jacobian
    function PhaseTangentParam(phase_dynamics, x0)
        new(phase_dynamics,
            similar(x0, (length(x0), length(x0))))
    end
end

function tangent_dynamics(t, u, param, du)
    ForwardDiff.jacobian!(
        param.Jacobian,
        (y, x) -> param.phase_dynamics(t, x, y),
        (@view du[:, 1]),       # phase space derivative goes here
        (@view u[:, 1]),
    )
    A_mul_B!((@view du[:, 2:end]), param.Jacobian, (@view u[:, 2:end]))
end

function keepgoing!(sol)
    sol.prob.u0 = sol(sol.prob.tspan[end])
    solve(sol.prob)
end

"""
S_n = ((n-1)/n) S_{n-1} + r_n / n
"""
@inline function lyap_add_R!(n, lyap, R)
    a = (n - 1) / n
    b = 1 - a
    for i = 1:length(lyap)
        lyap[i] = a * lyap[i] + b * log(R[i, i])
    end
end

""" A = A * diag(sgn) """
@inline function A_mul_sign!(A, sgn)
    for i = 1:size(A)[2]
        if !sgn[i]
            A[:, i] *= -1
        end
    end
    A
end

""" A = diag(sgn) * A """
@inline function sign_mul_A!(sgn, A)
    for i = 1:size(A)[1]
        if !sgn[i]
            A[i, :] *= -1
        end
    end
    A
end


""" sgn = sign(diag(A))  """
@inline function sign_diag!(sgn, A)
    for i = 1:size(A)[1]
        sgn[i] = A[i, i] >= 0
    end
    sgn
end

"""
`lyapunov_exponents(phase_dynamics, u0, t_chunk, num_tran, num_attr; ...)`

### Positional Arguments

* `phase_dynamics`: Definition of the phase space dynamics in the
  inplace `(t, u, du)` format (as in `ODEProblem`).
* `u0`: Initial state for the phase space dynamics.
* `t_chunk`: Length of numerical integration between orthonormalization.
* `num_tran`: Number of iterations to through away to get rid of the
  transient dynamics.
* `num_attr`: Number of orthonormalization steps for Lyapunov exponent
  calculation (which is presumably done inside an attractor).

### Keyword Arguments

* `dim_lyap`: Number of Lyapunov exponents to be calculated.
  Default to the full system dimension.
* `Q0`: The initial guess of the Gram-Schmidt "Lyapunov vectors".
  Default to the identity matrix.
"""
function lyapunov_exponents(
        phase_dynamics,
        u0,
        t_chunk,
        num_tran,
        num_attr;
        dim_lyap=length(u0),
        Q0=eye(length(u0), dim_lyap))

    # ODEProblem for dynamics in phase space
    pprob = ODEProblem(phase_dynamics, u0, (0.0, t_chunk))

    psol = solve(pprob)
    @showprogress 1 "Transient dynamics..." for _ in 2:num_tran
        psol = keepgoing!(psol)
    end

    # ODEProblem for dynamics in phase & tangent spaces
    pt0 = similar(u0, (length(u0), dim_lyap + 1))
    pt0[:, 1] = psol(psol.prob.tspan[end])  # phase space initial condition
    pt0[:, 2:end] = Q0                      # tangent space ...
    ptparam = PhaseTangentParam(phase_dynamics, u0)
    tprob = ODEProblem(
        ParameterizedFunction(tangent_dynamics, ptparam),
        pt0,
        (0.0, t_chunk),
    )

    # H2 method in Geist, Parlitz, Lauterborn (1990)
    lyap = similar(u0, dim_lyap)
    lehist = similar(u0, (dim_lyap, num_attr))
    signR = Array(Bool, dim_lyap)

    local tsol
    xt = tprob.u0[:, 1]
    Q = tprob.u0[:, 2:end]
    @showprogress 1 "Computing Lyapunov exponents..." for n in 1:num_attr
        tprob.u0[:, 1] = xt
        tprob.u0[:, 2:end] = Q
        tsol = solve(tprob)
        uend = tsol(tsol.prob.tspan[end])
        xt = uend[:, 1]
        P = uend[:, 2:end]

        F = qrfact!(P)
        Q = F[:Q][:, 1:dim_lyap]
        R = F[:R]

        sign_diag!(signR, R)        # signR = diagm(sign(diag(R)))
        A_mul_sign!(Q, signR)       # Q = Q * signR
        sign_mul_A!(signR, R)       # R = signR * R

        lyap_add_R!(n, lyap, R)
        lehist[:, n] = lyap
    end

    lehist /= t_chunk
    LyapunovExponentsResult(
        psol,
        tsol,
        lehist,
    )
end

""" Plot `LyapunovExponentsResult` via `RecipesBase`."""
@recipe function f(res::LyapunovExponentsResult;
                   known=nothing)
    dim_lyap = size(res.values)[1]
    layout --> (dim_lyap, 1)
    xscale --> :log10

    ylims = [[minimum(res.values[i, :]),
              maximum(res.values[i, :])]
             for i = 1:dim_lyap]
    if known != nothing
        for i = 1:dim_lyap
            ylims[i][1] = min(ylims[i][1], known[i])
            ylims[i][2] = max(ylims[i][2], known[i])
        end
    end
    for i = 1:dim_lyap
        ymin, ymax = ylims[i]
        dy = ymax - ymin
        ylims[i] = [ymin - dy * 0.05,
                    ymax + dy * 0.05]
    end

    for i in 1:dim_lyap
        @series begin
            subplot := i
            label --> ""
            ylabel := "LE$i"
            ylim --> ylims[i]
            res.values[i, :]
        end
        if known != nothing
            @series begin
                subplot := i
                linetype := :hline
                label --> ""

                # repeating ylabel/ylim; otherwise they are ignored
                ylabel := "LE$i"
                ylim --> ylims[i]

                [known[i]]
            end
        end
    end
end

""" Some example dynamical systems and their known Lyapunov exponents. """
module Examples

"""
Lorenz system.

* https://en.wikipedia.org/wiki/Lorenz_system
* http://sprott.physics.wisc.edu/chaos/comchaos.htm
* E. N. Lorenz, J. Atmos. Sci. 20, 130-141 (1963)
"""
function lorenz(t, u, du)
    du[1] = 10.0(u[2]-u[1])
    du[2] = u[1]*(28.0-u[3]) - u[2]
    du[3] = u[1]*u[2] - (8/3)*u[3]
end
lorenz_les_known = [0.9056, 0, -14.5723]

"""
Simplest piecewise linear dissipative chaotic flow.

* http://sprott.physics.wisc.edu/chaos/comchaos.htm
* S. J. Linz and J. C. Sprott, Phys. Lett. A 259, 240-245 (1999)
"""
function piecewise_linear(t, u, du)
    du[1] = u[2]
    du[2] = u[3]
    du[3] = -0.6 * u[3] - u[2] - (u[1] > 0 ? u[1] : -u[1]) + 1
end
piecewise_linear_les_known = [0.0362, 0, -0.6362]

end

end

Hey, thanks for the interest. I just noticed your comment.

Hmm... I have no idea. As you said, it's possible that the version of the packages I used could be old by now. FYI, I was using Julia 0.5 and Plots 0.10.3.