# Problems making bifurcation diagram for damped driven pendulum

tl:dr Bifurcation diagram has weird gaps and displaced points,

We made some code in mathematica to plot a bifurcation diagram for the driven damped pendulum. We worked out the position for a given value of the controlling variable g using NDsolve. After that, we took the points past the middle to eliminate transients.

Using a while loop, we iterated over g. Plotting this gave us

If you look at the period doubling (0.62 - 0.68) you see that many points are missing in the plot. There are random gaps. And in many places it looks like some of the points have been displaced.

Trying to plot lower values of g gave a single line which increased the position gradually till a big drop to the negatives. It was discontinuous. Plotting the velocity does the same thing.

We have no clue why any of this is happening.

Here's the code in raw text.

a = 0.61;
xList = {};
n = 10;

While[a <= 0.75,
{

ClearAll[sols, Data];

g = a;

coupledDiffEq = {ω'[t] == -(1/q) ω[t] -

Sin[θ[t]] + g*Cos[ϕ[t]], θ'[t] == ω[

t], ϕ'[t] == drive};

sols = Block[{q = 3.9, drive = 2/3},

First@Last@

Reap@NDSolve[{coupledDiffEq, θ[0] == 0, ω[0] ==

0, ϕ[1] == 2 π,

WhenEvent[Mod[ϕ[t], 2 π] == 0,

Sow[{Mod[θ[t] + π, 2 π] - π, ω[

t]}]]}, {θ[t], ω[t]}, {t, 0, 20000},

MaxSteps -> ∞]];

Data = Table[

sols[[i, 1]], {i, Length[sols]/2 + Mod[Length[sols], 2]/2,

Length[sols]}];

list = {g};

alist = {};

For[i = 1, i <= n, i++,

list = Append[list, Data[[i]] ];

alist = Append[alist, list];

list = {g};

];

xList = Join[xList, alist];

Print[a],

a = a + 0.0005}

]

• I think there is nothing wrong with your result. Look at the bottom of this page thphys.uni-heidelberg.de/~gasenzer/… – ZaMoC Apr 2 '17 at 20:39
• Ok, I see our graph is not dissimilar to the ones I've seen in the replies. I don't understand why that's going on. Shouldn't the line at the beginning, that is before chaos, be continuous? – Mr insert name here R Apr 3 '17 at 20:38
• there is chaos in the line, too! and everywhere around it. and it's beautiful. – ZaMoC Apr 3 '17 at 20:48

This is some old (so could be made more elegant and efficient) code of mine that I will just drop here in case anyone finds it useful. Note the timings - one might want to consider a bigger step first.

q = 2;
ωD = 2/3;
step = 0.001;
tmax = 1000;
sol = Table[
y[t] /. NDSolve[{ x'[t] == y[t],
y'[t] == -1/q y[t] - Sin[x[t]] + g Cos[z[t]],
z'[t] == ωD, x[0] == 1, y[0] == 1, z[0] == 0}, {x[t],
y[t], z[t]}, {t, 100, tmax}, MaxSteps -> Infinity], {g, 1, 1.5,
step}]; // AbsoluteTiming


{23.0122, Null}

bif = Table[
Table[sol[[i, 1]], {t, 100, tmax, 3 π}], {i, 1,
Length[sol]}]; // AbsoluteTiming


{78.7297, Null}

gs = Table[i, {i, 1, 1.5, step}];
data = Table[
Table[{gs[[j]], bif[[j, i]]}, {i, 1, Length[bif[[j]]]}], {j, 1,
Length[bif]}];
ListPlot[Flatten[data, 1], PlotRange -> All,
PlotStyle -> {Black, PointSize[Tiny]},
AxesLabel -> {g, OverDot[y][t]}, ImageSize -> 600]