4
$\begingroup$

My group partner and I are having trouble plotting a bifurcation plot of two recursive functions.

Our functions are as follows:

x[1 + n] == x[n] + z[n]
z[1 + n] == a z[n] - b Cos[x[n] + z[n]]

We want to plot x vs b but we aren't quite sure how to rearrange the formulas to allow for proper input for plotting, correctly use RSolve to generate functions which we can plot, or use parametric plotting.

We have used the Recursive table function to generate data but this was when our data was for a fixed b value. Other examples we have found use differential equations or nest functions, where our examples do not. When I input the formulas using RSolve in the format of another question of this site I get a value for x but just the input restated for y. What's the issue here?

h1 = RSolve[{x[p] == x[p - 1] + y[p - 1], x[0] == 1}, {x[p]}, p]
h2 = RSolve[{y[p] == a*y[p - 1] - b*Cos[x[p - 1] + y[p - 1]], y[0] == 1}, {y[p]}, p]

Again our ultimate goal is a bifurcation plot.

$\endgroup$
6
  • $\begingroup$ Could you double check the equations in your question? They seem to have multiple typos. $\endgroup$ – Chris K Nov 26 '16 at 20:25
  • $\begingroup$ What's an exemplary value of a? $\endgroup$ – corey979 Nov 26 '16 at 20:35
  • $\begingroup$ a=.5 in our case. The second sections of equations is just a shifted index, where x+1=p @corey979 $\endgroup$ – snowday1004 Nov 26 '16 at 20:36
  • $\begingroup$ And a range of b? $\endgroup$ – corey979 Nov 26 '16 at 20:44
  • $\begingroup$ we are trying to show once b is larger than 2 we start to get chaotic behaviour so from 0 to 6 @corey979 $\endgroup$ – snowday1004 Nov 26 '16 at 20:47
4
$\begingroup$

This answer only augments corey979's answer.

There may be a little more going on than you expect. I've amplified corey979's function (and made a possibly pointless stab at Compile[]ing the inner operation) to iterate over initial x, initial z, a, and b values, returning only the results of the last iteration in the form a big list of {b,x} pairs. (Aside: This really feels like a job for CUDA or OpenCL, which is on my "I should learn about those some time" queue.)

f[{xMin_, xMax_, dx_}, {zMin_, zMax_, dz_}, 
  {aMin_, aMax_, da_}, {bMin_, bMax_, db_}, iter_] :=
Module[{retVal, pairsStart, pairs, trans},
  retVal = {};
  pairsStart = 
    Flatten[
      N[Table[{x, z}, {x, xMin, xMax, dx}, {z, zMin, zMax, dz}]]
      ,1];
  Table[
    pairs = pairsStart;
    trans = Compile[{},
      Module[{i},
        For[i = 1, i <= iter, i++,
          pairs = Transpose[{
            pairs[[All, 1]] + pairs[[All, 2]], 
            a pairs[[All, 2]] - b Cos[pairs[[All, 1]] + pairs[[All, 2]]]
          }];
        ];
        retVal = retVal ~Join~ ({b, #} & /@ Flatten[pairs[[All, 1]]]);
        ]
      ];
    trans[];
    ,
    {a, aMin, aMax, da},
    {b, bMin, bMax, db}
  ];
  retVal
]

Example usage for initial x and z values of 0, 1, 2, 3, 4, 5, and 6, for a = 0.5, and for b = 3.5, taking 10 iterations:

f[{0, 2 π, 1}, {0, 2 π, 1}, {0.5, 0.5, 1}, {3.5, 3.5, 1}, 10]

(*{{3.5, -0.800026}, {3.5, -10.0644},   {3.5,  5.66061}, {3.5, 10.2941}, 
   {3.5, 11.8645},   {3.5,   5.63848},  {3.5,  5.83255}, {3.5, -0.899943}, 
   {3.5,  3.83456},  {3.5,   9.88333},  {3.5, 16.6377},  {3.5,  4.6175}, 
   {3.5,  5.66706},  {3.5,  -3.04366},  {3.5, 13.6011},  {3.5,  5.53968}, 
   {3.5,  5.58129},  {3.5,   5.64841},  {3.5,  5.69737}, {3.5,  5.36022}, 
   {3.5,  9.98684},  {3.5,   3.68096},  {3.5,  5.76446}, {3.5,  5.69186}, 
   {3.5,  3.81613},  {3.5,   4.00639},  {3.5,  9.9787},  {3.5, 30.5241}, 
   {3.5,  5.61322},  {3.5,   3.82061},  {3.5,  3.77938}, {3.5, 11.4152}, 
   {3.5, 11.0206},   {3.5,  18.0734},   {3.5, 66.4164},  {3.5,  3.75469}, 
   {3.5,  5.00735},  {3.5,  -0.824495}, {3.5, 18.2025},  {3.5, 18.217}, 
   {3.5, 51.748},    {3.5,  16.287},    {3.5,  4.95266}, {3.5,  5.70206}, 
   {3.5, 24.1745},   {3.5,  24.2642},   {3.5, 30.7148},  {3.5, 11.647}, 
   {3.5, 10.029}}  *)

Letting initial x and z range from 0 to in steps of 0.1, fixing a at 0.5, and letting b range from 0 to 6 in steps of 0.01, taking 100 iterates and keeping only those points with -π <= x <= 7π gives (after waiting a few seconds or a small number of minutes, depending on details of your computing environment)

ListPlot[
  Select[
    f[{0, 2 π, .1}, {0, 2 π, .1}, {0.5, 0.5, 1}, {0, 6, 0.01}, 
      100],
    -π <= #[[2]] <= 7 π &]
]

Mathematica graphics

It is perhaps unsurprising that your data is periodic in x with period . The fragments appearing for 3 <= x <= 3.2 don't appear to be a "small iterations" artifact. (Only takes a few times longer to run...)

ListPlot[
  Select[
    f[{2, 9, .05}, {0, 2 π, 0.1}, {0.5, 0.5, 1}, {2.95, 3.25, 
      0.001}, 500],
    0 <= #[[2]] <= 3 π &]
]

Mathematica graphics

and that feature for 3.02 <= b <=3.22 and 3<=x<=4 seems very familiar (if upside down) and seems to have its own fragments around b=3.165. This suggests you may have islands of chaos and some astoundingly twitchy metastable behaviour just past your first bifurcation. Zooming in a bit on the fragment just below (and using a slightly different f which reduces x modulo ,

ListPlot[
  Select[
    f[{3.1, 3.8, .01}, {0, 2 π, 0.05}, {0.5, 0.5, 1}, {3.05, 3.225, 
      0.0006}, 1000],
    3 <= #[[2]] <= 4 &]
]

Mathematica graphics

which is a very familiar object.

$\endgroup$
2
  • $\begingroup$ Really neat how you can see the chaos with the meta stable behaviour like you said! Thank you @Eric Towers $\endgroup$ – snowday1004 Nov 27 '16 at 19:14
  • 2
    $\begingroup$ @snowday1004 : BTW: I wouldn't "accept" my answer. There's nothing here that corey979 didn't demonstrate. $\endgroup$ – Eric Towers Nov 27 '16 at 20:58
4
$\begingroup$
f[a_, b_, max_, trans_] := Module[{x, z},
  x[0]  := 1;
  z[0]  := 1;
  x[n_] := x[n] = x[n - 1] + z[n - 1];
  z[n_] := z[n] = a z[n - 1] - b Cos[x[n - 1] + z[n - 1]];
  Drop[Table[{b, x[i]}, {i, 0, max}], trans]
  ]

f takes the values of a and b, and max+1 is the number of iterations of the map (including the initial conditions); trans is the number of transient iterations to be skipped. E.g.,

f[0.5, 0.6, 10, 0]

{{3.6, 1}, {3.6, 2}, {3.6, 3.99813}, {3.6, 7.3554}, {3.6, 7.31261}, {3.6, 5.43608}, {3.6, 2.11406}, {3.6, 2.31402}, {3.6, 4.85}, {3.6, 5.62416}, {3.6, 3.16512}}

The first number in each pair is the value of b (via ...{b, x[i]}... in the definition of f) to easily make a bifurcation diagram with

ListPlot[#, Frame -> True, PlotStyle -> {Black, PointSize[Tiny]}]& @
   Flatten[#, 1]& @ Table[f[0.5, i, 500, 400], {i, 0, 6, 0.01}]

enter image description here

You may want to play with the parameters, like PlotRange, PointSize, the steps, number of transients etc.

$\endgroup$
4
  • $\begingroup$ Nice! You might want to avoid including the transient dynamics by capturing only the last 50 points with Table[{b, x[i]}, {i, Max[0, max - 50], max}] in your definition of f. $\endgroup$ – Chris K Nov 26 '16 at 20:59
  • $\begingroup$ Thank you so much! @corey979 You just saved us a couple of hours banging our heads off the wall! A quick question: why would we have some of our variables turn black when we are entering them? They are originally blue when we enter them and they turn black, does that affect our code? $\endgroup$ – snowday1004 Nov 26 '16 at 21:33
  • $\begingroup$ thank you @ChrisK for your input as well! $\endgroup$ – snowday1004 Nov 26 '16 at 21:35
  • 1
    $\begingroup$ Blue is when the MMA doesn't know the variable, e.g. in a fresh session type a and it will be blue. When you assign a value, a=1, it turns black because MMA knows now what it is. All built-in functions/commands are black, obviously. In case of functions, defined with SetDelayed (:=), like f in the above answer, they will turn black after defining, but re-defining usually does not work - you have to Clear[f] beforehand. E.g., g[x_]:=x^2 and then g[x_]:=Sin[x] - g won't change, it will still remember it was an x^2. $\endgroup$ – corey979 Nov 26 '16 at 21:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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