Nothing
#' Maximum lyapunov exponent
#' @description
#' Functions for estimating the maximal Lyapunov exponent of a dynamical system
#' from 1-dimensional time series using Takens' vectors.
#' @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
#' \eqn{\lambda}{lambda}). If \eqn{\delta(0)}{delta(0)} is the distance between
#' two Takens' vectors in the embedding.dim-dimensional space, we expect that
#' the distance after a time \eqn{t} between the two trajectories arising
#' from this two vectors fulfills:
#' \deqn{\delta (n) \sim \delta (0)\cdot exp(\lambda \cdot t)}{\delta (n) is.approximately \delta (0) exp(\lambda *t).}
#' The lyapunov exponent is estimated using the slope obtained by performing a
#' linear regression of
#' \eqn{S(t)=\lambda \cdot t \sim log(\delta (t)/\delta (0))}{S(t)=\lambda *t is.approximately log(\delta (t)/\delta (0))}
#' on \eqn{t}. \eqn{S(t)} will be estimated by averaging the divergence of
#' several reference points.
#'
#' The user should plot \eqn{S(t) Vs t} when looking for the maximal lyapunov
#' exponent and, if for some temporal range \eqn{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 \eqn{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.
#' @param time.series The original time series from which the maximal Lyapunov
#' exponent will be estimated.
#' @param min.embedding.dim Integer denoting the minimum dimension in which we
#' shall embed the time.series (see \link{buildTakens}).
#' @param max.embedding.dim Integer denoting the maximum dimension in which we
#' shall embed the time.series (see \link{buildTakens}).Thus,
#' we shall estimate the Lyapunov exponent between \emph{min.embedding.dim}
#' and \emph{max.embedding.dim}.
#' @param time.lag Integer denoting the number of time steps that will be use
#' to construct the Takens' vectors (see \link{buildTakens}).
#' @param radius Maximum distance in which will look for nearby trajectories.
#' @param theiler.window Integer denoting the Theiler window: Two Takens'
#' vectors must be separated by more than \emph{theiler.window} time steps in
#' order to be considered neighbours. By using a Theiler window, we exclude
#' temporally correlated vectors from our estimations.
#' @param min.neighs Minimum number of neighbours that a Takens' vector must
#' have to be considered a reference point.
#' @param min.ref.points Number of reference points that the routine will try
#' to use. The routine stops when it finds \emph{min.ref.points} reference
#' points, saving computation time.
#' @param max.time.steps Integer denoting the number of time steps marking the
#' end of the linear region.
#' @param number.boxes Number of boxes that will be used in the box assisted
#' algorithm (see \link{neighbourSearch}).
#' @param sampling.period Sampling period of the time series. When dealing
#' with a discrete system, the \emph{sampling.period} should be set to 1.
#' @param do.plot Logical value. If TRUE (default value), a plot of \eqn{S(t)}
#' Vs \eqn{t} is shown.
#' @param ... Additional plotting parameters.
#' @return A list with three components named \eqn{time} and \eqn{s.function}.
#' \eqn{time} is a vector containing the temporal interval where the system
#' evolves. It ranges from 0 to
#' \emph{\eqn{max.time.steps \cdot sampling.period}{max.time.steps * sampling.period}}.
#' \eqn{s.function} is a matrix containing the
#' values of the \eqn{S(t)} for each t in the time vector(the columns) and each
#' embedding dimension (the rows).
#' @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 \dontrun{
#' ## 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")
#' }
#' @author Constantino A. Garcia
#' @rdname maxLyapunov
#' @export maxLyapunov
#' @exportClass maxLyapunov
#' @useDynLib nonlinearTseries
maxLyapunov = function(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,
...) {
if (is.null(number.boxes)) {
number.boxes = estimateNumberBoxes(time.series, radius)
}
# s.function = matrix(0, ncol = (max.time.steps+1), nrow = n.embeddings)
s.function = .Call('_nonlinearTseries_lyapunov_exponent',
PACKAGE = 'nonlinearTseries',
time.series, min.embedding.dim, max.embedding.dim,
time.lag, radius, theiler.window, min.neighs,
min.ref.points, max.time.steps, number.boxes)
dimnames(s.function) = list(min.embedding.dim:max.embedding.dim,
0:max.time.steps)
time = (0:max.time.steps) * sampling.period
max.lyapunov.structure = list(
time = time,
s.function = s.function,
embedding.dims = min.embedding.dim:max.embedding.dim
)
class(max.lyapunov.structure) = "maxLyapunov"
id = deparse(substitute(time.series))
attr(max.lyapunov.structure,"id") = id
attr(max.lyapunov.structure,"time.lag") = time.lag
attr(max.lyapunov.structure,"theiler.window") = theiler.window
attr(max.lyapunov.structure,"min.neighs") = min.neighs
attr(max.lyapunov.structure,"min.ref.points") = min.ref.points
if (do.plot) {
tryCatch(
plot(max.lyapunov.structure,...),
error = function(e) {
warning("Error while trying to plot the Lyapunov object.")
})
}
max.lyapunov.structure
}
#' Returns the time in which the divergence of close trajectories was computed
#' in order to estimate the maximum Lyapunov exponent.
#' @param x A \emph{maxLyapunov} object.
#' @return A numeric vector representing the time in which the divergence of
#' close trajectories was computed.
#' @seealso \code{\link{maxLyapunov}}
#' @export divTime
divTime = function(x) {
UseMethod("divTime")
}
#' @return The \emph{divTime} function returns the
#' time in which the divergence of close trajectories was computed.
#' @rdname maxLyapunov
#' @export
divTime.maxLyapunov = function(x) {
x$time
}
#' @return The \emph{embeddingDims} function returns the
#' embeddings in which the divergence of close trajectories was computed
#' @rdname maxLyapunov
#' @export
embeddingDims.maxLyapunov = function(x) {
embeddingDims.default(x)
}
#' Returns the rate of divergence of close trajectories needed for the maximum
#' Lyapunov exponent estimation.
#' @param x A \emph{maxLyapunov} object.
#' @return A numeric matrix representing the time in which the divergence of
#' close trajectories was computed. Each row represents an embedding dimension
#' whereas that each column represents an specific moment of time.
#' @seealso \code{\link{maxLyapunov}}
#' @export divergence
divergence = function(x) {
UseMethod("divergence")
}
#' @return The \emph{divergence} function returns the
#' rate of divergence of close trajectories needed for the maximum Lyapunov
#' exponent estimation.
#' @rdname maxLyapunov
#' @export
divergence.maxLyapunov = function(x) {
x$s.function
}
#' @rdname maxLyapunov
#' @param main A title for the plot.
#' @param xlab A title for the x axis.
#' @param ylab A title for the y axis.
#' @param type Type of plot (see \code{\link[graphics]{plot}}).
#' @param col Vector of colors for each of the dimensions of the plot.
#' @param pch Vector of symbols for each of the dimensions of the plot.
#' @param add.legend add a legend to the plot?
#' @export
plot.maxLyapunov = function(x, main = "Estimating maximal Lyapunov exponent",
xlab = "time t", ylab = "S(t)", type = "p",
col = NULL, pch = NULL, add.legend = T, ...) {
embeddings = x$embedding.dims
n.embeddings = length(embeddings)
# obtain default parameter vectors if not specified
col = vectorizePar(col,n.embeddings)
pch = vectorizePar(pch,n.embeddings)
for (i in 1:n.embeddings) {
if (i != 1) {
lines(x$time[-1], x$s.function[as.character(embeddings[[i]]), -1],
type = type, col = col[[i]], pch = pch[[i]], ...)
} else {
# first time: use plot
plot(x$time[-1], x$s.function[as.character(embeddings[[i]]), -1],
xlab = xlab, ylab = ylab, main = main, type = type, col = col[[1]],
pch = pch[[1]], ...)
}
}
if (add.legend) {
legend("bottomright",
col = col[1:n.embeddings],
pch = pch[1:n.embeddings], lty = rep(1, n.embeddings),
lwd = rep(2.5, n.embeddings), bty = "n",
legend = embeddings, title = "Embedding dimension"
)
}
}
#' @return In order to obtain an estimation of the Lyapunov exponent the user
#' can use the \emph{estimate} function. The \eqn{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
#' \emph{use.embeddings} parameter. The slopes are determined by performing a
#' linear regression over the radius' range specified in \emph{regression.range}
#' @param ylim Numeric vector of length 2, giving the y coordinates range.
#' @param x A \emph{maxLyapunov} object.
#' @param regression.range Vector with 2 components denoting the range where
#' the function will perform linear regression.
#' @param use.embeddings A numeric vector specifying which embedding dimensions
#' should the \emph{estimate} function use to compute the Lyapunov exponent.
#' @param fit.col A vector of colors to plot the regression lines.
#' @param fit.lty The type of line to plot the regression lines.
#' @param fit.lwd The width of the line for the regression lines.
#' @rdname maxLyapunov
#' @export
estimate.maxLyapunov = function(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, ...) {
if (is.null(regression.range)) {
# the first position is always 0
min.time = x$time[[2]]
max.time = tail(x$time, 1)
} else {
min.time = regression.range[[1]]
max.time = regression.range[[2]]
}
if (!is.null(use.embeddings)) {
s.function = x$s.function[as.character(use.embeddings),]
} else {
use.embeddings = x$embedding.dims
s.function = x$s.function
}
indx = which(x$time >= min.time & x$time <= max.time )
x.values = x$time[indx]
lyapunov.estimate = c()
n.embeddings = length(use.embeddings)
# color and symbols for the plots if needed
if (do.plot) {
# obtain vector of graphical parameters if not specified
col = vectorizePar(col,n.embeddings)
pch = vectorizePar(pch,n.embeddings)
fit.col = vectorizePar(fit.col,n.embeddings,col)
}
for (i in 1:n.embeddings) {
current.embedding = use.embeddings[[i]]
y.values = s.function[as.character(current.embedding),indx]
fit = lm(y.values ~ x.values)
lyapunov.estimate = c(lyapunov.estimate,fit$coefficients[[2]])
if (do.plot) {
if (i != 1) {
lines(x$time, s.function[as.character(current.embedding), ],
type = type, col = col[[i]], pch = pch[[i]], ...)
} else {
if (is.null(ylim)) {
ylim = range(s.function)
}
plot(x$time, s.function[as.character(current.embedding), ],
xlab = xlab, ylab = ylab, main = main, type = type, col = col[[1]],
pch = pch[[1]], ylim = ylim, ...)
}
lines(x.values, fit$fitted.values, col = fit.col[[i]],
lty = fit.lty, lwd = fit.lwd)
}
}
lyapunov.estimate = mean(lyapunov.estimate)
if (do.plot && add.legend) {
legend("bottomright", col = col[1:n.embeddings],
pch = pch[1:n.embeddings], lty = rep(1, n.embeddings), bty = "n",
lwd = rep(2.5, n.embeddings), legend = use.embeddings,
title = "Embedding dimension")
}
lyapunov.estimate
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.