**A blog from Sydney**, and kindly contributed to R-bloggers)

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

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

`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.

*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.

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

**leave a comment**for the author, please follow the link and comment on their blog:

**A blog from Sydney**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...