# help to plot Poincaré section for double pendulum

I am reading a book about classical mechanics. In the chapter about chaos, it gives the simplified and scaled equations for double pendulum as

$$\frac{d}{dt}\left[ \begin{matrix} \alpha \\[3mm] l_\alpha \\[3mm] \beta \\[3mm] l_\beta \end{matrix} \right] = \left[ \begin{matrix} 2\frac{l_\alpha - (1+\cos\beta)l_\beta}{3-\cos 2\beta} \\[2mm] -2\sin\alpha - \sin(\alpha + \beta) \\[2mm] 2\frac{-(1+\cos\beta)l_\alpha + (3+2\cos\beta)l_\beta}{3-\cos2\beta}\\[2mm] -\sin(\alpha+\beta) - 2\sin\beta\frac{(l_\alpha-l_\beta)l_\beta}{3-\cos2\beta} + 2\sin2\beta\frac{l_\alpha^2-2(1+\cos\beta)l_\alpha l_\beta + (3+2\cos\beta)l_\beta^2}{(3-\cos2\beta)^2} \end{matrix} \right]$$

I am going to plot the Poincaré section, which in this case is the collection of $(\beta, l_\beta)$ when $\alpha=0$ and $d\alpha/dt >0$. I follow the code found in the example of Mathematica for Poincaré section and write the following code

abc = {a'[t] == 2*(la[t] - (1 + Cos[b[t]])*lb[t])/(3 - Cos[2*b[t]]),
la'[t] == -2*Sin[a[t]] - Sin[a[t] + b[t]],
b'[t] == 2*(-(1 + Cos[b[t]])*la[t] + (3 + 2*Cos[b[t]])*lb[t])/(3 - Cos[2*b[t]]),
lb'[t] == -Sin[a[t] + b[t]] - 2*Sin[b[t]]*((la[t] - lb[t])*lb[t])/(3 - Cos[2*b[t]]) +
2*Sin[2*b[t]]*(la[t]^2 - 2*(1 + Cos[b[t]])*la[t]*lb[t] +
(3 + 2*Cos[b[t]])*lb[t]^2)/(3 - Cos[2*b[t]])^2};

psect[{a0_, la0_, b0_ , lb0_}] :=
Reap[NDSolve[{abc, a[0] == a0, la[0] == la0, b[0] == b0, lb[0] == lb0,
WhenEvent[a[t] == 0 && la'[t]>0, Sow[{b[t], lb[t]}]]}, {},
{t, 0, 1000},  MaxSteps -> ∞]][[-1, 1]]

abcdata = Map[psect, {{0.01, 0.01, 0.01, 0}, {0.01, 0.01, 0.01, 0.01},
{0.01, 0.01, 0.01, -0.01}}];

ListPlot[abcdata, ImageSize -> Medium]


Here a is $\alpha$, b is $\beta$, la is $l_\alpha$ and lb is $l_\beta$. I don't know why, they code has two problems. First the condition "la'[t]>0" doesn't work. Secondly, even I remove the condition "&& la'[t]>0", it sill reports that error code like "NDSolve::deqn: Equation or list of equations expected instead of 0.01 in the first argument". I spent quite a long time to figure out the problem but I still cannot find out the reason.

Here is plot for the Poincaré section from the book I am reading (I redraw with image editor).

• One mistake: You wrote lb[0] = lb0 which is given away by the fact that lb is coloured black, as well as by Mathematica seeing a number where there should be an equation. Use ==. – Szabolcs Jan 9 '14 at 21:32
• Thanks for pointing out the mistake. I correct that but I still don't know how to extract the section when $\alpha=0$ and $d\alpha/dt>0$ – user1285419 Jan 9 '14 at 21:52

WhenEvents[] is a tricky beast (thanks to @Szabolcs for pointing out a better way to write the event detection part)

psect[{a0_, la0_, b0_, lb0_}] :=
Reap[NDSolve[{abc, a[0] == a0, la[0] == la0, b[0] == b0, lb[0] == lb0,
WhenEvent[a[t] + 0*b'[t] + 0*la[t] + 0*lb[t] == 0,
If[la[t] >= 0, Sow[{b[t], lb[t]}]]]},
{a, la, b, lb}, {t, 0, 200}, MaxSteps -> \[Infinity],
MaxStepSize -> .0005]][[2, 1]];
ps = psect /@ Table[{0.01, 0.005, 0.01, i}, {i, -.01, .01, .00125}];
ListLinePlot[#[[FindCurvePath[#][[1]]]] & /@ ps, Mesh -> All, PlotRange -> All]


• This didn't work for me when I tried the same ... When I tried to actually Sow la'[t], I noticed that it came out evaluated ... let me paste your code ... – Szabolcs Jan 9 '14 at 21:55
• OK, I see the trick you did now: you included la' in the event part of WhenEvent to force its evaluation in the action part. This is a bit confusing though. Suppose you'd want la'[t] >= 0 instead of la'[t] > 0. How would you avoid la[t] == 0 && a[t] != 0 being detected? Also, do you also think that what we see is a bug? – Szabolcs Jan 9 '14 at 22:02
• @Szabolcs I found a lot of WhenEvent[] quirks like this one lately. I assume it's just immature (but useful) – Dr. belisarius Jan 9 '14 at 22:04
• This seems to give an incorrect result, because of the reason I mentioned. Try Sowing a[t] and la'[t] instead to see their values. It will detect points where la'[t] is zero up to a Chop but la'[t]>0 is also true. We can modify it slightly to make it work, and we don't even need the If part. The point is to mention la'[t] in the event somewhere but make absolutely sure that it does't make a difference. note that WhenEvent is HoldAll! So we can do this: – Szabolcs Jan 9 '14 at 22:17
• WhenEvent[a[t] + 0*la'[t] == 0 && la'[t] > 0, Sow[{whatever}]]. This will give a correct result and multiplying la'[t] by 0 ensures that it does not make a difference in the event. – Szabolcs Jan 9 '14 at 22:18

Thanks to Szabolcs for pointing out the difference between f and pred in the form

WhenEvent[f == 0 && pred,...]


The function f is allowed to be either a function or its derivative that appears in the system of DEs (for example lb'[t]), but apparently pred is not. Judging from the error messages, only values of t, a[t], b[t], la[t], lb[t] are passed to pred.

The difference, I suppose, is that lb'[t] is not a state variable. Only a[t], b[t], la[t], lb[t] are; those, as well as t, are replaced by their values in WhenEvent. The way around is to use the fact that the DE defines what lb'[t] etc. are in terms of the state variables. We can use the equation in the DE system instead of lb'[t]. We'll also mark both crossings, where lb'[t] changes from negative to positive and vice versa (two WhenEvents), otherwise only the right half of the figure is drawn.

psect[{a0_, la0_, b0_, lb0_}] :=
Reap[NDSolve[{abc, a[0] == a0, la[0] == la0, b[0] == b0, lb[0] == lb0,
WhenEvent[a[t] == 0 && -2*Sin[a[t]] - Sin[a[t] + b[t]] < 0, Sow[{b[t], lb[t]}]],
WhenEvent[a[t] == 0 && -2*Sin[a[t]] - Sin[a[t] + b[t]] > 0, Sow[{b[t], lb[t]}]]},
{a[t], b[t], la[t], lb[t]}, {t, 0, 1000},
MaxSteps -> Infinity]][[-1, 1]]

abcdata = Map[psect,
{{0.01, 0.01, 0.01, 0},
{0.01, 0.01, 0.01, 0.01},
{0.01, 0.01, 0.01, -0.01}}];

ListPlot[abcdata, ImageSize -> Medium]


Using @belisarius' data:

ps = psect /@ Table[{0.01, 0.005, 0.01, i}, {i, -.01, .01, .00125}];
ListPlot[ps, ImageSize -> Medium]


Update

The dimension of the phase space is four, {a[t], b[t], la[t], lb[t]} and the Poincaré section a[t] == 0 is three dimensional. What we see above is a projection onto the {b, lb} plane, and projections can make things appear to cross. However we can plot the iterations of the return map in 3D. We should Sow[{b[t], la[t], lb[t]}] to sow the all nontrivial coordinates of the section:

psect3[{a0_, la0_, b0_, lb0_}] :=
Reap[NDSolve[{abc, a[0] == a0, la[0] == la0, b[0] == b0, lb[0] == lb0,
WhenEvent[a[t] == 0 && -2*Sin[a[t]] - Sin[a[t] + b[t]] < 0,
Sow[{b[t], la[t], lb[t]}]],
WhenEvent[a[t] == 0 && -2*Sin[a[t]] - Sin[a[t] + b[t]] > 0,
Sow[{b[t], la[t], lb[t]}]]},
{a[t], b[t], la[t], lb[t]}, {t, 0, 1000},
MaxSteps -> Infinity]][[-1, 1]]


This shows a 3D image from a direction that nearly matches the image above:

ps3 = psect3 /@ Table[{0.01, 0.01, 0.01, i}, {i, -.01, .01, .00125}];
ListPointPlot3D[ps3, ImageSize -> Medium]


We could also show a different angle & coloring:

Graphics3D[
MapIndexed[{Hue[First[#2]/Length[ps], 0.8, 0.8], Point[#]} &, ps],
ImageSize -> Medium, Axes -> True]


I would say there is nothing wrong with WhenEvent. However, the documentation might be a bit more complete. The reason I gave for the assertion that lb'[t] is not a valid variable for pred in WhenEvent is not clearly stated in the docs, but the docs are long and I may not have read carefully enough. However, it makes sense from a mathematical point of view. The integrator/solver is not going to be keeping track of derivatives in a first-order system. While Mathematica could recursively determine any order derivative from the differential equation, it is doubtful that it does.

The confusing thing is that some of the examples in the docs are of second-order equations. For those equations, the first derivative is a state variable and it may be used in WhenEvent.

• But if you do mention lb'[t] in the event part of WhenEvent (and not in the predicate coming after &&) then it does work and the value of lb'[t] suddenly becomes available. (See the comments below belisarius's answer.) This is why I thought it was a bug. If WhenEvent knows that lb'[t] needs to be kept track of when it's mentioned in the event, it should also do it when it's mentioned in the predicate. (I'm referring to the WhenEvent[f==0 && pred, ...] syntax mentioned in the docs.) Nice colour scheme in your list plot btw ;-) and +1 for actually generating the figure. – Szabolcs Jan 10 '14 at 2:45
• @Szabolcs Right. No wonder WhenEvent is so confusing. It seems clear (?) that the second argument to && is a function of only {t, a, b, la, lb}, judging from the error messages. – Michael E2 Jan 10 '14 at 2:55
• Thanks a lot for the neat code. I learn something from this also. But I just don't understand the result why there is some curves crossing while the book say it is not suppose to cross. – user1285419 Jan 10 '14 at 3:11
• @user1285419 I cannot explain the image in the book. It's close to what is produced above. Artistic license, perhaps? – Michael E2 Jan 10 '14 at 17:07
• Hi Michael, thanks for the update. The plots doesn't look about the same but there are one important thing not consistent. The image I attached in the post is about the same to the book. We see that there are two small islands, corresponding to the fix points, in the top and bottom. And they are not symmetric. In your plot, the plots looks identical and symmetric in the top and bottom. – user1285419 Jan 10 '14 at 17:27

You have one mistake: you wrote lb[0] = lb0 instead of lb[0] == lb0.

Once correcting that there's still the problem with including && la'[t]>0. I don't know what causes this. As far as I can tell, this is valid syntax (it's documented). I would contact WRI support and ask if this is a bug or an unsupported feature ... When you get a reply, it'd be nice if you could post here.

Fortunately I found one workaround:

Add lad[t] == la'[t] to the equations and write the condition on lad instead of la'. It'll work, and the results seem to be correct (as compared with a plot).

Here's some shorter code reproducing the same problem:

NDSolve[{a'[t] == b[t], b'[t] == -a[t], a[0] == 1, b[0] == 0,
WhenEvent[a[t] == 0.5 && a'[t] > 0, Sow[t]]},

{a, b}, {t, 0, 10}]
`