Time Evolution of Systems

Trajectory and Timeseries

The word "timeseries" can be confusing, because it can mean a univariate (also called scalar or one-dimensional) timeseries or a multivariate (also called multi-dimensional) timeseries. To resolve this confusion, in DynamicalSystems.jl we have the following convention: "timeseries" always refers to a one-dimensional vector of numbers, which exists with respect to some other one-dimensional vector of numbers that corresponds to a time-vector. On the other hand, the word "trajectory" is used to refer to a multi-dimensional timeseries, which is of course simply a group/set of one-dimensional timeseries.

Note that the data representation of a "trajectory" in Julia may vary: from a 2D Matrix to independent Vectors. In our package, a trajectory is always represented using a Dataset, which is a Vector of SVectors, and each SVector represents a data-point (the values of the variables at a given time-point).

DynamicalSystems.jl provides a convenient function for getting a trajectory of a system at equally spaced time points:

DynamicalSystemsBase.trajectory — Function

trajectory(ds::DynamicalSystem, T [, u]; kwargs...) -> dataset

Return a dataset that will contain the trajectory of the system, after evolving it for total time T, optionally starting from state u. See Dataset for info on how to use this object.

A W×D dataset is returned, with W = length(t0:dt:T) with t0:dt:T representing the time vector (not returned) and D the system dimension. For discrete systems both T and dt must be integers.

Keyword Arguments

dt : Time step of value output during the solving of the continuous system. For discrete systems it must be an integer. Defaults to 0.01 for continuous and 1 for discrete.
Ttr : Transient time to evolve the initial state before starting saving states.
diffeq... : Keyword arguments propagated into init of DifferentialEquations.jl. For example abstol = 1e-9. Only valid for continuous systems. If you want to specify a solver, do so by using the name alg, e.g.: alg = Tsit5(), maxiters = 1000. This requires you to have been first using OrdinaryDiffEq to access the solvers. See DynamicalSystemsBase.CDS_KWARGS for default values. These keywords can also include callback for event handling. Using a SavingCallback with trajectory will lead to unexpected behavior!

Notice that if you want to do repeated evolutions of different states of a continuous system, you should use the integrator interface instead.

Example

using DynamicalSystems
ds = Systems.towel()
tr = trajectory(ds, 100)

3-dimensional Dataset{Float64} with 101 points
 0.085     -0.121       0.075   
 0.285813  -0.0675286   0.238038
 0.76827   -0.038933    0.672094
 0.681871   0.0508933   0.825263
 0.837347   0.0372633   0.555269
 0.51969    0.0616256   0.940906
 0.966676  -0.00171595  0.2225  
 0.112748   0.0674955   0.653573
 0.386547  -0.0886542   0.869349
 0.910741  -0.0316828   0.411607
 ⋮                              
 0.729743   0.0566836   0.25032 
 0.739275   0.0308018   0.720691
 0.740845   0.047263    0.767058
 0.740185   0.0494093   0.684863
 0.738165   0.0466359   0.825703
 0.747373   0.0506512   0.553337
 0.719603   0.0437958   0.944377
 0.784243   0.0495777   0.20732 
 0.631288   0.0375437   0.631114

To get every 3-rd point of the trajectory, do

tr = trajectory(ds, 100; dt = 3)

3-dimensional Dataset{Float64} with 34 points
 0.085     -0.121       0.075   
 0.681871   0.0508933   0.825263
 0.966676  -0.00171595  0.2225  
 0.910741  -0.0316828   0.411607
 0.545332   0.0508239   0.819404
 0.544182  -0.0940112   0.270849
 0.690701  -0.0679162   0.727605
 0.939111   0.010943    0.736301
 0.841034   0.0322099   0.689976
 0.140788   0.0876455   0.918065
 ⋮                              
 0.223785   0.0902301   0.75948 
 0.510326   0.0726233   0.682717
 0.529268  -0.069801    0.962927
 0.435893  -0.0686173   0.92088 
 0.583619  -0.0769684   0.76526 
 0.729705  -0.0680752   0.508933
 0.739275   0.0308018   0.720691
 0.738165   0.0466359   0.825703
 0.784243   0.0495777   0.20732

Identical syntax is used for continuous systems

ds = Systems.lorenz()
tr = trajectory(ds, 10.0; dt = 0.01)

3-dimensional Dataset{Float64} with 1001 points
  0.0       10.0       0.0      
  0.951226  10.0352    0.0479002
  1.82768   10.3224    0.18683  
  2.6593    10.839     0.416592 
  3.47133   11.5678    0.744804 
  4.28518   12.496     1.18585  
  5.1191    13.6125    1.76047  
  5.98858   14.9063    2.49569  
  6.90648   16.3638    3.42487  
  7.88291   17.9659    4.58749  
  ⋮                             
 -8.3603    -9.01682  25.9478   
 -8.42664   -9.09574  26.0152   
 -8.49391   -9.16911  26.0931   
 -8.56145   -9.23601  26.1809   
 -8.62859   -9.29558  26.278    
 -8.69462   -9.34701  26.3833   
 -8.75883   -9.38954  26.4959   
 -8.82053   -9.4225   26.6146   
 -8.879     -9.44532  26.7382

And a final example controlling the integrator accuracy:

ds = Systems.lorenz()
tr = trajectory(ds, 10.0; dt = 0.1, atol = 1e-9, rtol = 1e-9)

3-dimensional Dataset{Float64} with 101 points
  0.0      10.0       0.0    
  8.92478  19.6843    6.02839
 20.0553   24.8278   39.9785 
  9.89834  -7.67213  42.312  
 -2.67297  -9.86486  30.4531 
 -6.73013  -8.71064  27.4841 
 -7.91     -8.6127   26.6731 
 -8.44188  -8.84913  26.7275 
 -8.72491  -8.86503  27.1909 
 -8.70356  -8.53165  27.5098 
  ⋮                          
 -7.90031  -8.12215  25.9823 
 -8.35258  -8.96853  25.9951 
 -8.94423  -9.40349  26.9656 
 -9.09671  -8.90393  27.9988 
 -8.62374  -7.96528  28.0069 
 -8.01075  -7.54849  27.038  
 -7.8329   -7.96259  26.0287 
 -8.23296  -8.84597  25.8467 
 -8.879    -9.44531  26.7382

Solution precision for continuous systems

A numerical solution of an ODE is not the "true" solution, uniquely defined by a (well-defined) ODE and an initial condition. Especially for chaotic systems, where deviations are amplified exponentially, one is left worried if the numerical solutions truly are part of the system and can truly give insight in understanding the system.

DifferentialEquations.jl offers a tool, called Uncertainty Quantification, which allows users to asses up to what time-scales the numerical solution is close to the "true" solution. For example, using the default solving parameters of DynamicalSystems.jl, the Lorenz system is accurate up to time t = 50.0.

However, fortunately for us, there is not too much worry about the numerical solution diverging from the true solution. That is because of the shadowing theorem (or shadowing lemma):

Shadowing Theorem

Although a numerically computed chaotic trajectory diverges exponentially from the true trajectory with the same initial coordinates, there exists an errorless trajectory with a slightly different initial condition that stays near ("shadows") the numerically computed one.

This simply means that one can always numerically study chaos not only qualitatively but also quantitatively. For more information, see the book Chaos in Dynamical Systems by E. Ott, or the scholarpedia entry.