- #1
Jay Climaco
- 1
- 0
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:
[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 :)
[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: