The strange attraction of the logistic map

Question

The purpose of the challenge is to approximately plot the attractor of the logistic map as a function of its parameter r (also called bifurcation diagram), or a subregion of it. The appearance of the graph can be seen in the following image from Wikipedia:

Background

The logistic map is a mathematical function that takes an input x_k and maps it to an output x_k+1 defined as

             x_k+1 = r x_k (1−x_k)

where r is the parameter of the map, assumed to lie in the interval [0, 4].

Given r in [0,4], and an initial value x₀ in the interval [0,1], it is interesting to repeatedly apply the function for a large number N of iterations, producing a final value x_N. Note that x_N will necessarily lie in [0,1] too.

As an example, consider r = 3.2, N = 1000. The initial value x₀ = 0.01 gives x₁₀₀₀ = 0.5130. For x₀ = 0.02 the result is x₀ = 0.7995. For any other initial values x₀ the final values x₁₀₀₀ are extremely close to either 0.5130 or 0.7995. This is seen in the graph as the height of the two lines at horizontal position r = 3.2.

This does not mean that for r = 3.2 each sequence converges to one of those two values. In fact, for the two initial values considered above, the sequences are (note the oscillating behaviour):

             x₀ = 0.01, ..., x₁₀₀₀ = 0.5130, x₁₀₀₁ = 0.7995, x₁₀₀₂ = 0.5130, ...
             x₀ = 0.02, ..., x₁₀₀₀ = 0.7995, x₁₀₀₁ = 0.5130, x₁₀₀₂ = 0.7995, ...

What is true is that for sufficiently large N, and for almost all initial values x₀, the term x_N will be close to one of the elements of the set {0.5130, 0.7995}. This set is called the attractor for this specific r.

For other values of the parameter r the size of the atractor set, or its elements, will change. The graph plots the elements in the attractor for each r.

The attractor for a specific r can be estimated by

1. testing a wide range of initial values x₀;
2. letting the system evolve for a large number N of iterations; and
3. taking note of the final values x_N that are obtained.

The challenge

Inputs

• N: number of iterations.

• r₁, r₂ and s. These define the set R of values of r, namely R = {r₁, r₁ + s, r₁ + 2 s, ..., r₂}.

Procedure

The set X of initial values x₀ is fixed: X = {0.01, 0.02, ..., 0,99}. Optionally, 0 and 1 may also be included in X.

For each r in R and each x₀ in X, iterate the logistic map N times to produce x_N. Record the obtained tuples (r, x_N).

Output

Plot each tuple (r, x_N) as a point in the plane with r as horizontal axis and x_N as vertical axis. Output should be graphic (not ASCII art).

Additional rules

• The indicated procedure defines the required result, but is not enforced. Any other procedure that procudes the same set of (r, x_N) tuples can be used.
• Input is flexible as usual.
• Floating point errors won't be held against the answerer.
• Graphic output is required, in any of the accepted formats. In particular, output may be displayed on screen, or a graphics file may be produced, or an array of RGB values may be output. If outputting a file or an array, please post an example of what it looks like when displayed.
• Graphics may be vector or raster. For raster graphics, the size of the image should be at least 400×400 pixels.
• Each point should be shown as a single pixel, or as a mark with size of the order of one pixel (otherwise the graph quickly gets cluttered).
• Axis range should be [0,4] for r (horizontal axis) and [0,1] for x_N (vertical axis); or it may be smaller as long as it includes all obtained points.
• Axis scales are arbitrary. In particular, the scale need not be the same for both axes.
• Grid lines, axis labels, colors and similar elements are acceptable, but not required.
• Shortest code in bytes wins.

Test cases

Click on each image for a high-resolution version.

N = 1000; r1 = 2.4; r2 = 4; s = 0.001;

N = 2000; r1 = 3.4; r2 = 3.8; s = 0.0002;

N = 10000; r1 = 3.56; r2 = 3.59; s = 0.00002;

Acknowledgment

Thanks to @FryAmTheEggman and @AndrasDeak for their helpful comments while the challenge was in the sandbox.

Answer 1

MATL, 32 30 28 27 bytes

4 bytes saved thanks to @Luis

3$:0:.01:1!i:"tU-y*]'.'3$XG

The input format is r1, s, r2, and N

Try it at MATL Online

Explanation

        % Implicitly grab the first three inputs
3$:     % Take these three inputs and create the array [r1, r1+s, ...]
0:.01:1 % [0, 0.01, 0.02, ... 1]
!       % Transpose this array
i       % Implicitly grab the input, N
:"      % For each iteration
  tU    % Duplicate and square the X matrix
  -     % Subtract from the X matrix (X - X^2) == x * (1 - x)
  y     % Make a copy of R array
  *     % Multiply the R array by the (X - X^2) matrix to yield the new X matrix
]       % End of for loop
'.'    % Push the string literal '.' to the stack (specifies that we want
        % dots as markers)
3$XG    % Call the 3-input version of PLOT to create the dot plot

Answer 2

Mathematica, 65 bytes

Graphics@Table[Point@{r,Nest[r#(1-#)&,x,#]},{x,0,1,.01},{r,##2}]&

Pure function taking the arguments N, r1, r2, s in that order. Nest[r#(1-#)&,x,N] iterates the logistic function r#(1-#)& a total of N times starting at x; here the first argument to the function (#) is the N in question; Point@{r,...} produces a Point that Graphics will be happy to plot. Table[...,{x,0,1,.01},{r,##2}] creates a whole bunch of these points, with the x value running from 0 to 1 in increments of .01; the ##2 in {r,##2} denotes all of the original function arguments starting from the second one, and so {r,##2} expands to {r,r1,r2,s} which correctly sets the range and increment for r.

Sample output, on the second test case: the input

Graphics@Table[Point@{r,Nest[r#(1-#)&,x,#]},{x,0,1,.01},{r,##2}]&[2000,3.4,3.8,0.0002]

yields the graphics below.

Answer 3

Mathematica, 65 bytes

I used some of Greg Martin's tricks and this is my version without using Graphics

ListPlot@Table[{r,NestList[#(1-#)r&,.5,#][[-i]]},{i,99},{r,##2}]&

input

[1000, 2.4, 4, 0.001]

output

input

[2000, 3.4, 3.8, 0.0002]

output

Answer 4

TI-Basic, 85 bytes

Prompt P,Q,S,N
P→Xmin:Q→Xmax
0→Ymin:1→Ymax
For(W,.01,1,.01
For(R,P,Q,S
W→X
For(U,1,N
R*X*(1-X→X
End
Pt-On(R,X
End
End

A complete TI-Basic program which takes input in the order r1,r2,s,N and then shows the output in real time on the graph screen. Note that this tends to be incredibly slow.

Here is an incomplete sample output generated after about 2.5 hours for the input 3,4,0.01,100:

Answer 5

ProcessingJS, 125 123 120 bytes

Thanks to Kritixi Lithos for saving 3 bytes.

var f(n,q,r,s){size(4e3,1e3);for(i=0;i<1;i+=.01)for(p=q;p<=r;p+=s){x=i;for(j=0;j<n;j++)x*=p-p*x;point(p*1e3,1e3-x*1e3)}}

Try it online! Call using f(N, r_1, r_2, s);

Answer 6

GEL, 158 bytes

`(N,r,t,s)=(LinePlotWindow=[r,t,0,1];for i=r to t by s do(p=.;for w=0to 1by 0.01do(x=w;for a=0to N do(x=i*x*(1-x););p=[p;q=[i,x]];);LinePlotDrawPoints(p);););

It may not be the shortest, but it draws in real time, although it can be incredibly slow with huge inputs. Anyways, this is an anonymous function which takes input in the format (N,r1,r2,s) and outputs the plot in a new window. Note that this must be run with the GNOME version of Genius.

Answer 7

R, 159 147 bytes

pryr::f({plot(NA,xlim=c(a,b),ylim=0:1);q=function(r,n,x=1:99/100){for(i in 1:n)x=r*x*(1-x);x};for(i in seq(a,b,s))points(rep(i,99),q(i,n),cex=.1)})

Which productes the function

function (a, b, n, s) 
{
    plot(NA, xlim = c(a, b), ylim = 0:1)
    q = function(r, n, x = 1:99/100) {
        for (i in 1:n) x = r * x * (1 - x)
        x
    }
    for (i in seq(a, b, s)) points(rep(i, 99), q(i, n), cex = 0.1)
}

plot(NA,...) creates an empty canvas that has the correct dimensions. q is the function that does the iterating. It takes a value of r, and then does n iterations for all starting points between 0.01 and 0.99. It then returns the resulting vector.

The for-loop applies the function q to the sequence a to b with step s. Instead of returning the values, it adds them as points to the plot. If the attraction point is one value, all the points would just overlap and show as one point. cex=.1 is a necessary addition to make the points as small as possible.