A short post inspired by

this question. Please also see this literature:

1 ,

2 ,

3. Duffing oscillator is one of the simplest system that exhibit chaos. It is a periodically forced oscillator with a nonlinear elasticity. We will use the following from:

eq = x''[t] + ? x'[t] - x[t] + x[t]^3 == ? Cos[ t];

eq // TraditionalForm

The unforced undamped case with gamma and delta equal to zero is a Hamiltonian system with potential energy representing a double well potential. So basically a general case is when a driving force bounces a point between to potential minima with friction present. Usually it is associated with a model of a periodically forced steel beam which is deflected toward the two magnets:

When the system in non-chaotic regime it is synchronized perfectly with the periodic force. In the phase space of x and x' the trajectory of the system is passing always through the same point and is a closed lopped trajectory over a period 2 PI of the force. When we are in chaotic regie the trajectory will split and visit many close points during different periods - with points forming an attractor which can be visualized via a Poincaré section:

Imagine you will move the Poincaré section in above picture around the complete circle - how would attractor change shape? It is inefficient to run NDSolve for every new Poincaré section and we better get all the sections within one single run of NDSolve. To do that we need 2 things:

1) Specify events to track by NDSolve with WhenEvent

2) Differentiate between events so at the end they won't pile up in mixed data

Here is a simple case of tracking 6 events distant at Pi/3 - all events on one figure colored differently:

evs = Mod[t, 2 \[Pi]] == # & /@ (Range[0, 2 Pi - #, #] &@(2 Pi/6))

Out[] = {Mod[t, 2 ?] == 0, Mod[t, 2 ?] == ?/3, Mod[t, 2 ?] == (2 ?)/3,

Mod[t, 2 ?] == ?, Mod[t, 2 ?] == (4 ?)/3, Mod[t, 2 ?] == (5 ?)/3}

Note that Sow function has a 2nd argument - which is a tag that allows to Reap events separately.

data = Block[{? = 0.15, ? = 0.3},

Reap[NDSolve[{eq, x[0] == 0, x'[0] == 0, WhenEvent[Evaluate@evs,

Sow[{x[t], x'[t]}, Round[100 Mod[t, 2 ?]]]]}, {}, {t, 0, 200000}, MaxSteps -> ?]]];

ListPlot[data[[2]], PlotStyle -> Directive[Opacity[.5], PointSize[.002]],

PlotRange -> {1.8 {-1, 1}, 1.3 {-1, 1}}, AspectRatio -> 1, Frame -> True, ImageSize -> 450]

And here is a more complicated cases of 50 events shown as animation that runs while we move Poincaré section around full circle:

evs = Mod[t, 2 ?] == # & /@ (Range[0, 2 Pi - #, #] &@(2 Pi/50));

data = Block[{? = 0.15, ? = 0.3},

Reap[NDSolve[{eq, x[0] == 0, x'[0] == 0, WhenEvent[Evaluate@evs,

Sow[{x[t], x'[t]}, Round[100 Mod[t, 2 ?]]]]}, {}, {t, 0, 200000}, MaxSteps -> ?]]];

Manipulate[

ListPlot[data[[2]][[k]], PlotStyle -> PointSize[0],

PlotRange -> {1.8 {-1, 1}, 1.3 {-1, 1}}, AspectRatio -> 1,

Frame -> True, ImageSize -> 450], {k, 1, 50, 1}]