help to plot Poincaré section for double pendulum

Question

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).

Answer 1

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]

Answer 2

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.

Answer 3

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}]