Chaos and bifurcations with R

using lapply and ggplot2…

Alexandros Kouretsis

2019-03-09

Starting from Poincare’s research on instabilities of celestial trajectories in the end of 19th century, chaos theory has emerged as an important field of study for most of the natural sciences. The catalitic paper in the development of the field was written by Paul Lorentz more than 50 years ago. His research involved weather forcasting with statistical analysis of weather data. He experimented with various systems of differential equations and he made a simple yet crucial observation. That the evolution of a system is very sensitive to its initial conditions. He presented a system of three coupled differential equations modeling the fluid dynamics of weather within a simplistic framework. Plotting the trajectories of the system for particular initial conditions revealed not only the sensitivity of the system but also that there was an attractive region that within the system was bound to evolve ⊕In CRAN you can find excellent wrappers of low level functions like for example deSolve which is an interface to rock solid FORTRAN amd C solvers.

Yet was more to come, a few years later Hewlett Packard released the programmable pocket calculator HP-65, a cutting edge technology for its time. It was the first programmable handheld calculator in outer space, but also it was one of the first pocket programmable devices that “witnessed” chaos. Mitchell Feigenbaum started calculating the cosinus function iteratively and for any initial angle HP-65 after some iterations always returned ~0.739. Of course this discussion went much further and HP-65 could do much more complicated tasks that opened the discussion to the notion of universality in chaos, but lets remain at this crucial point and do the experiment in R, using the apply function. ⊕The sapply function starts an iteration along the values we provide in the first argument. A new environment is created in each step for the function we defined. If we don’t sent the calculated value to the parent environment it will be lost when the new environment will be invoked in the next step.

#get random angle
set.seed(1)
random_angle  <- sample(seq(1, 2*pi, length.out = 100), 1)
#start from here
get_value <- cos(random_angle)
n_steps <- seq(100)
#perfome iteration 150 times
bizarre_sequence <- 
  sapply(n_steps, function(x){
    get_value <<- cos(get_value)
  }) 
#get convergence/last value
tail(bizarre_sequence, 1)

## [1] 0.7390851

What we just did was to calculate the sequence \[x_{n+1} = cos x_{n}\] for 100 steps. Be cautious to the double arrow assignment to the parent environmet of the lapply fuction to update the value used by cosinus. We need the double arrow because the environments of the function calls done in each step of the lapply do not communicate but they are all children of the parent frame, which is whithin their scope. We can do this also with a classic and simple for loop but when you go the “lapply” way its fun and pays off when more complicated operations are involved.

The latter relation belongs to the family of one-dimensional maps and we can describe its behavior as a damped oscillation that converges to the stable value, ~0.739. This value is the solution to the transcendental equation \(x = con x\). A visualisation will make things morSe clear

After some iterations the system eventualy rests at ~0.739

ggplot() +
  geom_line(aes(x = n_steps, y = bizarre_sequence)) +
  ylab(expression('cos x'[n])) +
  xlab("n")

From the previous plot it is straightforward to see that once the iteration reaches ~0.739 it stabilises to this bizzare value. Now lets spice things up and add a free parameter. We have an angle so lets try to add a radius-like parameter, namely \[x_{n+1} = r cos x_n\]

Let’s do some calculations for a range of the \(r\) parameter \(1.5 \le r \le 2.9\) and see how the system behaves. In particular, we will calculate a sequence for nine distinct values of r.

set.seed(11)
getbizarre_sequences <- function(r, n_steps = 100){
  #get random angle
  random_angle  <- sample(seq(1, 2*pi, length.out = 100), 1)
  #get cos value
  get_value <- cos(random_angle)
  dt <- 
    sapply(seq(n_steps), function(x){
      get_value <<- r * cos(get_value)
    }) %>% data.table
  
  dt[ , n := seq(n_steps)]

}

r_range <- seq(from = 1.5, to = 2.9, length.out = 9)

#calculating sequencies
bizarre_list <- lapply(r_range, getbizarre_sequences)
names(bizarre_list) <- paste("r =", r_range)
#tidying data to long range
bizarre_sequences <- rbindlist(bizarre_list, idcol = TRUE)

We can inspect the data by exposing the first three values of the sequence for the first two values of r with the following chunk

dt <- bizarre_sequences[ , head(.SD, 3), by = .id]
kable(dt %>% head)

To inspect the behavior of the system for the nine different values of r we will plot the values of the sequences faceting on r.

ggplot(bizarre_sequences, aes(x = n, y = ., group = .id)) + 
  geom_line() +
  ylab(expression('r cos x'[n])) +
  xlab("n") + 
  facet_wrap( ~ .id, ncol = 3)

The above diagram depicts the transition of the system between chaotic and ordered phases as r changes. From the previous analysis we can safely conclude thet the system may rest to a stable fixed point (r = 1), may enter a cyclic behavior (r = 1.5) or a chaotic phase (r = 2.75). The critical values of r that the system changes state between stable, static and chaotic behavior due to changes in the parameters of the system is known as bifurcations.

To identify this state transition we can plot a Feigenbaun diagram. To get this diagram we let the system evolve for a considerable amount of iterations and collect the tail of the sequence generated. This process is repeated for a range of r.

r_range <- seq(from = 1, to = 3, length.out = 1000)

bizarre_list <- 
  lapply(r_range, function(x){
    tail(getbizarre_sequences(x, n_steps = 10000), 3000)
  })
names(bizarre_list) <- r_range
bizarre_sequences <- rbindlist(bizarre_list, idcol = TRUE)

Then, we plot r in the x-axis and \(r cos x\) in the y - axis.

bizarre_sequences[, .(r = as.numeric(.id), y = .)] %>%
  ggplot(aes(x = r, y = y)) +
  geom_point(size = 0.01, alpha = 0.01) +
  theme_light()

Voilà, this is one of the best visualisations to communicate the concept of bifurcations. Inspecting from the left to the right, the system rests after some iterations to a single value (a stable point), then enters a cyclic phase that repeats every two points (period-2 cycle), then enters a period-4 cycle etc. The system eventualy enters a chaotic region with acyclic behavior and for \(r \approx 2.3\) it enters again an ordered state. In conclusion, we can safely claim that as r grows the system enters and exits chaotic regions for finite windows of \(r\).

In the following posts I will continue the analysis of nonlinear dynamics and chaos using the R language. This particular topic is well rounded and will let us discuss topics from high performance computing, to data structures, multi-threading, visualisations and animation.