# PoincareSection for a driven damped pendulum is not generating a Poincaré section at all, why?

So I have the general code from the PoincareSection documentation that is changed up for a Driven Damped Pendulum:

coupledDiffEq =
{ω'[t] == -(1/q) ω[t] - Sin[θ[t]] +
g*Cos[ϕ[t]],
θ'[t] == ω[t],
ϕ'[t] ==  Drive};

data = Block[{q = 3.9, g = 1.5, Drive = 1},
Reap[NDSolve[{coupledDiffEq,
θ == 0, ω == 0, ϕ == 2*Pi,
WhenEvent[Mod[ϕ[t], (2*Pi)] == 0,
Sow[{θ[t], ω[t]}]]}, {}, {t, 0, 100000},
MaxSteps -> ∞]]][[-1, 1]];
ListLinePlot[data, PlotRange -> {{-3, 3}, All}]


Then I ListPlot data and it comes out as not a PoincareSection. Often it looks like a phase plot more than anything else. I have attempted to change the

Mod[ϕ[t], DriveFrequency]


where DriveFrequency should be Drive / (2*Pi), but that often does nothing. The system of equations evaluates and produces points, they just are random though. How do I get a nice Poincaré section with this? I am getting this when I should be getting • Yes, when I use the Parameters (q = 3.9, g =1.5, Drive = .6667) I get something that looks more like a phase diagram. Imager Link to result: imgur.com/8natnP0 Nov 7 '15 at 23:12

Update

J.M. has some good comments, so I'm going to implement them.

coupledDiffEq = {θ''[t] == -(1/q) θ'[t] -
Sin[θ[t]] + g*Cos[ϕ[t]], ϕ'[t] == drive};

data = Block[{q = 3.9, g = 1.5, drive = 2/3}
, First@Last@Reap@NDSolve[{
coupledDiffEq, θ == 0, θ' == 0, ϕ == 2 π
, WhenEvent[Mod[ϕ[t], 2 π] == 0, Sow[{Mod[θ[t], 2 π, -π], θ'[t]}]]
}, {}, {t, 0, 100000}, MaxSteps -> ∞]
];
ListPlot[Drop[data, 100]]


I dropped the first 100 points to make sure that we aren't capturing any of the transients.

I played around with your code, and couldn't get it to do what I wanted, so I went this route instead. Note that you need to move θ[t] back into the interval -π to π in order to get the attractor, and you should not use ListLinePlot, because you will get a mess.

ClearAll[coupledDiffEq, ω, q, θ, g, ϕ, drive]
coupledDiffEq := {ω'[t] == -(1/q) ω[t] - Sin[θ[t]] + g*Cos[ϕ[t]], θ'[t] == ω[t], ϕ'[t] == drive};

sols = Block[{q = 3.9, g = 1.5, drive = 2/3}
, First@NDSolve[
{coupledDiffEq, θ == 0, ω == 0, ϕ == 2*π}
, {θ[t], ω[t]}
, {t, 0, 100000}
, MaxSteps -> ∞
]
];

data = Table[{Mod[θ[t] + π, 2 π] - π, ω[t]} /. sols, {t, 20*(2 π)/(2/3), 100000, (2 π)/(2/3)}];
ListPlot[data, PlotRange -> All] Update

This code should work too, more along the lines of what you did:

coupledDiffEq = {ω'[t] == -(1/q) ω[t] - Sin[θ[t]] + g*Cos[ϕ[t]], θ'[t] == ω[t], ϕ'[t] == drive};

sols = Block[{q = 3.9, g = 1.5, drive = 2/3}
, First@Last@Reap@NDSolve[{coupledDiffEq, θ == 0, ω == 0, ϕ == 2 π
, WhenEvent[Mod[ϕ[t], 2 π] == 0, Sow[{Mod[θ[t] + π, 2 π] - π, ω[t]}]]
}, {θ[t], ω[t]}, {t, 0, 100000}, MaxSteps -> ∞]];

ListPlot[sols] • Two notes: 1. As I noted in this answer, you can pass an empty list as the second argument of NDSolve[] so that it doesn't save the InterpolatingFunction[] objects, potentially saving memory. OP seems to have tried this, but could not make it work. 2. Mod[θ[t] + π, 2 π] - π is slightly more compactly done as Mod[θ[t], 2 π, -π]. Nov 8 '15 at 8:32
• You are amazing! Its quite interesting how you have to play around with Sow[...] and θ[t] to get this to work! I really appreciate your help! Nov 8 '15 at 14:20
• @Escap3faith, BTW, you don't really need ω for your system (assuming you only used it as a temporary intermediate); NDSolve[] can handle θ''[t] == -θ'[t]/q - Sin[θ[t]] + g Cos[ϕ[t]] perfectly, and your Sow[] becomes Sow[{Mod[θ[t], 2 π, -π], θ'[t]}]. Nov 8 '15 at 17:28
• @J.M. Thanks for comments and links! I actually couldn't think of a good reason why not to actually solve for the functions, but not returning a memory-intensive InterpolatingFunction is a good one. I will modify the post with your comments, I think, to make the answer more worthwhile. Nov 8 '15 at 17:38
• @J.M. I don't think I do, but I was reading a paper and they broke it down into 3 equations due to the necessity of ϕ[t], and I just followed out of frustration for the project. I'm 100% positive I don't need all 3 coupled equations because Mathematica has the power to deal with them, but rather just the 2. Paper: thphys.uni-heidelberg.de/~gasenzer/… Nov 8 '15 at 19:50