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=".")
}

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 closed 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 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].