# How to plot a bifurcation diagram of a system of recursive formulas in terms of two variables

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.

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

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_] :=
retVal = {};
Flatten[
N[Table[{x, z}, {x, xMin, xMax, dx}, {z, zMin, zMax, dz}]]
,1];
Table[
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 2π 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 π &]
]


It is perhaps unsurprising that your data is periodic in x with period 2π. 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 π &]
]


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 2π,

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 &]
]


which is a very familiar object.

• Really neat how you can see the chaos with the meta stable behaviour like you said! Thank you @Eric Towers Nov 27, 2016 at 19:14
• @snowday1004 : BTW: I wouldn't "accept" my answer. There's nothing here that corey979 didn't demonstrate. Nov 27, 2016 at 20:58
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}]


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

• 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. Nov 26, 2016 at 20:59
• 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? Nov 26, 2016 at 21:33
• thank you @ChrisK for your input as well! Nov 26, 2016 at 21:35
• 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. Nov 26, 2016 at 21:39