Skip to content

Instantly share code, notes, and snippets.

@tkf
Last active June 3, 2019 01:06
Show Gist options
  • Save tkf/014787f386e99b7b8b42683cfbd8da01 to your computer and use it in GitHub Desktop.
Save tkf/014787f386e99b7b8b42683cfbd8da01 to your computer and use it in GitHub Desktop.
Calculating Lyapunov Exponents with ForwardDiff.jl and DifferentialEquations.jl
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
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
@finmod
Copy link

finmod commented May 10, 2017

Very nice notebook on the use of DifferentialEquations and ForwardDiff. I have a slight syntax problem with the two:

plot(sol, vars=(1,2,3), linewidth=0.5, label="")

Must be simple but Plots is constantly under revision. I get this error:

MethodError: Cannot convert an object of type Tuple{Float64} to an object of type Tuple{}
This may have arisen from a call to the constructor Tuple{}(...),
since type constructors fall back to convert methods.

in macro expansion at .\broadcast.jl:129 [inlined]
in macro expansion at .\simdloop.jl:73 [inlined]
in macro expansion at .\broadcast.jl:123 [inlined]
in _broadcast!(::DiffEqBase.##71#72, ::Array{Tuple{Float64,Float64},1}, ::Tuple{Tuple{Bool},Tuple{Bool},Tuple{Bool}}, ::Tuple{Tuple{Int64},Tuple{Int64},Tuple{Int64}}, ::Tuple{Array{Float64,1},Array{Float64,1},Array{Float64,1}}, ::Type{Val{3}}) at .\broadcast.jl:117
in broadcast!(::Function, ::Array{Tuple{Float64,Float64},1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}) at .\broadcast.jl:172
in broadcast_t(::Function, ::Type{T}, ::Array{Float64,1}, ::Vararg{Array{Float64,1},N}) at .\broadcast.jl:228
in broadcast(::Function, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Vararg{Array{Float64,1},N}) at .\broadcast.jl:230
in solplot_vecs_and_labels(::Int64, ::Array{Tuple{DiffEqBase.##71#72,Int64,Int64,Int64},1}, ::Array{Array{Float64,1},1}, ::Array{Float64,1}, ::DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,#lorenz,Void,UniformScaling{Int64}},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#lorenz,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}}, ::Bool, ::Void) at C:\Users\Denis.julia\v0.5\DiffEqBase\src\solutions\solution_interface.jl:300
in macro expansion at C:\Users\Denis.julia\v0.5\DiffEqBase\src\solutions\solution_interface.jl:94 [inlined]
in apply_recipe(::Dict{Symbol,Any}, ::DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,#lorenz,Void,UniformScaling{Int64}},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#lorenz,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}}) at C:\Users\Denis.julia\v0.5\RecipesBase\src\RecipesBase.jl:238
in _process_userrecipes(::Plots.Plot{Plots.PyPlotBackend}, ::Dict{Symbol,Any}, ::Tuple{DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,#lorenz,Void,UniformScaling{Int64}},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#lorenz,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}}}) at C:\Users\Denis.julia\v0.5\Plots\src\pipeline.jl:73
in _plot!(::Plots.Plot{Plots.PyPlotBackend}, ::Dict{Symbol,Any}, ::Tuple{DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,#lorenz,Void,UniformScaling{Int64}},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#lorenz,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}}}) at C:\Users\Denis.julia\v0.5\Plots\src\plot.jl:171
in (::Plots.#kw##plot)(::Array{Any,1}, ::Plots.#plot, ::DiffEqBase.ODESolution{Float64,2,Array{Array{Float64,1},1},Void,Void,Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,#lorenz,Void,UniformScaling{Int64}},OrdinaryDiffEq.Tsit5,OrdinaryDiffEq.InterpolationData{#lorenz,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}}}) at .<missing>:0

@tkf
Copy link
Author

tkf commented May 28, 2017

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment