1. Not finding help here? Sign up for a free 30min tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Poincaré Sections of the double pendulum with Mathematica

  1. Jul 10, 2014 #1
    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 lenght s. The equations of motions are:

    [itex]
    \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)}
    [/itex]

    [itex]
    \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)]}
    [/itex]

    [itex]
    \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)]}
    [/itex]

    [itex]
    \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}
    [/itex]

    [itex]
    \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}
    [/itex]

    [itex]
    \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)
    [/itex]

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

    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 :)
     
    Last edited: Jul 10, 2014
  2. jcsd
  3. Jul 11, 2014 #2
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Poincaré Sections of the double pendulum with Mathematica
  1. Double Pendulum ! (Replies: 1)

Loading...