“Chaos is the science of surprises, of the nonlinear and the unpredictable. It teaches us to expect the unexpected.” - James Gleick
The world is full of natural visual patterns, from spots on a leopard to spirals of a fiddle-head fern. Some patterns are as small as the molecular arrangement of crystals and as big as the massive spiral pattern of the Milky Way Galaxy. Patterns can be found everywhere in nature. While one might think of patterns as uniform and regular, some patterns appear more random yet consistent. The definition of a pattern in nature is a consistent form, design, or expression that is not random. In the image below you see four strange yet mesmerizing patterns.
Some areas of these patterns are extremely bright, while others are rather dark. The bright areas seem to follow well-defined lines and curves. But these lines and curves don’t seem to follow any obvious logic like a circle or a spiral. The patterns have something “organic” about them and when seen from a distance they look like stardust in the night sky.
What makes these patterns even stranger is the fact that they resemble rather simple numerical systems of equations. These systems of equations are called strange attractors, which are some of the most fascinating structures in all of mathematics. The main aim of this article is to share fundamental and frontier ideas and techniques for modern science and technology that can stimulate research interest in the exploration of mathematical applications and directly pass the new knowledge to the next generation.
Starting from Poincare’s research on the instabilities of celestial trajectories at the end of the 19th century, chaos theory has emerged as an important field of study for most of the natural sciences.
The catalytic paper in the development of the field was written by Paul Lorentz more than 50 years ago. His research involved weather forecasting 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.
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 the chaos. Mitchell Feigenbaum started calculating the cosine 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.
Let’s remain at this crucial point and do the experiment in R, using the sapply( ) 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 send the calculated value to the parent environment it will be lost when the new environment will be invoked in the next step.
We want to calculate the sequence \(x_{n+1} =cos(x_n)\) for 100 steps.
#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.7390851We need the double arrow because the environments of the function call done in each step of the sapply do not communicate but they are all children of the parent frame, which is within their scope. We can do this also with a classic and simple for loop but when you go the “lapply” way it’s 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, of ~\(0.739\). This value is the solution to the transcendental equation \(x=cos(x)\) . A visualization will make things clear.
library(tidyverse)
ggplot() +
geom_line(aes(x = n_steps, y = bizarre_sequence)) +
ylab(expression('cos x'[n])) +
xlab("n")
In this plot, it is straightforward to see that once the iteration reaches ~\(0.739\) it stabilizes to this bizarre value. Now let us spice things up and add a free parameter. We have an angle so let us 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≤r≤2.9\) and see how the system behaves. In particular, we will calculate a sequence for nine distinct values of \(r\).
library(dplyr)
library(data.table)
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 sequences
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)
#View(bizarre_sequences)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
library(kableExtra)
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 that 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 are known as bifurcations.
As we have just seen changes in the parameters, produce some strange patterns for \(1.5≤r≤2.9\). Let’s study this phenomenon more rigorously. To do so we will study hundreds of different values for \(r\), all between \(1.5\) and \(2.9\), and examine which values of \(x\) the system eventually settles. To identify this state transition, 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)Now we plot \(r\) in the \(x\)-axis and \(rcos(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()
This resulting plot is called ‘Bifurcation Diagram’. First, starting from \(1.0\) there is a simple upward-sloping curve. After \(1.0\) somewhere and before it meets at \(1.5\) the curve splits in two. Just short of \(2.0\) these two curves split again. And somewhere after \(r=2.0\) the image starts to look quite bizarre. Let’s start at the left of the diagram, where the image is a simple curve. This curve means that our systems converge to a steady state. If \(r=1.0\) this steady state is somewhere in between \(0\) and \(1\). If we increase \(r\), this steady state also increases. This is why the curve increases. However, as we have seen before if \(r\) grows large enough, the system develops a cyclic behavior. Instead of reaching a unique steady state, the system switches back and forth between two particular values. This is why the initial curve in the bifurcation diagram splits into two somewhere after \(r>1.0\). The system now settles on two instead of a single point. These two points grow further apart if r increases. Once \(r\) is large enough, somewhere before \(2.0\), the system starts to settle on four different points. If \(r\) grows even further, the system becomes increasingly chaotic and constantly switches between many different points which are, however, all between \(-\ 3\) and \(3\). We are now in a chaotic region. The system is entirely deterministic, i.e. it’s not random, but it does not reach a fixed point solution or stable cyclic solution either, nor does the system explode. The system is chaotic.
So far we have only concerned ourselves with a single variable \(x\). We already found some very interesting types of chaotic behaviors. But things get a whole lot more interesting when we consider more than one variable. So let’s move from a single equation to a two-dimensional system of equations with two variables \(x\) and \(y\). In particular, we will focus on what are called two-dimensional quadratic maps, which in their most general form look like this:
\[ x_{t+1} = \alpha_1+\alpha_2x_t\ + \alpha_3x_t^2\ + \alpha_4x_ty_t \ + \alpha_5y_t+\alpha_6 y_t^2 \ \ \ \ \ \ \ \ \ \ ...(1) \]
\[ y_{t+1}= \alpha_{7}+ \alpha_8x_t \ + \alpha_9x_t^2+\alpha_{10}x_{t}y_{t} \ + \alpha_{11}y_t+\alpha_{12}y^2_t \ \ \ \ \ \ ...(2)\]
The two variables \(x\) and \(y\) depend on their previous period values in levels and squares as well as the combination of the previous period \(x\) and \(y\). Let’s write a function that takes a starting point \((x_0,y_0)\) and traces its position through each of \(n\) many iterations. Moreover, we will also write a function to plot these points. In this function, we will include an option to drop the first couple of points for which the attractor might not have yet reached its basin of attraction, which is the region where once reached, all future points will wander around in.
# function to run two-dimensional system
eqSys2d <- function(a=rep(0,12), n=1e5, p0=rep(1e-6,2)){
# matrix to store points
p <- matrix(NA,ncol=2,nrow=n)
# starting point
p[1,] <- p0
# iterations
for(i in 2:n){
x = p[i-1,1]
y = p[i-1,2]
p[i,1] <- a[1] + a[2]*x + a[3]*x^2 + a[4]*x*y + a[5]*y + a[6]*y^2
p[i,2] <- a[7] + a[8]*x + a[9]*x^2 + a[10]*x*y + a[11]*y + a[12]*y^2
}
# return
return(list("parameters"=a,"start"=p0,"points"=p))
}
# function to plot two-dimensional system
plotSys2d <- function(p, dropfirst=T,
pch=".", cex=1,
color="white"){
if(dropfirst) p <- p[((nrow(p) %/% 9) : nrow(p)), ]
par(mar=rep(1,4), bg="black", bty="n")
plot(p, cex=cex,
pch=pch, col=color,
xlab="", ylab="",
xaxt="n", yaxt="n")
}
# example
eqSys2d(a=1:12, n=10)Unsurprisingly, setting \(\alpha_1\) through \(\alpha_2\) equal to \(1\) through \(12\) causes the system to explode as indicated by the rapidly increasing values in the output that quickly surpass R’s range of feasible numerical values. So plotting this attractor would throw us an error (or at least a warning). The question now is which values for the \(12\) coefficients produce strange attractors.
One of the most prominent strange attractors based on a quadratic map is the so-called Hénon map named after French mathematician Michel Hénon. The Hénon map is typically written as
\[x_{t+1}=1-ax^2_t+y_t\]
\[y_{t+1}=bx_t\]
with \(a=0.9\) and \(b=0.3\). Another well-known chaotic quadratic map is the Tinkerbell map, which is commonly written as,
\[x_{t+1} = ax_t+by_t+x^2_t-y^2_t\]
\[ y_{t+1} = cx_t+dy_t+2x_ty_t \]
with \(a=0.9, b=−0.6013, c=2\ and\ d=0.5\). Now, all we need to do is to substitute these coefficients into equations (1),(2) and we are ready to plot the first of many hypnotic images. As we programmed different coloring options into our plot function, we can even produce plots with different colors.
# henon map
henon <- eqSys2d(a=c(1,0,-1.4,0,1,0,0,.3,0,0,0,0),n=1e4)
# tinkerbell map
tinkerbell <- eqSys2d(a=c(0,.9,1,0,-.6013,-1,0,2,0,2,.5,0),n=1e5)
# plots
par(mfrow=c(2,2))
plotSys2d(henon$points)
plotSys2d(tinkerbell$points)
# plots with color
plotSys2d(henon$points, color="maroon")
plotSys2d(tinkerbell$points, color="limegreen")
We are now ready to find our own strange attractors. The following code snippet searches and plots ten strange attractors. It does so by checking whether random combinations of our twelve coefficients constitute strange attractors.
n_attractors <- 10
counter <- 0
while(T){
# run new system of equations for 1000 steps
a <- runif(12,-1,1)
p <- eqSys2d(a=a, n=1e4)$points
# check if all points are unique
if(nrow(p) != nrow(unique(p))) next
# check if trajectory does not explode
if(any(is.infinite(p))) next
if(any(is.nan(p))) next
# check if lyapunov exponent is positive
d <- rep(NA,nrow(p)-1)
for(i in 2:nrow(p)) d[i-1] <- dist(p[(i-1):i,])
l <- sum(log2(d[-1]/d[-length(d)]))/nrow(p)
if(is.na(l) | l <= 0) next
# run again with more steps and plot
attractor <- eqSys2d(a=a, n=1e5)
plotSys2d(attractor$points)
# breaking condition
counter <- counter + 1
if(counter == n_attractors) break
}
rm(a,p,d,l)Every time you run this code you will get a whole new set of strange attractors. Some might look a little boring, others will produce the most fascinating patterns. Here are some examples:
The Hénon map and the Tinkerbell map are both situated in two-dimensional space. Many other famous strange attractors are, however, three-dimensional. Moreover, they are often formulated as differential equations:
\[\dfrac{dx}{dt}=f(x,y,z)\]
\[\dfrac{dy}{dt}=f(x,y,z)\]
\[\dfrac{dz}{dt}=f(x,y,z)\]
If we maintain our focus on quadratic maps, we only need to make a few small changes to our existing function to also simulate and visualize these three-dimensional attractors:
# function to run three-dimensional system
eqSys3d <- function(a=rep(0,30), n=1e5, p0=rep(1e-5,3), dt=1e-3){
p <- matrix(NA,ncol=3,nrow=n)
p[1,] <- p0
for(i in 2:n){
x = p[i-1,1]
y = p[i-1,2]
z = p[i-1,3]
dx <- a[1] +a[2]*x +a[3]*y +a[4]*z +a[5]*x^2 +a[6]*y^2 +a[7]*z^2 +a[8]*x*y +a[9]*x*z+a[10]*y*z
dy <- a[11]+a[12]*x+a[13]*y+a[14]*z+a[15]*x^2+a[16]*y^2+a[17]*z^2+a[18]*x*y+a[19]*x*z+a[20]*y*z
dz <- a[21]+a[22]*x+a[23]*y+a[24]*z+a[25]*x^2+a[26]*y^2+a[27]*z^2+a[28]*x*y+a[29]*x*z+a[30]*y*z
p[i,1] = p[i-1,1] + dt*dx
p[i,2] = p[i-1,2] + dt*dy
p[i,3] = p[i-1,3] + dt*dz
}
return(list("parameters"=a,"points"=p,start="p0"))
}
# function to plot three-dimensional systems
library(rgl)
plotSys3d <- function(p, dropfirst=T,
pch=".",
color=rainbow(nrow(p)),
theta=20, phi=10){
# drop initial observations
if(dropfirst) p <- p[((nrow(p) %/% 9) : nrow(p)), ]
# plot
library(plot3D)
par(bg = "black", bty="n")
scatter3D(p[,1],p[,2],p[,3],
pch=pch, bty="n",
colkey=F, axes=F,
theta=theta, phi=phi,
col=color
)
}With these tools at hand, we are now able to visualize some of the most famous attractors ever, including the butterfly-like Lorenz attractor or the Rössler attractor.
# lorenz attractor
lorenz <- eqSys3d(a=c(0, -10, 10, 0, 0, 0, 0, 0, 0, 0,
0, 28, -1, 0, 0, 0, 0, 0, -1, 0,
0, 0, 0, -8/3, 0, 0, 0, 1, 0, 0))
plotSys3d(lorenz$points)
# roessler attractor
roessler <- eqSys3d(a=c(0, 0, -1, -1, 0, 0, 0, 0, 0, 0,
0, 1, 0.2, 0, 0, 0, 0, 0, 0, 0,
0.2, 0, 0, -5.8, 0, 0, 0, 0, 1, 0))
plotSys3d(roessler$points, color="goldenrod2")
This is one of the best visualizations to communicate the concept of bifurcations. We will continue the analysis of nonlinear dynamics and chaos using the R language in the next article Chaos and Butterfly Effect.
Strogatz, Steven H. (1994). Nonlinear Dynamics and Chaos. Addison-Wesley. p. 262. ISBN 0-201-54344-3.
Luo, Dingjun (1997). Bifurcation Theory and Methods of Dynamical Systems. World Scientific. p. 26. ISBN 981-02-2094-4.
Afrajmovich, V. S.; Arnold, V. I.; et al. (1994). Bifurcation Theory and Catastrophe Theory. ISBN 978-3-540-65379-0.