Logistic map: Feigenbaum diagram in R
The other day I found some old basic code I had written about 15 years ago on a Mac Classic II to plot the Feigenbaum diagram for the logistic map. I remember, it took the little computer the whole night to produce the bifurcation chart.
logistic.map <- function(r, x, N, M){
## r: bifurcation parameter
## 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] <- r *z[i] * (1 - z[i])
}
## Return the last M iterations
z[c((N-M):N)]
}
## Set scanning range for bifurcation parameter r
my.r <- seq(2.5, 4, by=0.003)
system.time(Orbit <- sapply(my.r, logistic.map, x=0.1, N=1000, M=300))
## user system elapsed (on a 2.4GHz Core2Duo)
## 2.910 0.018 2.919
Orbit <- as.vector(Orbit)
r <- sort(rep(my.r, 301))
plot(Orbit ~ r, pch=".")
Let's not forget when Mitchell Feigenbaum started this work in 1975 he did this on his little calculator!
Update, 18 March 2012
The comment from Berend has helped to speedup the code by a factor of about four, thanks to byte compiling (using the same parameters as above), and Owe got me thinking about the alpha value of the plotting colour. Here is the updated result, with the R code below:
library(compiler) ## requires R >= 2.13.0
logistic.map <- cmpfun(logistic.map) # same function as above
my.r <- seq(2.5, 4, by=0.001)
N <- 2000; M <- 500; start.x <- 0.1
orbit <- sapply(my.r, logistic.map, x=start.x, N=N, M=M)
Orbit <- as.vector(orbit)
r <- sort(rep(my.r, (M+1)))
plot(Orbit ~ r, pch=".", col=rgb(0,0,0,0.05))
8 comments :
Way cool!
There is a small problem in your code. M is not defined outside the function logistic.map, so the call to sort does not work. Substituting 300 to M in the call everyting goes well.
Thanks for pointing this out. Fixed.
You can get a large speedup of this code if you use the compiler package.
Insert at top of file
library(compiler)
and after the definition of logistic.map insert
logistic.map <- cmpfun(logistic.map)
and run.
Thanks for sharing! So this is what byte compiling is about.
Here is my speedup measured on the same machine:
> library(compiler)
> logistic.map <- cmpfun(logistic.map)
> system.time(Orbit <- sapply(my.r, logistic.map, x=0.001, N=1000, M=300))
user system elapsed
0.641 0.011 0.648
Thanks, this is very nice. I tried to get the look of the plot in wikipedia (http://de.wikipedia.org/w/index.php?title=Datei:LogisticMap_BifurcationDiagram.png&filetimestamp=20050914001120) with the ggplot2 function, but I can't get the points small enough :(
I will try to change the alpha value of the colour. Maybe this helps.
Nice logistic maps. great to see your blog. Thanks....
http://www.suainlogistics.com
This code, with the same problem about M, was posted on Twitter/R Bloggers site today.
Post a Comment