logisticmap.jl

using Base.Test
using Gadfly

logistic_growth = (x, r) -> r*x*(1-x)

function logistic_map(x::Number, r::Float64, n::Int64)
   @test x > 0
   @test x <= 1
   orbit = zeros(Float64, n)
   orbit[1] = x
   for i in 2:n
      orbit[i] = logistic_growth(orbit[i-1], r)
   end
   return orbit
end

## Single diagram

draw(PDF("map.pdf", 6inch, 6inch),
   plot(x = 1:50, y = logistic_map(0.3, 3.5, 30))
)

## Bifurcation diagram

r_space = linspace(1.0, 4.0, 20000)
n = 2000
x0 = 0.4

all_r = Float64[]
all_x = Float64[]

@time for r in r_space
   orbit = logistic_map(x0, r, n)
   orbit = logistic_map(orbit[end], r, n)
   uniques = unique(orbit[1500:end])
   for u in uniques
      append!(all_r, [r])
      append!(all_x, [u])
   end
end

draw(PNG("bifurcation.png", 100cm, 70cm),
   plot(x = all_r, y = all_x,
      Theme(default_point_size=.1mm, default_color=color("black"), highlight_width=0mm, panel_fill=color("white"))
   )
)