maxLyapunov: Maximum lyapunov exponent

Description

Functions for estimating the maximal Lyapunov exponent of a dynamical system from 1-dimensional time series using Takens' vectors.

Usage

maxLyapunov(
  time.series,
  min.embedding.dim = 2,
  max.embedding.dim = min.embedding.dim,
  time.lag = 1,
  radius,
  theiler.window = 1,
  min.neighs = 5,
  min.ref.points = 500,
  max.time.steps = 10,
  number.boxes = NULL,
  sampling.period = 1,
  do.plot = TRUE,
  ...
)
# S3 method for maxLyapunov
divTime(x)
# S3 method for maxLyapunov
embeddingDims(x)
# S3 method for maxLyapunov
divergence(x)
# S3 method for maxLyapunov
plot(
  x,
  main = "Estimating maximal Lyapunov exponent",
  xlab = "time t",
  ylab = "S(t)",
  type = "p",
  col = NULL,
  pch = NULL,
  add.legend = T,
  ...
)
# S3 method for maxLyapunov
estimate(
  x,
  regression.range = NULL,
  do.plot = FALSE,
  use.embeddings = NULL,
  main = "Estimating maximal Lyapunov exponent",
  xlab = "time t",
  ylab = "S(t)",
  type = "p",
  col = NULL,
  pch = NULL,
  ylim = NULL,
  fit.col = NULL,
  fit.lty = 2,
  fit.lwd = 2,
  add.legend = T,
  ...
)

Value

A list with three components named $time$ and $s.function$. $time$ is a vector containing the temporal interval where the system evolves. It ranges from 0 to $max.time.steps \cdot sampling.period$. $s.function$ is a matrix containing the values of the $S(t)$ for each t in the time vector(the columns) and each embedding dimension (the rows).

The divTime function returns the time in which the divergence of close trajectories was computed.

The embeddingDims function returns the embeddings in which the divergence of close trajectories was computed

The divergence function returns the rate of divergence of close trajectories needed for the maximum Lyapunov exponent estimation.

In order to obtain an estimation of the Lyapunov exponent the user can use the estimate function. The $estimate$ function allows the user to obtain an estimation of the maximal Lyapunov exponent by averaging the slopes of the embedding dimensions specified in the use.embeddings parameter. The slopes are determined by performing a linear regression over the radius' range specified in regression.range

Arguments

time.series: The original time series from which the maximal Lyapunov exponent will be estimated.
min.embedding.dim: Integer denoting the minimum dimension in which we shall embed the time.series (see buildTakens).
max.embedding.dim: Integer denoting the maximum dimension in which we shall embed the time.series (see buildTakens).Thus, we shall estimate the Lyapunov exponent between min.embedding.dim and max.embedding.dim.
time.lag: Integer denoting the number of time steps that will be use to construct the Takens' vectors (see buildTakens).
radius: Maximum distance in which will look for nearby trajectories.
theiler.window: Integer denoting the Theiler window: Two Takens' vectors must be separated by more than theiler.window time steps in order to be considered neighbours. By using a Theiler window, we exclude temporally correlated vectors from our estimations.
min.neighs: Minimum number of neighbours that a Takens' vector must have to be considered a reference point.
min.ref.points: Number of reference points that the routine will try to use. The routine stops when it finds min.ref.points reference points, saving computation time.
max.time.steps: Integer denoting the number of time steps marking the end of the linear region.
number.boxes: Number of boxes that will be used in the box assisted algorithm (see neighbourSearch).
sampling.period: Sampling period of the time series. When dealing with a discrete system, the sampling.period should be set to 1.
do.plot: Logical value. If TRUE (default value), a plot of $S(t)$ Vs $t$ is shown.
...: Additional plotting parameters.
x: A maxLyapunov object.
main: A title for the plot.
xlab: A title for the x axis.
ylab: A title for the y axis.
type: Type of plot (see ?plot).
col: Vector of colors for each of the dimensions of the plot.
pch: Vector of symbols for each of the dimensions of the plot.
add.legend: add a legend to the plot?
regression.range: Vector with 2 components denoting the range where the function will perform linear regression.
use.embeddings: A numeric vector specifying which embedding dimensions should the estimate function use to compute the Lyapunov exponent.
ylim: Numeric vector of length 2, giving the y coordinates range.
fit.col: A vector of colors to plot the regression lines.
fit.lty: The type of line to plot the regression lines.
fit.lwd: The width of the line for the regression lines.

Author

Constantino A. Garcia

Details

It is a well-known fact that close trajectories diverge exponentially fast in a chaotic system. The averaged exponent that determines the divergence rate is called the Lyapunov exponent (usually denoted with $\lambda$). If $\delta(0)$ is the distance between two Takens' vectors in the embedding.dim-dimensional space, we expect that the distance after a time $t$ between the two trajectories arising from this two vectors fulfills: $$\delta (n) \sim \delta (0)\cdot exp(\lambda \cdot t)$$ The lyapunov exponent is estimated using the slope obtained by performing a linear regression of $S(t)=\lambda \cdot t \sim log(\delta (t)/\delta (0))$ on $t$. $S(t)$ will be estimated by averaging the divergence of several reference points.

The user should plot $S(t) Vs t$ when looking for the maximal lyapunov exponent and, if for some temporal range $S(t)$ shows a linear behaviour, its slope is an estimate of the maximal Lyapunov exponent per unit of time. The estimate routine allows the user to get always an estimate of the maximal Lyapunov exponent, but the user must check that there is a linear region in the $S(t) Vs t$. If such a region does not exist, the estimation should be discarded. The computations should be performed for several embedding dimensions in order to check that the Lyapunov exponent does not depend on the embedding dimension.

References

Eckmann, Jean-Pierre and Kamphorst, S Oliffson and Ruelle, David and Ciliberto, S and others. Liapunov exponents from time series. Physical Review A, 34-6, 4971--4979, (1986).

Rosenstein, Michael T and Collins, James J and De Luca, Carlo J.A practical method for calculating largest Lyapunov exponents from small data sets. Physica D: Nonlinear Phenomena, 65-1, 117--134, (1993).

Examples

Run this code

if (FALSE) {
## Henon System
h=henon(n.sample=  5000,n.transient= 100, a = 1.4, b = 0.3, 
        start = c(0.63954883, 0.04772637), do.plot = FALSE) 
my.ts=h$x 
ml=maxLyapunov(time.series=my.ts,
               min.embedding.dim=2,
               max.embedding.dim=5,
               time.lag=1,
               radius=0.001,theiler.window=4,
               min.neighs=2,min.ref.points=500,
               max.time.steps=40,do.plot=FALSE)
plot(ml)
ml.estimation = estimate(ml,regression.range = c(0,15),
                         use.embeddings=4:5,
                         do.plot = TRUE)
# The max Lyapunov exponent of the Henon system is 0.41
cat("expected: ",0.41," calculated: ",ml.estimation,"\n")

## Rossler system
r=rossler(a=0.15,b=0.2,w=10,start=c(0,0,0), time=seq(0,1000,0.1), 
          do.plot=FALSE)
my.ts=r$x
use.cols = c("#999999","#E69F00","#56B4E9")
ml=maxLyapunov(time.series=my.ts,min.embedding.dim=5,max.embedding.dim = 7,
               time.lag=12,radius=0.1,theiler.window=50,
               min.neighs=5,min.ref.points=length(r),
               max.time.steps=300,number.boxes=NULL,
               sampling.period=0.1,do.plot=TRUE,
               col=use.cols)
#  The max Lyapunov exponent of the Rossler system is 0.09
ml.est=estimate(ml,col=use.cols,do.plot=T,
                fit.lty=1,
                fit.lwd=5)
cat("expected: ",0.090," calculated: ",ml.est,"\n")
}

Run the code above in your browser using DataLab