Mathematica code for Bifurcation Diagram

Question

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.

Answer 1

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},
  If[Head[cp] === Real,
    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]