The logistic map is a well-known and simple growth model that is defined by the iterative equation
\[x_{t+1} = 4rx_t(1-t_t)\],
where \(r\) is a parameter that can be thought of as a fertility and reproduction rate of the population. The allowed values of \(x\) range between 0 an 1 inclusively, where 0 means the population is extinct. The maximum of 1 can be interpreted as the ecological carrying capacity of the system. \(r\) is also allowed to vary between 0 and 1.
Let’s choose a value for \(r\) and see the development over time.
If \(r\) is not bigger than \(1/4\) than the population is deemed to extinct, as all numbers in the formula are lower than 1 which reduces their product in the next iteration.
Let’s check that in R. I will use the logistic map function from Markus Gesman’s blog.`
logistic_map <- function(r, x, N, M){
## r: growth rate
## x: initial value
## N: number of iteration
## M: number of iteration points to be returned
z <- 1:N
z[1] <- x
for(i in c(1:(N-1))){
z[i+1] <- 4 * r * z[i] * (1 - z[i])
}
## Return the last M iterations
z[c((N-M):N)]
}A package for plotting, text processing (glue) and data handling:
library(tidyverse)
library(glue)
library(gridExtra)Let’s define a work-horse function:
log_map_plot <- function(data = df, y, r, x, time = "time"){
out_plot <- ggplot(data) +
aes_string(x = time, y = y) +
geom_point(alpha = .5, size = .5) +
geom_line() +
coord_cartesian(ylim = c(0,1)) +
labs(title = paste0("r = ",r,", x = ",x),
y = "x",
x = "t")
print(out_plot)
return(out_plot)
}Now apply this function with some start values:
df <- data.frame(time = 1:100)
r = .20
x = .5
df$log_map1 <- logistic_map(r = r,
x = x,
N = 100,
M = 100)
p1 <- log_map_plot(data = df,
r = r,
x = x,
y = "log_map1")
p1
Doomed to extinguish. Runs depressingly to zero:
tail(df$log_map1)
#> [1] 9.23e-11 7.39e-11 5.91e-11 4.73e-11 3.78e-11 3.03e-11Now let’s try \(r=.33\).
r = .33
x = .5
df$log_map2 <- logistic_map(r = r,
x = x,
N = 100,
M = 100)
p2 <- log_map_plot(data = df,
r = r,
x = x,
y = "log_map2")
p2
Having crossed the border of \(.25\), we see that the behavior changed: Now the attractor is above zero, but still constant (and boring).
This behavior will continue until \(.75\), let’s check \(.7\) for instance:
r = .7
x = .5
df$log_map3 <- logistic_map(r = r,
x = x,
N = 100,
M = 100)
p3 <- log_map_plot(data = df,
r = r,
x = x,
y = "log_map3")
p3
Interestingly, some severe bouncing appeared but eventually, the curve flats out. Now let’s increase \(r\) to \(.8\):
r = .8
x = .5
df$log_map4 <- logistic_map(r = r,
x = x,
N = 100,
M = 100)
p4 <- log_map_plot(data = df,
r = r,
x = x,
y = "log_map4")
p4
Now the systems settles into a period-2 limit cycle oscillating between two values. A slight increase to \(r=.88\) will again change the system behavior:
r = .88
x = .5
df$log_map5 <- logistic_map(r = r,
x = x,
N = 100,
M = 100)
p5 <- log_map_plot(data = df,
r = r,
x = x,
y = "log_map5")
p5
Eventually, let’s see what happens if we set \(r=1\):
r = 1
x = .49999999
df$log_map6 <- logistic_map(r = r,
x = x,
N = 100,
M = 100)
p6 <- log_map_plot(data = df,
r = r,
x = x,
y = "log_map6")
p6
Now let’s compare the dependence on start values for \(r\approx 1\). Let’s choose a very similar start values, say \(x = .987654321\).
r = 1
x = .987654321
df$log_map7 <- logistic_map(r = r,
x = x,
N = 100,
M = 100)
p7 <- log_map_plot(data = df,
r = r,
x = x,
y = "log_map7")
p7
The pattern over time is completely different! Let’s try another again slightly different start value, say \(x = .987654320\).
r = 1
x = .987654320
df$log_map8 <- logistic_map(r = r,
x = x,
N = 100,
M = 100)
p8 <- log_map_plot(data = df,
r = r,
x = x,
y = "log_map8")
p8
Let’s compute the difference between log_map8 and log_map7.
df %>%
mutate(delta_map7_8 = abs(log_map8 - log_map7)) -> df
p9 <- log_map_plot(data = df,
r = .999999,
x = "Delta",
y = "delta_map7_8")
p9
Quite different.
Let’s plot two trajectories in one plot to better compare them.
p10 <- df %>%
select(time, log_map7, log_map8) %>%
gather(key = start_value, value = value, -time) %>%
ggplot() +
aes(x = time, y = value, color = start_value, shape = start_value) +
geom_point() +
geom_line() +
theme(legend.position = "none") +
labs(x = "t",
y = "x")
p10
At this point, let’s halt for a moment and put what we have done a little in context, or maybe just in one picture.
grid.arrange(p3, p4, p5,
p6, p7, p8,
nrow = 2)
As we see, slight differences in start values have dramatic effect on the development of the system behavior. Now, let’s analyse that a bit more systematically. Let’s compute the system development for a sequence of small differences for the start value (\(r\approx 1\) in each case).
start_values <- seq(from = 0, to = 1, by = .01)
start_values %>%
map_dfc(~logistic_map(r = .999999, x = ., N = 100, M = 100)) -> feigenbaum
feigenbaum <- feigenbaum[-1]What’s the correlation between the start values?
cor(feigenbaum$V2, feigenbaum$V3)
#> [1] -0.0915
cor(feigenbaum) -> cor_feigenbaum
cor_feigenbaum_upper.tri <- cor_feigenbaum[upper.tri(cor_feigenbaum)]library(corrplot)
corrplot(cor_feigenbaum)
What we can see that there’s basically not much correlation besides the main diagonal and, somewhat more interestingly, the “back main diagonal”. So, for instance, V2 and V99 seems to be nearly identical:
start_values
#> [1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12 0.13
#> [15] 0.14 0.15 0.16 0.17 0.18 0.19
#> [ reached getOption("max.print") -- omitted 81 entries ]
start_values[c(2,99)]
#> [1] 0.01 0.98Thus, there’s a symmetry around \(r=.5\). Let’s leave that aside and only consider the values between \(r=.01\) and \(r=.49\).
cor_feigenbaum2 <- cor_feigenbaum[1:49, 1:49]
cor_feigenbaum_upper.tri2 <- cor_feigenbaum2[upper.tri(cor_feigenbaum2)]
mean(cor_feigenbaum_upper.tri2)
#> [1] 0.00911ggplot(data.frame(cor_feigenbaum_upper.tri2)) +
geom_histogram(aes(x = cor_feigenbaum_upper.tri2))
So the correlation is essentially zero. Let’s repeat that: there’s no association between start values and the population development.
For instance, let’s take to neighboring start values and see what their population values are at \(t_{100}\):
feigenbaum %>%
filter(row_number() == 100) %>%
gather() %>%
ggplot(aes(x = key, y = value)) +
geom_line(group = 1) +
geom_point()
Let’s take two very close starting position (\(.00001\)), let the population run and see what happens:
start_values_close <- seq(from = 0, by = .00001, length.out = 100)
start_values %>%
map_dfc(~logistic_map(r = .999999, x = ., N = 100, M = 100)) -> feigenbaum2
feigenbaum2 <- feigenbaum2[-1]feigenbaum2 %>%
filter(row_number() == 100) %>%
gather() %>%
ggplot(aes(x = key, y = value)) +
geom_line(group = 1) +
geom_point()
Even with minimal differences, the system diverges dramatically.
If so simply systems can behave so unpredictable, how can we we hope to predict human behavior?
Finally, some fun plotting at the end. Let’s plot the Feigenbaum map.
r_values <- seq(from = .25, to = 1, by = .01)
r_values %>%
map_dfc(~logistic_map(r = ., x = .1, N = 1000, M = 100)) -> feigenbaum_r
names(feigenbaum_r) <- paste0("r_", r_values)
feigenbaum_r %>%
gather() %>%
separate(key, into = c("r_const", "r_value"), sep = "_") -> feigenbaum_r_long
feigenbaum_r_long %>%
mutate(r_value_num = as.numeric(r_value)) %>%
filter(between(r_value_num, 0.01, .999)) -> feigenbaum_r_longfeigenbaum_r_long %>%
ggplot(aes(x = r_value_num, y = value)) +
geom_point(size = .1)