Hi all. I have been working on a package on my free time that I have mentioned in a previous post in this category. Now it is at a decent state and has been accepted as a Julia package. It’s called DynamicalSystems.jl! With this post I would like to do the following:
- Introduce the package
- Hear constructive criticism about what could be done better regarding the internals
(mainly, should all data be represented as a vector of static arrays?) - Invite relevant groups to create JuliaDynamics organization with me!
Point 1
The main goal is for it to be a useful companion for students and scientists treading on the field of Chaos, nonlinear dynamics and dynamical systems in general. I am really hopping that it will expand to become something big, assuming enough people find it useful and in turn would like to contribute. Therefore, if you feel like it, don’t be shy! 
So far the package is very basic however it does offer a lot of fundamental stuff:
- Interface for defining dynamical systems and smoothly interacting with DifferentialEquations (for continuous systems)
- Computation of lyapunov spectra
- Computation of entropies in a highly efficient manner
- Automated computation of attractor dimension (by finding the linear region automatically)
- Delay coordinates reconstruction
- Very detailed documentation strings and in general documentations (I actually cite papers where the method is taken from because one of the goals for the package is to be as transparent as possible)
**If you want to test the package right now, please use checkout: Pkg.add("DynamicalSystems"); Pkg.checkout("DynamicalSystems") because there have been some major changes in automated attractor dimension estimation.
For example, to get the lyapunov spectrum of a map, and then the attractor dimension you would simply do:
using DynamicalSystems, StaticArrays
@inline @inbounds eom_henon(x) = SVector{2}(1.0 - a*x[1]^2 + x[2], b*x[1])
@inline @inbounds jacob_henon(x) = @SMatrix [-2*a*x[1] 1.0; b 0.0]
he = DiscreteDS(zeros(2), eom_henon, jacob_henon)
this defines a system, and if the Jacobian is missing, it is calculated automatically using ForwardDiff. Now,
ly = lyapunovs(ds, 100000)
# result: [0.420, -1.62]
gives the lyapunov spectrum by applying QR method 100000 times. Also, you can see that it is not very slow:
using BenchmarkTools
@btime lyapunovs(he, 100000);
# result: 298.302 ms (2100006 allocations: 119.02 MiB)
but can probably be even better. Moving on, it is also very simple to calculate generalized entropies and thus dimensions.
ts = timeseries(he, 100000) # a 100000 x 2 matrix that contains the data
α = 1; #generalized entropy exponent
ε = 0.001 # boxsize
g = genentropy(α, ε, ts) #generalized entropy
D = generalized_dim(α, ts) #generalized dimension (aka Renyi dimension)
# Result: D = 1.2291586862831074
@btime generalized_dim(1, $ts);
# Result: 1.701 s (25138684 allocations: 522.62 MiB)
If you ever use the function generalized_dim() please give me feedback on how well it does, because it performs a lot of automated steps (see the documentation string). In essence it calculates the generalized entropy for many different boxsizes.
All these functions that I have used here are generic and will work for your systems as well, not only the Henon map!
Point 2
Data inside the package is represented as an N x D matrix with N the amount of data points and D the system dimension. There have been comments that the better off approach would be to completely change this and make all data be represented as a vector of StaticArrays. This would not have much of an impact on Discrete systems, since they already operate with static arrays, but it will be a significant change for Continuous systems that operate using in-place functions.
The “best” interface for data representation could be to use RecursiveArrayTools and have the main data structure be a recursive array of static arrays. (notice that for my package almost all data is in the form of timeseries)
@ChrisRackauckas and @kristoffer.carlsson had some comments on that. I would seriously appreciate it if you had a chance to only give a glance at the source code (see e.g. the function timeseries in the files discrete.jl and continuous.jl) and simply confirm if you feel this is best.
Besides this, if anybody that uses the package has ideas for improvement please go ahead and say them, or even better contribute them directly.
Point 3
I have made 2 packages so far; this one and DynamicalBilliards.jl. Both are about some kind of dynamical systems. I am sure there are other packages out there that fit this theme. Interest was already shown from my previous post in this topic. Some of the interested persons were @tkoolen with his package about rigid body dynamics.
However, there isn’t yet any Julia organization that unites these packages. I want to change that, and make this organization. However, I would really, really prefer to not start it by myself with only my 2 packages. So, if there is anybody that wants to participate in this, please let me know asap because I will probably do this organization thingy before the end of september.
