# Mathematica code for Bifurcation Diagram

At the moment I am trying to construct a bifurcation diagram of the iterative function $f(x)=$ $ax-1.1975x^3$. I've scoured the internet for pre-made bifurcation diagrams and found many (mostly of the logistic map). However, as the code is quite complicated I am not sure how to edit the code so that it deals with my function instead of the logistic one. Would anyone have a general template for the code to create a bifurcation diagram of a function? Ideally, I would like to have $a$ on the x-axis and equilibrium values on the y-axis.

• Commented Oct 18, 2012 at 20:50
• @chris Thank you for your suggestion. Unfortunately, I tried using that code and it's not quite what I'm looking for. I require x and y axes, but that code you suggested only provides a graphic.
– Mel
Commented Oct 18, 2012 at 21:01
• You can add the x and y axes by adding Frame->True or Axes->True Commented Oct 18, 2012 at 21:19
• You've seen Stan Wagon's book? Commented Oct 18, 2012 at 22:20
• mathematica.stackexchange.com/questions/5123 Commented Oct 19, 2012 at 8:41

There are two aspects of this question that distinguish it from previous questions:

1. The request for a general template, as opposed to just a single example.
2. The fact that the given example is a polynomial of degree three, whereas as opposed to the quadratic examples which appear in many places, including Stan's book.

To deal with the first issue, in part, let's simply define the function and then refer only to that definition throughout the code.

f[a_][x_] := a*x - 1.19 x^3;


A well known and important fact in dynamics is that each attractive orbit must attract at least one critical point. Thus, to detect attractive behavior for a given $a$, we should iterate from each critical point. Let's find the critical points in terms of $a$.

cps[a_] = x /. Quiet[Solve[D[f[a][x], x] == 0, x],
Solve::ratnz]


Next, given a parameter value $a$ and critical point cp, we need to find the resulting critical orbit, after dropping some transient behavior. Then we write a function points which does this for each critical point.

criticalOrbits[a_, cp_] := Module[{try},
try = NestWhileList[f[a], cp, Abs[#] < 100 &, 1, 500];
If[Abs[Last[try]] >= 100, try = {},
try = Drop[{a, #} & /@ try, 100]
],{}]];
points[k_] := Partition[Flatten[Table[criticalOrbits[a, cps[a][[k]]],
{a, -2, 4, 0.002}]], 2]


A little experimentation shows that a natural range for the parameter $a$ would be $0$ to $3$. I've allowed $a$ to range from $-2$ to $4$ to illustrate the fact that the code takes care to exit gracefully if given a divergent orbit or non-real critical point is input - necessary, if we would like this to work with a variety of functions.

Finally, we generate the image using color to differentiate the orbits of the critical points.

Graphics[{Opacity[0.02], PointSize[0.002],
Table[{ColorData[1, k], Point[points[k]]},
{k, 1, Length[cps[a]]}]}, Frame -> True]


• From page 207 of Wagon's book: "Note that BifurcationPlot was written so that we can use any function in place of $x(1-x)$..." Commented Oct 19, 2012 at 13:15
• @J.M. I don't have the text handy but I actually collaborated with Stan on that chapter and related material, although I'm only credited in the complex dynamics chapter. I don't recall the text dealing with multiple critical points, which it must do to work as advertised. I certainly could be mistaken, though, and will look at a copy a bit later today. Commented Oct 19, 2012 at 13:32
• Yes, I saw your name in the bibliography. :) If memory serves, the version in the second edition indeed is unable to cope with things more complicated than the good ol' logistic map, but the new edition seems to imply that things have changed... Commented Oct 19, 2012 at 13:39
• @J.M. Having looked into this further, I stand by my assertion that Stan's code is not ready to accept a cubic with parameter. I quite literally sat down with him, opened up the notebook file that forms that chapter in his book, and tried it. Boom. Commented Oct 19, 2012 at 22:41
• @MarkMcClure Yes, I was using the cloud to test the code already. Now, using your notebook, I can see the output. I guess, there was a problem with my copy-pasting. Thanks for the responses!
– Our
Commented Jan 2, 2021 at 12:42