Skip to content

Instantly share code, notes, and snippets.

@aniline
Created February 21, 2020 14:39
Show Gist options
  • Save aniline/a2899bb98c792fe1de5b512d4b6821a8 to your computer and use it in GitHub Desktop.
Save aniline/a2899bb98c792fe1de5b512d4b6821a8 to your computer and use it in GitHub Desktop.
Render bifurcation diagram (octave script)
## Period doubling bifurcation
## https://en.wikipedia.org/wiki/Logistic_map
1;
function ret = popfn(prev, r)
ret = r * (prev) * (1-prev);
endfunction
# Iterate till n-iterations
function ret = popfniter(numiter, r, initial)
nv = ones(1, numiter+1) .* initial;
for iter = 2:numiter+1
nv(iter) = popfn(nv(iter-1), r);
if nv(iter) > 10000,
nv(iter) = 10000;
elseif nv(iter) < -10000,
nv(iter) = -10000;
endif;
endfor;
ret = arrayfun(@(x) round(x*10000)/10000, nv);
endfunction
## List of converging values
function ret = popfconv(numiter, r, initial)
ll = popfniter(numiter, r, initial);
[u, ~, j] = unique(ll);
aa = [accumarray(j', 1), u'];
aa = sortrows(aa, 1);
numelems = round((numiter+1)/aa(end,1))-1;
ret = aa(end-numelems:end);
endfunction;
## list of points.
function ret = poplot(numiter, rrange, initial)
points = [];
for r = rrange
m = popfconv(numiter, r, initial);
points = [points; [ones(1, length(m)) * r; m]'];
endfor
ret = points;
endfunction;
m = poplot(2000, 1:0.02:4, 0.5);
scatter(m(:,1)', m(:,2)', '.')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment