Problems making bifurcation diagram for damped driven pendulum

Question

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}

 ]

Answer 1

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]