Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
1  
related? mathematica.stackexchange.com/a/5131/1089 –  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
1  
You can add the x and y axes by adding Frame->True or Axes->True –  Gabriel Oct 18 '12 at 21:19
2  
You've seen Stan Wagon's book? –  J. M. Oct 18 '12 at 22:20
1  
show 1 more comment

1 Answer 1

up vote 6 down vote accepted

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]

enter image description here

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]

enter image description here

share|improve this answer
    
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. 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. Oct 19 '12 at 13:39
    
@J.M. Ahh, I'm thinking third edition. Again, I'll check in a bit. –  Mark McClure Oct 19 '12 at 13:46
    
@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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.