# 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.

• – chris Oct 18 '12 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 Oct 18 '12 at 21:01
• You can add the x and y axes by adding Frame->True or Axes->True – Gabriel Oct 18 '12 at 21:19
• You've seen Stan Wagon's book? – J. M.'s ennui Oct 18 '12 at 22:20
• mathematica.stackexchange.com/questions/5123 – Vitaliy Kaurov Oct 19 '12 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)$..." – J. M.'s ennui Oct 19 '12 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. – Mark McClure Oct 19 '12 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... – J. M.'s ennui Oct 19 '12 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. – Mark McClure Oct 19 '12 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 Jan 2 at 12:42