-
Notifications
You must be signed in to change notification settings - Fork 90
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Local lyapunov exponents map #29
Comments
I implemented the naive way using function compute_lyapunov(prob::ODEProblem; d0=1e-9, threshold=10^4*d0, dt = 0.1,
diff_eq_kwargs = Dict(:abstol=>d0, :reltol=>d0))
threshold <= d0 && throw(ArgumentError("Threshold must be bigger than d0!"))
if haskey(diff_eq_kwargs, :solver)
solver = diff_eq_kwargs[:solver]
else
println("Using DPRKN12 as default solver")
solver = DPRKN12()
end
# Initialize
st1 = prob.u0
integ1 = init(prob, solver; diff_eq_kwargs..., save_first=false, save_everystep=false)
integ1.opts.advance_to_tstop = true
prob.u0 = st1 .+ d0
integ2 = init(prob, solver; diff_eq_kwargs..., save_first=false, save_everystep=false)
integ2.opts.advance_to_tstop = true
prob.u0 = st1
lyapunov(integ1, integ2, prob.tspan[2]; d0=d0, threshold=threshold, dt=dt)
end Then I use the following function for the parallel λ computations function compute_λs(q0list, p0list, tmax, d0, dt, tr, kwargs, prefix)
n = size(q0list, 1) # number of initial conditions
λs = SharedArray{Float64}(n)
pmap(i->(prob = defProb(q0list[i,:], p0list[i,:], (0., tmax));
λs[i] = compute_lyapunov(prob, d0=d0, dt=dt, threshold=tr,
diff_eq_kwargs=kwargs)),
Progress(n),
1:n)
save("$prefix/lyapunov.jld", "λs", λs, "d0", d0, "dt", dt, "tr", tr,
"tmax", tmax, "n", n);
return λs
end in a function like this addprocs(4);
using JLD
using ArgParse
using Plots, LaTeXStrings
using ProgressMeter
using PmapProgressMeter
@everywhere begin
using OrdinaryDiffEq, DiffEqCallbacks
using ParallelDataTransfer
end
function main()
# Hamiltonian parameters
A, B, D = readdlm("param.dat")
E_list, tmax, d0, dt, tr, solver, kwargs, defProb = input_param()
# Broadcast parameters to all workers
sendto(workers(), A=A, B=B, D=D, defProb=defProb)
for E in E_list
sendto(workers(), E=E)
prefix = "../output/B$B D$D E$E"
if !isdir(prefix)
mkpath(prefix)
end
if isfile("$prefix/z0.jld")
q0list, p0list = load("$prefix/z0.jld", "q0list", "p0list")
else
error("$prefix/z0.jld not found! Generate the initial conditions.")
end
λs = compute_λs(q0list, p0list, tmax, d0, dt, tr, kwargs, prefix)
end
end where the initial conditions are read from a file. |
Thanks! There may be better approaches or algorithms I am not aware of, but on the very least your code gives a concrete start for a general parallelized interface in case a different approach is not found! |
This issue was moved to JuliaDynamics/ChaosTools.jl#5 |
I am talking about a function that would partition the phase space, and for each box of the phase space would give a value for a lyapunov exponent.
The naive way to do this is to run the already existing function
lyapunov
for each box of the phase-space, but maybe there are much better ways to do it. The above would be seriously slow, since the partitioning would have 3-6 orders of magnitude amount of initial conditions.@RainerEngelken may have some comments for this.
The text was updated successfully, but these errors were encountered: