Chaos, bifurcation diagrams and Lyapunov exponents with R (2)

(The first part of this article can be read here)

Iteration of one-dimensional maps can generate stunning complexity and famed examples of chaotic behavior.  R can be used to get the flavor of this richness and reproduce some of the most famous pictures in the history of science, such as the bifurcation diagram of the logistic map or the representation of its Lyapunov exponents.

Given a one dimensional map depending on a parameter, a bifurcation diagram shows the stable structures (fixed point, cycles, attractors) visited by the dynamics for each value of the so called "bifurcation parameter".  Resisting the temptation to use again and again the beloved logistic map, consider instead the following dynamical system

x(t+1)=a cos x(t),  

taken from exercise 10.2.6, page 389 of Strogatz, "Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering".  The following script defines a function to get the bifurcation diagram for a given map (I modified some existing code, but I was not able to find and credit the exact source: https://web.stanford.edu/group/heeh/cgi-bin/web/node/59, say, has the same ideas. I defined a function and removed the use of "expression" and "evalf").  See also the post at https://www.r-bloggers.com/dynamical-systems-mapping-chaos-with-r/

bif_diagram <- function(f=function(x,a) (a*x*(1-x)),alow=2.5,ahigh=4,
                        thinness=1000, transient=200, collect=200){
        # f function, parameter must be named a
    n <- 1
    R <- seq(alow,ahigh,length=thinness)
    data <- matrix(0,collect,thinness+1)

    for(a in R){
      x <- runif(1) # random initial condition
      ## first converge to attractor
      for(i in 1:transient){
        x <- f(x,a)
      } # collect points on attractor
      for(i in 1:collect){
        x <- f(x,a)
        data[i,n] <- x
      }
  n <- n+1
}

data <- data[,1:thinness]
yrange <- range(data)+c(-0.1,0.1)
plot(R,data[1,], pch=".", xlab="a", ylab="States",ylim=yrange)
for(i in 2:collect) points(R,data[i,],pch=".")
}

Now, enter

f <- function(x,a) a* cos(x)
bif_diagram(f,alow=0.5,ahigh=4)

Bifurcation diagram for f(x,a)=a cos x, when a is the range [0.5,4].

You can see that, for low values of the parameter a, there are unique fixed points or simple cycles.  Then, through  a series of (quite typical) period-doubling bifurcations, chaos appears and suddenly disappears when the parameter crosses 3.  bif_diagram has defaults parameters that can be adjusted to specify the parameter's interval of interest or increase/decrease the resolution and legibility of the resulting plot.

One of the most interesting signature of chaos is the divergence of orbits arbitrarily close in their initial conditions.  Broadly speaking, orbits are stretched apart by some systems and a positive stretching rate is signaling the presence of chaos.  The Lyapunov Exponents (LE) is the average (exponential) growth rate of the divergence of initially nearby orbits.  LEs can be computed, for any value of the bifurcation parameter using the lyap function below, where is assumed that you have properly defined the dynamics f in advance.

lyap <- function(a,trans=300,num=1000){
    x0 <- runif(1)
    for(time in 1:trans){
        x1 <- f(x0,a=a);x0 <- x1
    }
    sl <- 0
    for(time in 1:num){
        x1 <- f(x0,a=a);x0 <- x1
        sl <- sl+log(abs(grad(f,x1,a=a)))
    }
    sl/num
}

In order to graph LEs define a sequence of values of the parameters, use sapply to get the sequence of corresponding exponents by applying lyap to the elements of the sequence and finally plot.  With the same map given above:

a <- seq(0.5,4,length=100)
ly <- sapply(a,lyap)
plot(a,ly,t="l")

Lyapunov exponents for f(x,a)=a cos x, when a is the range [0.5,4].

It can be seen, say, that when a=2, the LE is positive and chaos is in action [Check the bifurcation diagram to get the same intuition for that value of a].  Entering lyap(2) would indeed produce 0.15 (approximately).  In rough terms, this means that divergence of close initial points will be on average amplified by about 15% per iteration of the map over its domain.

To conclude:

• it is fun to try the same things with a slightly modified map: x(t+1)=cos(a x(t)) (the parameter is "inside" the cos);
• if you want to replicate two of the most well-known pictures related to the logistic map, type

  f <- function(x,a) a*x*(1-x)
  bif_diagram(f,0,4)
  

• the graph of LE, seen on page 369 of Strogatz's book is readily obtained through

a <- seq(3,4,len=100)
f <- function(x,a) a*x*(1-x)
ly <- sapply(a,lyap,num=2000)
plot(a,ly,t="l",ylim=c(-1,1));abline(h=0)

Lyapunov exponens for the logistic map, when a is the range [3,4].