Why isn't this plot looking correct?

1 view (last 30 days)
Tyler
Tyler on 8 Oct 2014
Answered: Mischa Kim on 8 Oct 2014
Calculate P250 for 1000 different r values, ranging between 3 ≤ r ≤ 4 (linspace would be a good tool to use here), and save them as you go in a big vector. Plot your r values against your P250 values. What do you see? How are the points changing as r increases?
% Plotting a dynamical system
r = linspace(3,4,1000); %creates vector of length 1000 with evenly spaced values between 3 and 4
seed = linspace(0.1,0.9,100);
P_250 = zeros(1,1000); %creates zero vector of length 1000
hold on
for s = 1:100
P = zeros(1,250); %creates zero vector of length 25
seed_index = seed(s);
P(1) = seed_index; %seed value
for r_index = 1:1000
r_current = r(r_index);
for n = 1:250
P(n + 1) = r_current*P(n)*(1 - P(n)); %dynamical system
end
P_250(r_index) = P(250); %stores value of P(250) for each value of
end
end
plot(r,P_250,'.');
The plot just doesn't look right when this runs. I can't find the problem.

Answers (1)

Mischa Kim
Mischa Kim on 8 Oct 2014
Tyler, I am not sure what you are looking for. However, it seems like you did get what the exercise is asking for. Note, that the exercise is not asking you to create the entire bifurcation diagram. If you increase the number of data points (e.g. 100000) and add another bifurcation branch the plot starts looking more like the diagram you would expect for the logistic map.
r = linspace(3,4,100000); %creates vector of length 1000 with evenly spaced values between 3 and 4
seed = linspace(0.1,0.9,100);
P_250 = zeros(1,100000); %creates zero vector of length 1000
P_249 = zeros(1,100000); %creates zero vector of length 1000
hold on
for s = 1:100
P = zeros(1,250); %creates zero vector of length 25
seed_index = seed(s);
P(1) = seed_index; %seed value
for r_index = 1:100000
r_current = r(r_index);
for n = 1:250
P(n + 1) = r_current*P(n)*(1 - P(n)); %dynamical system
end
P_249(r_index) = P(249); %stores value of P(250) for each value of
P_250(r_index) = P(250); %stores value of P(250) for each value of
end
end
plot(r,[P_250;P_249],'b.');

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!