# Trouble with plotting a bifurcation diagram

I have the function $f(x)=x^3 + a x$ and I have to plot a bifurcation diagram for values of $a$ between 0 and 2. I have used the code below but it doesn't seem to work for me:

ListPlot[ParallelTable[
Thread[{a, Nest[x^3 + ax &, Range[0, 1, 0.01], 1000]}], {a, 0, 2}],
PlotStyle -> PointSize[0]]


Was wondering if anything is wrong with the code?

• 1) ax is not the same as a*x. 2) You need to use # instead of x in Nest. 3) Is this for the differential equation x'[t]==f[x[t]] or the difference equation x[t+1]==f[x[t]]? – Chris K Nov 30 '16 at 0:37
• You've got terrible overflows; see what's happening with NestList[#^3 + 0.1 # &, 0.5, 500]. Additionally, you should reference where you took the piece of code you're trying to adapt. – corey979 Nov 30 '16 at 0:43
• Many thanks for your help! :) – NDZS Nov 30 '16 at 22:16

The longterm behaviour of this map is either going to zero or blowing up to infinity. Hence the bifurcation diagram is just a set of points that hit "walls" progressively leftward as a increases.

To gain insight:

f[a_, x_] := x^3 + a  x
roots[a_] := {x, x} /. Solve[f[a, x] == x && x >= 0, x]
nf[a_, x_, n_] :=
Catenate[{{##}, {#2, #2}} & @@@
NestList[{#[[2]], f[a, #[[2]]]} &, {x, f[a, x]}, n]]
Manipulate[
Module[{arrows = {Red, Arrowheads[{{0, 0}, {0.03, 0.6}, {0, 1}}],
Arrow[##]} & /@ Partition[nf[a, c, n], 2, 1]},
ListAnimate[
Table[Plot[{x, f[a, x]}, {x, 0, 2.5}, PlotRange -> Table[{0, 1}, 2],
Epilog -> {Black, PointSize[0.03], Point[{c, f[a, c]}],
Point[roots[a]], Green, PointSize[0.02], Point[roots[a]],
arrows[[1 ;; step]]}, Frame -> True, AspectRatio -> Automatic,
PlotLabel -> Row[{"a = ", a}]], {step, 1, Length[arrows]}],
AnimationRunning -> False]]
, {{c, 0.1, "\!$$\*SubscriptBox[\(x$$, $$0$$]\)"}, 0, 1}, {n,
Range[5]}, {a, 0, 2}]


Inspired by bills

The following modifies the function to $x^3- a x$ and illustrates behaviour for parameter values $a$ 0 to 3. To run clear variables:

f[a_, x_] := x^3 - a x
roots[a_] := {x, x} /. Solve[f[a, x] == x, x]
nf[a_, x_, n_] :=
Catenate[{{##}, {#2, #2}} & @@@
NestList[{#[[2]], f[a, #[[2]]]} &, {x, f[a, x]}, n]]
Manipulate[
Module[{arrows = {Red, Arrowheads[{{0, 0}, {0.03, 0.6}, {0, 1}}],
Arrow[##]} & /@ Partition[nf[a, c, n], 2, 1]},
ListAnimate[
Table[Plot[{x, f[a, x]}, {x, -2, 2},
PlotRange -> Table[{-2, 2}, 2],
Epilog -> {Black, PointSize[0.03], Point[{c, f[a, c]}],
Point[roots[a]], Green, PointSize[0.02], Point[roots[a]],
arrows[[1 ;; step]]}, Frame -> True, AspectRatio -> Automatic,
PlotLabel -> Row[{"a = ", a}]], {step, 1, Length[arrows]}],
AnimationRunning -> False]], {{c, 0.1,
"\!$$\*SubscriptBox[\(x$$, $$0$$]\)"}, -2, 2}, {n, Range[5]}, {a,
0, 3}]


and the bifurcation diagram (100 iterates):

tu = Tuples[{Range[1, 3, 0.01], Range[-1, 1, 0.1]}];
nst[a_, x_, n_] := Nest[f[a, #] &, x, n]
ta = {#1, nst[##, 100]} & @@@ tu;
ListPlot[ta]


or stepping starting value by 0.01:

You can draw the bifurcation diagram pretty straightforwardly. First, define the iteration:

Clear[x, a];
x[a_, n_] := x[a, n] = x[a, n - 1]^3 + a x[a, n - 1];
x[a_, 1] := RandomReal[{-0.5, 0.5}];;


With a little checking, the iteration has interesting dynamics for a between about -3 and 0, hence:

all = Flatten[Table[seq = x[aa, #] & /@ Range[1000];
seq1000 = Union[seq[[900 ;; 1000]]];
{aa, -3., -0.5, 0.001}], 1];
ListPlot[all, PlotRange -> All]


When a is positive, the iteration diverges, as you can see from a linearization argument.

• beautiful...+1. I posted my answer to illustrate why the behaviour in the OP parameter range was causing 'problems' :) – ubpdqn Nov 30 '16 at 2:32

A direct fix of the posted code: change x to #, change Range[0, 1, 0.01] to Range[-1, 1, 0.01] for completeness, change the range of a to be {a, -3, 0, 0.01}, and Catenate the sublists:

ListPlot[Catenate@
ParallelTable[
Thread[{a, Nest[#^3 + a # &, Range[-1, 1, 0.01], 1000]}], {a, -3,
0, 0.01}], PlotStyle -> PointSize[0]]