Creating a bifurcation diagram of the logistics map

George has already provided a solution, but let's try and build up to it piece by piece so we can understand what we are doing. Coming from the worlds of loops, Mathematica's programming style can be confusing.

First of all, we need a way of representing the function to be recursively applied. A pure function would be great for this. For instance:

r # (1 - #) &

We will then use NestList to apply this recursively a set number of times, in your case 1000 times, to the $0.01$ starting value:

NestList[r # (1 - #) &, 0.01, 1000]

This would generate a list of 1001 elements (including the starting point); we only want the last 100 elements of this list, so we use Part (i.e. [[...]]) and Span (i.e. ;;) to select those:

NestList[r # (1 - #) &, 0.01, 1000 - 1][[-100 ;;]]

We would like to keep track of the value of $r$ that generated each of those 100x values of the recursion, so we transform each element in the list into a sublist of the form {r, generatedvalue}, by mapping another simple pure function {r, #} & on the list obtained from NestList (see Map and its shorthand /@):

{r, #} & /@ NestList[r # (1 - #) &, 0.01, 1000][[-100 ;;]]

We now need to apply this recursion to each value of $r$ in your interval of interest. We can use Subdivide to generate such a list:

Subdivide[2.8, 4, 1200]

This will generate a list of $1201$ values, i.e. $1200$ equally spaced subdivisions as requested. We will use Table to run the recursion multiple times, each time with a different value of $r$ from that list above, and save the result in a variable called nestedlist:

nestedlist =
  Table[
    {r, #} & /@ NestList[r # (1 - #) &, 0.01, 1000][[-100 ;;]],
    {r, Subdivide[2.8, 4, 1200]}
  ];

Dimensions[nestedlist]
(*  Out: {1201, 100, 2} *)

The output of Dimensions shows us that our nestedlist contains 1201 lists, one per value of $r$, each of which contains 100 pairs of values, i.e. the {r, generatedvalue} pairs we wanted. This is too "nested"; we really only want a long list of such pairs, so we can "flatten out" the top level of this nested list, to obtain a list of $1201\times100=120100$ pairs. We use Flatten (together with a level specification, 1) for that, and we save the result in flatlist:

flatlist = Flatten[nestedlist, 1];
(* Out: {120100, 2} *)

Now we are ready to plot these points using ListPlot:

ListPlot[flatlist]

bifurcation


f[r_] := NestList[ r # (1 - #) &, .01, 1000][[-100 ;;]]
ListPlot[Flatten[Table[ {r, #} & /@ f[r] , {r, 2.8, 4, 1/1000}], 1]]

enter image description here

if you want to fix your loop approach it should look something like this:

L1 = Subdivide[2.8, 4, 1200];
ListPlot[
 Flatten[Table[
   x = Table[0, {1000}];
   x[[1]] = .01;
   Do[x[[i + 1]] = r*x[[i]] (1 - x[[i]]), {i, 999}];
   {r, #} & /@ x[[-100 ;;]], {r, L1}], 1] ]]

Tags:

Code Request