Hybrid systems are a special type of dynamical system that switches between being continuous and discrete. An example could be the bouncing ball.

It will probably be a tiny bit tricky to implement them, but it has to be made as general as possible.

One need to have 2 equations of motion, one for continuous and one for discrete. Then a callback type function has to be created that will stop the ODE solving of DifferentialEquations and switch to discrete evolution (until some other kind of callback stops or for a specified amount of steps... implementation will show the best way...).

I tried the Lorenz example:

using DynamicalSystems, StaticArrays
ds = Systems.lorenz()
lyapunovs(ds,10^4;dt=1.0, Ttr=10.0)
3-element Array{Float64,1}:
   0.922575
  -1.25806 
 -13.3176

but the results are incorrect. The second Lyapunov exponent should be zero because that is the neutral exponent corresponding to perturbations in the direction of time.

There is SymEngine.jl,
SymPy and ParameterizedFunctions (which relies on SymEngine but defines a convenient interface for an ode problem and it's jacobian)

The lyapunov function adds a displacement to the current state to generate a neighboring state simply by adding a displacement like stat2 = state1 + d0 where d0 is a real number.

This can be bad in the case of e.g. Hamiltonian systems if the user is really concerned about energy conservation or where nearby states may have different energy that may lead to e,g, escaping to infinity instead of staying in the attractor. @SebastianM-C maybe you want to say something here.

I suggest the following simple fix: add an extra argument to lyapunov called displacement. This argument must be a function f(state, d0) where given a state and d0 produces a new state that is (approximately) d0 far away from the given state but in the way the user wants. The default value can simply be displacement = +.

To elaborate more on this, there is some discussion in the computation of lyapunov exponent of trying to create a neighboring state where the difference vector is already oriented along the maximal expansion direction. I don't really think this is necessary as orientation happens anyway and very quickly.

Makes sense or not?

Please follow the below guidelines, depending on whether you are
reporting an issue or you are requesting a new feature.

Issue description OR description of requested feature

mplah mplah.

Expected behavior OR description of possible implementation

mplah.

MWE (only for issue reports)

MWE stands for Minimal Working Example.

Please include a piece of code that can replicate your problem in a simple example. This will help us identify the problem much faster.

First, congratulations on how good DynamicalSystems is becoming. This work in Julia is amazing.

In nlts.jl the extension of reconstruction to Broomhead and King is really convincing and fast.
Since the pratitioner's starting point is always a dataset, it would be good to determine if the three times series from a given synthetic dataset (in your example gissinger + noise) come from either the same model or from three different models. At present, you apply separate SSA to each time series and then determine that tau=30 in each case. Can you extend SSA to Multichannel SSA https://arxiv.org/pdf/1309.5050 to complete the demonstration and bring to Julia a diagnostic tool already available in R.

Some kind of interface has to be constructed for using Hamiltonian continuous systems and integrating them without changing energy.

Either as a callback or as splitting the system to 2 equations and using dynamical odes

dynamical odes:
http://docs.juliadiffeq.org/stable/types/dynamical_types.html
http://docs.juliadiffeq.org/stable/solvers/dynamical_solve.html

callbacks:
http://docs.juliadiffeq.org/stable/features/callback_library.html

I wondered why it's g in the kwargs but G in the function body?

function double_pendulum(u0=rand(4); g=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
...
du[2] = (M2*L1*state[2]*state[2]*sin(del_)*cos(del_) +
               M2*G*sin(state[3])*cos(del_) +
               M2*L2*state[4]*state[4]*sin(del_) -
               (M1 + M2)*G*sin(state[1]))/den1
...

For both the lyapunovs and lyapunov.

add note that users of this package are expected to have basic knowledge of nonlinear dynamics and chaos.

Put up 3 books as well: strogatz, ott and argyris.

The plotting functions are not loaded with PyPlot with the current Require.

Maybe the even fit much better in this repo... Consider either moving them here or ask Require to make it possible to work at the current state.

ATM you have to explicitly write using ChaosTools to get e.g. plot_linear_regions.

Just trying to run the nlts example for continuous time. The command:

using DynamicalSystems, PyPlot

ds = Systems.lorenz() # Max lyapunov is around 0.90
dt = 0.05
x = timeseries(ds, 1000.0; dt = dt)[:, 1]
τ1 = estimate_delay(x) #gives 7

Returns error:

UndefVarError: x not defined

Stacktrace:
[1] estimate_delay(::Array{Float64,1}) at C:\Users\Denis.julia\v0.6\DynamicalSystems\src\delay_coords.jl:186
[2] include_string(::String, ::String) at .\loading.jl:515

The methods that calculate lyapunovs of continuous systems are extremely slow.
So far I am not sure whether this is due to DifferentialEquations, the qr-decomposition or that I Have simply written bad code. Profiling them will reveal issues in my code (which probably sucks).

In any case, the have to be made faster!

One way that this can be improved is by taking advantage of the integrator and step! interfaces of DiffEq.

For the moment, I would suggest to concentrate on using Automatic Differentiation, rather than hand-coded Jacobians or symbolic.

Do you have examples where you have checked that hand-coded Jacobians are actually much faster than ForwardDiff.jl?

cc @jrevels

Cannot tag a new version "v0.1.0" preceding all existing versions.
cc: @Datseris

Ott (et al) came up with a definitive way to claim whether or not a system is chaotic. This was done for the 25th anniversary of the journal Chaos.

In the article Defining Chaos he and Brian R. Hunt define a quantity called "expansion entropy" which can characterize chaoticity.

It would be cool to have a function that computes that! The paper also includes a numerical algorithm to compute the quantity.

StaticArrays recently got a PR that implements QR decomposition for SMatrix.

JuliaArrays/StaticArrays.jl#328

We have to use it here and replace the Base.qr whenever possible (in all discrete system algorithms that use QR), because of super juice benchmarks,. (see the PR)

StaticArrays.jl is great, but not useful for large arrays that are associated with high-dimensional systems, Maybe one could use StaticArrays only for low dimensional systems (say up to 10 variables) and use usual Base.Arrays for high dimensional systems? Also, I guess it might be worth mentioning that the computational cost involved in direct estimating entropies grows exponential in the degrees of freedom, so this is also rather limited to low dimensional systems.

Please follow the below guidelines, depending on whether you are
reporting an issue or you are requesting a new feature.

Issue description OR description of requested feature

There is a missing DOI for one of the papers listed in your bibliography; I believe the appropriate one can be found here

Cannot tag a new version "v0.1.0" preceding all existing versions.
cc: @Datseris

Currently the introduction has only one sentence about what this package is about and then says many stuff about gitter and whatever.

The first thing that should be displayed is useful information!

Another commonly used entropy is the "KS-entropy".

I would be reluctant to point people to the original papers (Ya.G. Sinai, On the Notion of Entropy of a Dynamical System, Doklady of Russian Academy of Sciences, (1959), 124, 768-771. and A.N. Kolmogorov, Entropy per unit time as a metric invariant of automorphism, Doklady of Russian Academy of Sciences, (1959), 124, 754-755.) as they are veeery mathematical and were written before computers became a thing. I also do not know any paper that treats the subject (I haven't searched for one!).

I think @RainerEngelken has done this, maybe he can help us a bit here with some comments.

https://github.com/Datseris/DynamicalSystems.jl/blob/master/src/continuous.jl#L79

Why not have that just extend ODEProblem to add a constructor for your input type?

We use ForwardDiff for discrete systems, but I am very sloppy and do not do the pre-allocations that are described here:

http://www.juliadiff.org/ForwardDiff.jl/stable/user/api.html#Preallocating/Configuring-Work-Buffers-1

Doing them will increase the speed of systems with automatically computed Jacobians.

With the latest release of DynamicalSystems I think it is about time that this package as well DynamicalBilliards go into an organization. I just learned that there is an organization called JuliaPhysics. I personally feel like this general name actually suit these packages, because they are also supposed to be general, even though they are specialized. Plus, you don't want to flood the "market" with organizations.

I would like to ask others that have contributed here or at the DynamicalBilliards, how they feel like about it. I am tagging you, please let me know what you think.
@RainerEngelken
@SebastianM-C
@dapias
@ChrisRackauckas
@dpsanders
@dapias

Similarly how it is done for lyapunov exponents here : https://juliadynamics.github.io/DynamicalSystems.jl/latest/nlts/#example-of-numerical-lyapunov-computation

do the same in the preceding section but for the dimension information_dim(R) as well.

Also put a figure (link it through imgur).

The Smaller and Generalized Alignment Indices methods are robust techniques that can distinguish chaotic behavior from regular, or stable periodic orbits from unstable ones, efficiently.

These efficient chaos indicators are based on the evolution of two or more, initially distinct, deviation vectors from a studied orbit. The methods are in general faster than computing maximum Lyaopunov exponents

In my eyes, the best literature to get this method from is the book "Chaos Detection and Predictability", Lecture Notes In Physics 915, Springer (Chapter 5).

For reference, these are the original publications related to the Issue:
(The following are to the best of my knowledge):

The first publication that introduced the Alingment Indices method is https://www.google.de/search?q=Skokos+Ch+2001+J.+Phys.+A:+Math.+Gen.+34+10029&ie=utf-8&oe=utf-8&gws_rd=cr&dcr=0&ei=AGnSWcHzFcqKgAbY3rCoDg

However, most known it became by this paper: http://iopscience.iop.org/article/10.1088/0305-4470/37/24/006/meta

Whereas the most known GALI paper is: http://www.sciencedirect.com/science/article/pii/S0167278907001273

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.

I am guessing that this may fall to all the steps taken in the automatical procedure to produce a dimension based on the scaling of genentropy. If yes, maybe the default keyword arguments of all calls have to be revised.

In system_definiion.jl, running:

@inline @inbounds function eom_roessler!(du, u)
a = 0.2; b = 0.2; c = 5.7
du[1] = -u[2]-u[3]
du[2] = u[1] + au[2]
du[3] = b + u[3](u[1] - c)
end
@inline @inbounds function jacob_roessler(u)
i = one(eltype(u))
o = zero(eltype(u))
@smatrix [o -i -i;
i a o;
u[3] o u[1] - c]
end

ros = ContinuousDS(rand(3), eom_roessler!, jacob_roessler)

produces this error:

UndefVarError: ds not defined

Stacktrace:
[1] show(::IOContext{Base.AbstractIOBuffer{Array{UInt8,1}}}, ::DynamicalSystems.ContinuousDS{Array{Float64,1},#eom_roessler!,#jacob_roessler}) at C:\Users\Denis.julia\v0.6\DynamicalSystems\src\systems\continuous.jl:139
[2] limitstringmime(::MIME{Symbol("text/plain")}, ::DynamicalSystems.ContinuousDS{Array{Float64,1},#eom_roessler!,#jacob_roessler}) at C:\Users\Denis.julia\v0.6\IJulia\src\inline.jl:24
[3] display_dict(::DynamicalSystems.ContinuousDS{Array{Float64,1},#eom_roessler!,#jacob_roessler}) at C:\Users\Denis.julia\v0.6\IJulia\src\execute_request.jl:27
[4] execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\Denis.julia\v0.6\IJulia\src\execute_request.jl:182
[5] eventloop(::ZMQ.Socket) at C:\Users\Denis.julia\v0.6\IJulia\src\eventloop.jl:8
[6] (::IJulia.##14#17)() at .\task.jl:335

There are two ways to create an e.g. ContinuousDS. One is to use closures, like it is done for almost all systems:

function roessler(u0=rand(3); a = 0.2, b = 0.2, c = 5.7)
    @inline @inbounds function eom_roessler!(du, u)
        du[1] = -u[2]-u[3]
        du[2] = u[1] + a*u[2]
        du[3] = b + u[3]*(u[1] - c)
    end
    i = one(eltype(u0))
    o = zero(eltype(u0))
    J = zeros(eltype(u0), 3, 3)
    J[1,:] .= [o, -i,      -i]
    J[2,:] .= [i,  a,       o]
    J[3,:] .= [u0[3], o, u0[1] - c]
    @inline @inbounds function jacob_roessler!(J, u)
        J[3, 1] = u[3]; J[3,3] = u[1] - c
    end
    name = "Roessler76 system (a=$(a), b=$(b), c=$(c))"
  return ContinuousDS(u0, eom_roessler!, jacob_roessler!, J; name = name)
end

The thing with this approach is that when one wants to change parameters, one needs to create a completely new instance of ContinuousDS because the enclosing function has to be called:
ds = rossler(; a = 5). Maybe this is extremely inefficient.

A second way to do it is to make a Functor, e.g.:

truct Lorenz96{T <: Real} # Structure with the parameters of Lorenz96 system
  F::T
  N::Int
end
# Equations of motion 
@inbounds function (obj::Lorenz96{T})(dx, x) where T
  N, F = obj.N, obj.F
  # 3 edge cases
  dx[1] = (x[2] - x[N - 1]) * x[N] - x[1] + F
  dx[2] = (x[3] - x[N]) * x[1] - x[2] + F
  dx[N] = (x[1] - x[N - 2]) * x[N - 1] - x[N] + F
  # then the general case
  for n in 3:(N - 1)
    dx[n] = (x[n + 1] - x[n - 2]) * x[n - 1] - x[n] + F
  end
  return dx
end

"""
    lorenz96(N::Int, u0 = rand(M); F=0.01)
`N` is the chain length, `F` the forcing. Jacobian is created automatically.
"""
function lorenz96(N::Int, u0 = rand(N); F=0.01)
    name = "Lorenz96 system, chain of $N (F = $(F))"
    lor96 = Lorenz96(F, N) # create struct
    return ContinuousDS(u0, lor96; name = name)
end

Becauses Lorenz96 is a mutable type, after the first creation of ds = lorenz96(5; F = 0.2)
one could then simply do ds.eom!.F = 23.4 to change the F parameter, without needing to re-create a dynamical system instance.

In a recent update I have made all DynamicalSystem types to be immutable instead. But still, maybe it is worth it to make this change and use functors everywhere?

Maybe @ChrisRackauckas has something to add here? Have you had this dilemma somewhere in DiffEq where you wondered "closure vs functor" ? I know that Julia translates all closures into functors but that pretty much creates a new instance of the functor each time.

If a Dataset has data that is not SVectors, it is not passable to KDTree.

This has to be fixed of course. Either you make the dataset be only with SVectors, or you make it SVectors
before passing.

Of course, making the data a matrix works as well.

This is extremely simple:

The only changes necessary are:

1. in the methods add an optional argument (NOT a keyword) e.g. k::Int = dimension(ds)
2. initialize the matrices Q and K in the lyapunovs function to be eye(k)
3. For continuous systems make the function tangentbundle_setup also depend on k and initialize properly and define the e.o.m. properly
4. update the doc string
5. Pull request and BOOM! You contributed!

The renyi function that calculates the renyi information (or entropy) uses the Histogram function of StatsBase.

However, only the non-zero elements are required and also, there is no need to store the edges of the bins.
(see source code, function: non0hist).

A dedicated function that does not create the complicated Histogram (by calling fit) would be much faster.

This should be added to all functions, from evolve and timeseries to lyapunovs etc..

Also, it should be done in two ways: first for Juno using the innate progress bar that has been set up,
and also for the REPL, using ProgressMeter.jl.

I suggest to simply create a progress(n, N, message) function that displays the message and the progress based on n which should go up to N. This function should know/be configured to decide whether to use Juno or ProgressMeter.

After this is done it is easy to put it on all functions!

This will (probably) help other developers

Running lyapunovs.jl at:

using DynamicalSystems

lor = Systems.lorenz(ρ = 32.0) #this is not the original parameter!
issubtype(typeof(ds), ContinuousDynamicalSystem) # true

λλ = lyapunovs(lor, 10000,
dt = 0.1, diff_eq_kwargs = Dict(:abstol => 1e-9, :reltol => 1e-9))
 result:
[0.999176, 0.000774754, -14.6666]

gives the following error:

UndefVarError: ContinuousDynamicalSystem not defined

Stacktrace:
 [1] include_string(::String, ::String) at .\loading.jl:515

Somewhere in the introduction we should point out that we use 3 data structurs (and various derivatives) in our packages.

1. continuous systems, that use DiffEq
2. discrete systems.
3. numerical datasets

Hi,

Excited to find this package for Julia, however I seem to run into small errors when trying to use the package with the included logistic function, and finding of Lyapunov Exponents or Max Lyapunov Exponents (with lyapunovs and lyapunov respectively). These seem to work for all other default function I have tried.

This is quite simple to explain, so please excuse me for not using your error reporting format exactly.

Using such simple examples as below, the answer is returned perfectly fine:

using DynamicalSystems
ds = Systems.roessler()
λ = lyapunov(ds, 10000)

Or
ds = Systems.towel()
λ = lyapunov(ds, 10000)

Every default function I've tried seems to work with default conditions in this manner. However using the logistic function for example:

ds = Systems.logistic()
λ = lyapunov(ds, 10000)

Returns the error:
DimensionMismatch("Cannot multiply two vectors")

I do not understand what causes the difference, similarly errors are encountered when trying to use lyapunovs with the logistic function which are not encountered with other default functions. (Experimenting with parameters does not seem to help.)

Plot of the orbit diagram for the logistic function using the example code on the other hand works perfectly fine.

With the new release of JUNO the console documentation will allow display of math equations.

All systems should have their equations of motion defined in the documentation string, right under the call signature.
Eg.
```julia
lorenz(a,b,c) = ...
```
```math
\dot{x} = ...\
\dot{y} = ...
```
etc.

The documentation exporting files should also be configured so that these equations are also displayed in the official documentation page. So far I have been unsuccsesful in making it work.

See e.g.: https://qojulia.org/

I am super-jelly that the link given looks so sexy.

Can we also do something similar for DynamicalSystems.jl ?

This is only documentation, no new feature or anything...
See https://juliadynamics.github.io/DynamicalSystems.jl/latest/nlts/#case-of-a-continuous-system

Just use a timeseries from the Lorenz attractor :P

Documentation on: https://juliadynamics.github.io/DynamicalSystems.jl/latest/definition/general/
Tree expansion: Fundamentals -> Dynamical Systems -> Section: Equations of motion

The labels "iip" and "oop" in the bulleted list seem to be mis-matched with iip returning an SVector for xnew and oop changing xnew by reference, which is inconsistent with the summary and examples.

From here: docs.juliadiffeq.org/latest/analysis/uncertainty_quantification.html#Example-3:-Adaptive-ProbInts-on-the-Lorenz-Attractor-1

OrdinaryDiffEq got a new reinit!(integrator) interface that allows one to easily change e.g. the state of an integrator.

SciML/OrdinaryDiffEq.jl#218

We should use it whenever possible as it will probably improve speed.

In periodicity.jl, running:

using DynamicalSystems, PyPlot

ds = Systems.standardmap()
xs = linspace(0, 2π, 21); ys = copy(xs)
ics = [SVector{2}(x,y) for x in xs for y in ys]

gives the error:

UndefVarError: SVector not defined

Stacktrace:
[1] next at .\iterators.jl:708 [inlined]
[2] grow_to!(::Array{Union{},1}, ::Base.Iterators.Flatten{Base.Generator{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},##1#3}}, ::Tuple{Int64,Base.Generator{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},##2#4{Float64}},Int64}) at .\array.jl:498
[3] grow_to!(::Array{Any,1}, ::Base.Iterators.Flatten{Base.Generator{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},##1#3}}) at .\array.jl:491
[4] collect(::Base.Iterators.Flatten{Base.Generator{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},##1#3}}) at .\array.jl:397
[5] include_string(::String, ::String) at .\loading.jl:515

JuliaLang/Pkg.jl#176 (comment)

this is a dummy issue to see if RocketChat integration works!