A Matlab program to iterate the logistic map - Dana Kester (with slight modifications by me)

bifur.m

% ICC 2016-04-06
% This is my attempt at some minor improvements to BifurDiag.m
% create a separate function
% I tried to use Matrices a bit more with R and X instead of itx and rv
% Dana Kester is plotting `rv` and `itx` and I'm plotting `R(:)` and `X(:)`
% The original BifurDiag builds up a very long vector which is what plotting expects
% But Andrew Ng taught me to put everything in matrices and only unroll them for plotting
% It means one line less for figuring out the offset of where to put the next set of iterations
% I think this function is also a bit more generic - with a step parameter
function [R, X] = bifur(rmin, rmax, iter, step)

% http://www.po.gso.uri.edu/tracking/tracking/chaos/presentations/bifurcationmicroscope/
%BifurDiag.m is a Matlab program to construct a bifurcation diagram for the logistic map
%to allow exploring the period doubling route to Chaos.  Dana Kester, Oct. 2002

%  The logistic map is x[t+1] = r*x[t]*(1-x[t]) where 0 < r < 4 is a parameter
%  value and 0 < x[t] < 1 is a variable.  The user specifies the following:
%  rmin = the lowest value of r to three decimals (or maybe two decimals)
%  rmax = the highest value of r to three decimals (or maybe two decimals)
%  x0 = initial value of x[t]
%  iter = number of iterations

%  WARNING THE CURRENT PARAMETERS rmin, rmax, and iter, (as well as the derived rmn, rmx)
%  MAY RESULT IN A LONG EXECUTION TIME -- ON MY 4-YEAR OLD COMPUTER WITH A PENTIUM II
%  166 MHZ PROCESSOR IT TOOK ABOUT 7 HOURS!  YOU CAN REDUCE THE NUMBER OF CALCULATIONS
%  BY INCREMENTING R-VALUES BY 0.01 RATHER THAN 0.001 BY CHANGING rmin and rmax TO TWO
%  DECIMAL PLACES AND CHANGING 1000 TO 100 IN LINES 24 & 25; OR YOU CAN NARROW THE 
%  RANGE OF R-VALUES; OR YOU MAY BE ABLE TO REDUCE THE NUMBER OF ITERATIONS AT EACH 
%  R-VALUE.  IT ALL DEPENDS ON WHAT YOU WANT TO EXAMINE IN THE BIFURCATION DIAGRAM.

% step = 0.01;
% iter = 400;
rmin += step;
rmax -= step;
x0 = 0.6;
iter2 = iter/2;

rmn = round(rmin/step); % Matlab requires integer subscripts
rmx = round(rmax/step);
% rct gets incremented at the end of the for loop
rct = 1; % rct is a counter for the number of r-values interated
R = zeros(iter, rmn-rmx);
X = zeros(iter, rmn-rmx);
for ridx = rmn:rmx
    x(1) = x0; % set initial condition--Matlab requires subscript > 0
    r = ridx * step; % converts back to decimal r
    for n = 2:iter
        x(n) = r * (x(n-1)) * (1 - (x(n-1)));
        % only store the values above a certain number of iterations
        if n > iter2
            R(n,rct) = r;
            X(n,rct) = x(n); % solution after the first iter2 iterations
        end
    end
    rct++;
end

end

BifurDiag.m

% http://www.po.gso.uri.edu/tracking/tracking/chaos/presentations/bifurcationmicroscope/
%BifurDiag.m is a Matlab program to construct a bifurcation diagram for the logistic map
%to allow exploring the period doubling route to Chaos.  Dana Kester, Oct. 2002

%  The logistic map is x[t+1] = r*x[t]*(1-x[t]) where 0 < r < 4 is a parameter
%  value and 0 < x[t] < 1 is a variable.  The user specifies the following:
%  rmin = the lowest value of r to three decimals (or maybe two decimals)
%  rmax = the highest value of r to three decimals (or maybe two decimals)
%  xo = initial value of x[t]
%  num = number of iterations

%  WARNING THE CURRENT PARAMETERS rmin, rmax, and num, (as well as the derived rmn, rmx)
%  MAY RESULT IN A LONG EXECUTION TIME -- ON MY 4-YEAR OLD COMPUTER WITH A PENTIUM II
%  166 MHZ PROCESSOR IT TOOK ABOUT 7 HOURS!  YOU CAN REDUCE THE NUMBER OF CALCULATIONS
%  BY INCREMENTING R-VALUES BY 0.01 RATHER THAN 0.001 BY CHANGING rmin and rmax TO TWO
%  DECIMAL PLACES AND CHANGING 1000 TO 100 IN LINES 24 & 25; OR YOU CAN NARROW THE 
%  RANGE OF R-VALUES; OR YOU MAY BE ABLE TO REDUCE THE NUMBER OF ITERATIONS AT EACH 
%  R-VALUE.  IT ALL DEPENDS ON WHAT YOU WANT TO EXAMINE IN THE BIFURCATION DIAGRAM.

clear
rmin = 2.001
rmax = 3.999
xo = 0.6
num = 400
num2 = num/2;

rmn = round(rmin*1000); %Matlab requires integer subscripts
rmx = round(rmax*1000);
rct = 0; %rct is a counter for the number of r-values interated
for r = rmn:rmx
    x(1)=xo; %set initial condition--Matlab requires subscript > 0
        rdec = r/1000; % converts back to decimal r
    for n=2:num
        x(n)=rdec*(x(n-1))*(1-(x(n-1)));
        if n > num2
            ir = n+(num-2)*rct; %ir is a counter for total interations
            itx(ir) = x(n); %after the first num2 iterations
            rv(ir)=rdec;
        end
    end
    rct=rct+1;
end


plot(rv,itx,'.')
axis([floor(rmin),ceil(rmax),0,1])
numitr = num2str(num);
title(['Bifurcation Diagram for the Logistic Map with Interations = ',numitr])

%  Note: The plot generated by this program needs some interactive editing in the Matlab GUI.
%  The default point size is "6", and I have not found a way to change that in the program.
%  If you enter PLOT EDIT mode you can select the data plotted and change the point size from
%  6 to 0.5 to show the detail of the results.

BifurDiag100.m

% I think the only point of the 100 in the file name 
% is that the default iterations are 100
clear
rmin = 2;
rmax = 4;
% default this to 100 then you should see results quicker
% the original default is 400
iter = 100;
step = 0.001;

[R, X] = bifur(rmin, rmax, iter, step);

% unroll the matrices to plot them
plot(R(:),X(:),'.','MarkerSize',2);
axis([rmin,rmax,0,1]);
title(sprintf('Bifurcation Diagram for the Logistic Map with iterations = %d, step = %.3f.', iter, step));