Advertisement
Guest User

Untitled

a guest
Jul 10th, 2014
217
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.71 KB | None | 0 0
  1. m1 = 1; m2 = 1; r1 = 1; c = 1; k = 1; g = 1;
  2.  
  3. 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)};
  4.  
  5. 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]];
  6.  
  7. abcdata2 = Map[psectF2, {{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.1, 0.01}}];
  8.  
  9. ListPlot[abcdata2, ImageSize -> Medium]
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement