drawing a bifurcation diagram in Julia

bifurcation.jl

using Plots

# a vector of r values
Rs=collect(0.1:0.001:3)
T = 5000

N=zeros(length(Rs), T)

#Set t0 values to 1
N[:,1] .= 1

for (row, r) in enumerate(Rs), t in 2:T
    N[row, t] = N[row, t-1] + N[row, t-1] * r * ((K - N[row, t-1])/K)
end

w=100
all_Rs=repeat(Rs, inner = w)
all_Ns_array=[N[s, (T-(w-1)):T] for s in 1:size(N)[1]]

all_Ns=vcat(all_Ns_array...)

scatter(all_Rs, all_Ns,
 markercolor=:green,
  markerstrokecolor=:white,
  markersize=2,
  markerstrokewidth=0,legend=false,
  markeralpha = 0.1,
  xlabel = "Intrinsic rate of increase",
  ylabel = "Population size (100 final values)")
length(N)

Here is a slightly faster version:

rs = LinRange(0.1, 3.0, 2901)
K = 1e3
T = 5000
N = zeros(Float64, (T, length(rs)))

logmod(n, r, K) = n + n * r * (K-n)/K

for (j,r) in enumerate(rs)
    N[1,j] = 0.1
    for t in 2:T
        N[t,j] = logmod(N[t-1,j], r, K)
    end
end

F = N[end-99:end,:]
R = repeat(rs, inner=100)
points = unique(hcat(R, vec(F)); dims=1)

This is great and fun to play with; thanks!

Also, @aammd, K is not defined above.

An alternative with fewer functions

function orbit_log_map(xi,it) #set initial condition and number of iterations

f(x0,r) = r*x0*(1-x0) #logistic function
r = 2.8:0.001:4 #parameter interval
nr = length(r)
M = zeros(Float64, (it*nr+1,2)) #matrix

q = 2
M[q-1,1] = xi #initial condition
for ri in r   
    M[q-1,2] = ri
    p = q:1:(q+it-1)
    for t in p
        M[t,2] = ri
        M[t,1] = f(M[t-1,1], ri) #Solve function
    end
    q = q+it
end
return M
end

M = orbit_log_map(0.3,100) # reduce the number of iterations to speed the code

using Plots
scatter(M[:,2],M[:,1], ms=0.1, legend =:false, title="bifurcation diagram", 
    xlabel = "Paramenter (r)", ylabel = "x_{t+1}")