In this blog post I will show you some ideas on how to plot a bifurcation diagram. I will use the well-known logistic map:

x[n+1] = r*x[n]*(1-x[n])


For every r we can start with some value x(0) and calculate x(100000) by repeatedly applying the above function. Where will it end up? And how does it depend on r?

For the rest of the blog I will use x(0) = 0.25 as the starting point.
Let’s look at the behavior for different r. For r = 0.4 we have exponential decay towards 0:

For r=1.6 it converges to a fixed value:

For r=2.7 it goes to a fixed value but it ‘overshoots’ each iteration:

For r=3.2 it oscillates between two values:

For r=3.5 it oscillates between 4 values:

For r=3.7 it looks like complete chaos, no nice repetition can be seen in the signal.

Bifurcation diagram

If we want to analyse this system, we would like to know what the system’s response is as a function of r. Of course I could make lots of plots, like I did above, for many many r’s, but this soon becomes cumbersome to interpret. A better way is to make a bifurcation diagram.

As seen in the plots above; the first few iteration really depend on the initial value x[0]. So to study the system we will focus on the values after many iteration, i.e. we will omit the first 1000 iterations.

The idea is now to plot as a function of r, the possible outcomes for x. We will plot this as a density map. For every r, we will computer x(1000) and the following 150000 values of x. The more likely a value of x occurs the darker we will make it on the y axis. If it doesn’t appear we will make it white, if it shows up many times it will be darker.

Details

In this plot we will vary r from 2.5 to 4 in 0.001 increments. We start with x[0]=0.25 and calculate the next 1000 to get rid of the transients. Next we will calculate 150000 values. These 150000 values are binned in bins from 0 to 1 with 0.001. (I.e. we tally how many values fall in the range [0,0.001] and [0.001,0.002] .. and [0.999,1]). To get the densities right we have to do some adjustments:
First we multiply the occurrences of each bin by the number of bins that are non-zero, then we clip all the occurrences by 2.5*numberofsteps (so 375000). Now we divide all the occurrences with the maximum values, such that all the values are between zero and one. Now all of these values can be interpreted as a grayscale.

Results

If you program all this (see the notebook below to check out the details!) you will see something like this:

*Click for a big image*

Conclusions

From this diagram we now have some idea what the system will do for a particular r. For r smaller than roughly 3.5 the reaction is simple, either 1, 2 or 4 values show up. If, however, we increase r, to (say) 3.7 the dynamics of the system is complicated and lots of values are possible. We can also see from the plot the points where the system `splits in two’ (a bifurcation), giving rise to period-doubling.

If you want to learn more about the logistic map, go to e.g. Wikipedia and MathWorld. The idea of this post is to show you on how to plot a bifurcation diagram in Mathematica. You can download the Mathematica code below, and apply it to your equations! The code should give you a good idea on how to program it.

Download the notebook:

This entry was posted on Wednesday, June 27th, 2012 at 3:08 pm and is filed under Math. You can follow any responses to this entry through the RSS 2.0 feed. Both comments and pings are currently closed.