Poincaré Sections of the double pendulum with Mathematica

[Jay Climaco]

I need to plot a Poincaré map for a double pendulum where the string connecting one mass to the other is elastic with elasticity constant k and rest length s. The equations of motions are:

$\dot{\theta}_1= \frac{p_{\theta_1}}{m_1 r_1^2} - \frac{p_{\theta_2}}{m_1 r_1 r_2} \cos{(\theta _1- \theta_2)} + \frac{p_{r_2}}{m_1 r_1} \sin{(\theta _1- \theta_2)}$

$\dot{\theta}_2= \frac{p_{\theta_2}}{m_2 r_2^2} + \frac{p_{\theta_2}}{2 m_1 r_2^2} - \frac{ p_{\theta_1}}{m_1 r_1 r_2} \cos{(\theta _1- \theta_2)} + \frac{p_{\theta_2}}{2 m_1 r_2^2} \cos{[2 (\theta _1- \theta_2)]} - \frac{p_{r_2}}{2 m_1 r_1 r_2} \sin{[2 (\theta _1- \theta_2)]}$

$\dot{r}_2= \frac{p_{r_2}}{m_2} + \frac{p_{r_2}}{2 m_1} - \frac{p_{r_2}}{2 m_1} \cos{[2 (\theta _1- \theta_2)]} + \frac{ p_{\theta_1}}{ m_1 r_1} \sin{(\theta _1- \theta_2)} - \frac{p_{\theta_2}}{2 m_1 r_2} \sin{[2 (\theta _1- \theta_2)]}$

$\dot{p_{\theta_1}}= -\frac{p_{\theta_1} p_{\theta_2}}{m_1 r_1 r_2} \sin{(\theta _1- \theta_2)} + \frac{(p_{\theta_2}^2 - r_2^2 p_{r_2}^2)}{2 m_1 r_2^2} \sin{[2 (\theta _1- \theta_2)]} - \frac{p_{\theta_1} p_{r_2}}{m_1 r_1} \cos{(\theta _1- \theta_2)} + \frac{p_{\theta_2} p_{r_2}}{m_1 r_2} \cos{[2 (\theta _1- \theta_2)]} - g r_1 (m_1 + m_2) \sin{\theta_1}$

$\dot{p_{\theta_2}}= \frac{p_{\theta_1} p_{\theta_2}}{m_1 r_1 r_2} \sin{(\theta _1- \theta_2)} - \frac{(p_{\theta_2}^2 - r_2^2 p_{r_2}^2)}{2 m_1 r_2^2} \sin{[2 (\theta _1- \theta_2)]} + \frac{p_{\theta_1} p_{r_2}}{m_1 r_1} \cos{(\theta _1- \theta_2)} - \frac{p_{\theta_2} p_{r_2}}{m_1 r_2} \cos{[2 (\theta _1- \theta_2)]} - g m_2 r_2 \sin{\theta_2}$

$\dot{p_{r_2}}= \frac{p_{\theta_2}^2}{m_2 r_2^3} + \frac{p_{\theta_2}^2}{2 m_1 r_2^3} - \frac{p_{\theta_1} p_{\theta_2}}{m_1 r_1 r_2^2} \cos{(\theta _1- \theta_2)} + \frac{p_{\theta_2}^2}{2 m_1 r_2^3} \cos{[2 (\theta _1- \theta_2)]} + g m_2 \cos{\theta_2} - \frac{p_{\theta_2} p_{r_2}}{2 m_1 r_2^2} \sin{[2 (\theta _1- \theta_2)]} - k(r_2-s)$

So, I tried the following code: (had to change some variables because Mathematica woudn't accept my original choices. The new variables are $a= \theta_1;\ A=p_{\theta_1};\ b=\theta_2;\ B=p_{\theta_2};\ l=r_2;\ L=p_{\theta_2},c=s)$

m1 = 1;
m2 = 1;
r1 = 1;
c = 1;
k = 1;
g = 1;

abc2 = {a'[t] == A[t]/(m1 r1^2) - (B[t]/(m1 r1 l[t]))*Cos[a[t] - b[t]] + (L[t]/(m1 r1))*Sin[a[t] - b[t]], b'[t] == B[t]/(m2 (l[t])^2) + B[t]/(2 m1 (l[t])^2) - (A[t]/(m1 r1 l[t]))*Cos[a[t] - b[t]] + (B[t]/(2 m1 (l[t])^2))*Cos[2 (a[t] - b[t])] - (L[t]/(2 m1 r1 l[t]))*Sin[2 (a[t] - b[t])], l'[t] == L[t]/(m2) + L[t]/(2 m1) - (L[t]/(2 m1))*Cos[2 (a[t] - b[t])] + (A[t]/(m1 r1))*Sin[a[t] - b[t]] - (B[t]/(2 m1 l[t]))*Sin[2 (a[t] - b[t])], A'[t] == -((A[t] B[t])/(m1 r1 l[t]))*Sin[a[t] -b[t]] + (((B[t])^2 - (l[t])^2 (L[t])^2)/(2 m1 (l[t])^2))*Sin[2 (a[t] - b[t])] - ((A[t] L[t])/(m1 r1))*Cos[a[t] - b[t]] + ((B[t] L[t])/(m1 l[t]))*Cos[2 (a[t] - b[t])] - g r1 (m1 + m2)*Sin[a[t]], B'[t] == ((A[t] B[t])/(m1 r1 l[t]))*Sin[a[t] - b[t]] - (((B[t])^2 - (l[t])^2 (L[t])^2)/(2 m1 (l[t])^2))*Sin[2 (a[t] - b[t])] + ((A[t] L[t])/(m1 r1))*Cos[a[t] - b[t]] - ((B[t] L[t])/(m1 l[t]))*Cos[2 (a[t] - b[t])] - g m2 l[t]*Sin[b[t]], L'[t] == ((B[t])^2/(m2 (l[t])^3)) + ((B[t])^2/(2 m1 (l[t])^3)) - ((A[t] B[t])/(m1 r1 (l[t])^2))*Cos[a[t] - b[t]] + ((B[t])^2/(2 m1 (l[t])^3))*Cos[2 (a[t] - b[t])] + g m2*Cos[b[t]] - ((B[t] L[t])/(2 m1 (l[t])^2))*Sin[2 (a[t] - b[t])] - k (l[t] - c)};

psectF2[{a0_, A0_, b0_, B0_, l0_, L0_}] := Reap[NDSolve[{abc2, a[0] == a0, A[0] == A0, b[0] == b0, B[0] == B0, l[0] == l0, L[0] == L0, WhenEvent[a[t] == 0 && l[t] == c && L[t] > 0 && A[t] > 0, Sow[{b[t], B[t]}]], WhenEvent[a[t] == 0 && l[t] == c && L[t] > 0 && A > 0, Sow[{b[t], B[t]}]]}, {a[t], b[t], A[t], B[t], l[t], L[t]}, {t, 0, 1000}, MaxSteps -> Infinity]][[-1, 1]];

abcdata2 = Map[psectF2, {{0.01, 0.01, 0.01, 0.01, 0.01,0.01}, {0.01, 0.01, 0.01, 0.05, 0.01,0.01}, {0.01, 0.01, 0.01, 0.1, 0.01,0.01}}];

ListPlot[abcdata2, ImageSize -> Medium]

The following problems appears when I try to run abcdata2:

Part::partw: "Part 1 of {} does not exist."
General::stop: "Further output of Part::partw will be suppressed during this calculation."

I think the problem is at the WhenEvent command, but I'm not sure. So please, if anyone could help me correct this problem, I'd appreciate :)

 

[Bill Simpson]

1: Notice you have an A>0 in your second WhenEvent. That should be A[t]>0. I found that by taking much of this apart and looking at what happens when I evaluate your various functions one at a time.

2: Your error message does not appear to be because of your WhenEvent. I believe it is because of the [[-1,1]] that you have tacked onto the end of the Reap to try to extract individual items. I'm fiddling with this and haven't tracked down exactly when in the process that happens. If you eliminate the subscripting and don't immediately try to ListPlot the result and just look at exactly what is being returned by psectF2 then I think you can track this down. It looks like it is the way you are trying to extract and use one of the resulting InterpolatingFunction. Try building a Table using the result returned. That often makes it obvious what the error is when Plot just gives you an empty plot, it knows why it didn't display anything and it isn't telling.