Non-linear dynamics and chaos using R (1)

I have been teaching the course "Mathematics" for the PhD students at Ca' Foscari University for a few years. In the lectures I cover some "scattered" material that may prove useful, sooner or later, to develop quantitative models in Economics. One of the weeks of the course is devoted to non-linear dynamics, one-dimensional maps and chaos.
R can be extremely useful to explore the dynamics, compute fixed points or cycles and have a numerical look at the map. The mix of computational techniques and theoretical investigation is widely acknowledged to be fruitful. As Strogatz, "Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering", puts it:

The study of maps is still in its infancy, but exciting progress has been made in the last twenty year, thanks to the growing availability of calculators, then computers, and now computer graphics. Maps are easy and fast to simulate on digital computers where time is inherently discrete. Such computer experiments have revealed a number of unexpected and beautiful patterns…

Given a non-linear map \(f\) you may want to plot the time series of states \(x_{t+1}=f(x_t)\), depict the cobweb diagram and look at fixed points, to begin with. Let \(f\) be the immortal logistic map depending on the parameter \(r\).

f <- function(x,r) r*x*(1-x)
bounce <- function(f,init=4,n=10,cobweb=T,timeseries=F,dom=NULL,...){
iterates <- NULL
x0 <- init
for(t in 1:n){
  x1<- f(x0,...)
  iterates[t] <- x1
  x0 <- x1
  }

if(cobweb & !timeseries){
    if(is.null(dom)){
    a <- range(c(init,iterates))[1]-0.5
    b <- range(c(init,iterates))[2]+0.5} else
    {a <- dom[1];b <- dom[2]}
    curve(f(x,...),a,b);grid(col=1);abline(0,1,lty=2)
    lines(c(init,iterates), f(c(init,iterates),...),t="s")
    points(c(init,iterates), f(c(init,iterates),...))
}

if(timeseries){
    plot(0:n,c(init,iterates),t="l")
}

}
bounce(f,0.1,r=3.7,n=100)

bounce computes \(n\) iterates of the map starting from init an plots a cobweb diagram. Defaults may or may not work for a specific map but the plotting domain can be provided if the educated guess doesn't work (here, say, bounce(f,0.1,r=3.7,n=100,dom=c(0,1) would be probably better).

bounce(f,0.1,r=3.7,n=100)

The function bounce shows the 45-degrees line (dashed) and its intersections with the graph of the function: fixed points and their stability can be visually investigated (or found with uniroot). Indeed, the name comes from the repeated "bounces" on the line and on the graph. Setting timeseries=T displays the sequence of states of the dynamics (i.e., the time series) as shown below for \(r=3.8282\), a value that generates intermittency, see page 363 of Strogatz:

bounce(f,0.1,r=3.8282,n=100,timeseries=T)

In a coming post, we'll use R to draw bifurcation diagrams and Lyapunov exponents.