SHARE
TWEET

Double Pendulum

Matthen Jul 23rd, 2011 1,846 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. (* Drawing of diagram *)
  2. DrawPendula[\[Theta]1_, \[Theta]2_, m1_, m2_, l1_, l2_] :=
  3.   Module[{radius1, radius2}, (
  4.     radius1 = 0.2 (m1/Max[m1, m2])^(1/3);
  5.     radius2 = 0.2 (m2/Max[m1, m2])^(1/3);
  6.     Graphics[{
  7.       PointSize[Large],
  8.       Blue, Point[{0, 0}],
  9.       Orange,
  10.       Line[{{0,
  11.          0}, {l1 Sin[\[Theta]1], -l1 Cos[\[Theta]1]}, {l1 \
  12. Sin[\[Theta]1] + l2 Sin[\[Theta]2], -l1 Cos[\[Theta]1] -
  13.           l2 Cos[\[Theta]2]}}],
  14.       Red, Disk[{l1 Sin[\[Theta]1], -l1 Cos[\[Theta]1]}, radius1],
  15.       Disk[{l1 Sin[\[Theta]1] +
  16.          l2 Sin[\[Theta]2], -l1 Cos[\[Theta]1] - l2 Cos[\[Theta]2]},
  17.        radius2]
  18.       }, PlotRange -> Table[{-1 - l1 - l2, 1 + l1 + l2}, {2}],
  19.      ImageSize -> {500}]
  20.     )];
  21. (* Interface *)
  22. doublepend[\[Theta]1_, \[Theta]2_, d\[Theta]1_, d\[Theta]2_, m1_, m2_,
  23.     l1_, l2_, g_, \[Alpha]_, ta_] :=
  24.   Module[{eqns, soln, A, B, plot, \[Alpha]2},
  25.    \[Alpha]2 = \[Alpha] B'[t];
  26.    eqns = {m2 l2 B''[t] + m2 l1 A''[t] Cos[A[t] - B[t]] -
  27.        m2 l1 (A'[t]^2) Sin[A[t] - B[t]] +
  28.        m2 g Sin[B[t]] + \[Alpha]2 m2 l2 ==
  29.       0, (m1 + m2) l1 A''[t] + m2 l2 B''[t] Cos[A[t] - B[t]] +
  30.        m2 l2 (B'[t]^2) Sin[A[t] - B[t]] +
  31.        g (m1 + m2) Sin[A[t]] + \[Alpha] A'[t] (m1 + m2) == 0,
  32.      A[0] == \[Theta]1, A'[0] == d\[Theta]1, B[0] == \[Theta]2,
  33.      B'[0] == d\[Theta]2};
  34.    soln =
  35.     NDSolve[eqns, {A, B}, {t, ta - ta/5 - 0.6, ta},
  36.      MaxSteps -> Infinity, PrecisionGoal -> 10];
  37.    If[ta < 25,
  38.     plot =
  39.      ParametricPlot[{{l1 Sin[A[t]], -l1 Cos[A[t]]}, {l1 Sin[A[t]] +
  40.           l2 Sin[B[t]], -l1 Cos[A[t]] - l2 Cos[B[t]]}} /. soln, {t,
  41.        ta - ta/5 - 0.5, ta}];
  42.     Show[First[DrawPendula[A[ta], B[ta], m1, m2, l1, l2] /. soln],
  43.      plot],
  44.     Graphics[Text[Style["(Solution now unstable, probs)", Large]],
  45.      ImageSize -> {500}]
  46.     ]
  47.    ];
  48.  
  49. Manipulate[
  50.  doublepend[\[Theta]1, \[Theta]2, d\[Theta]1, d\[Theta]2, m1, m2, l1,
  51.   l2, g, \[Alpha], ta]
  52.  ,
  53.  {{\[Theta]1, 7 Pi/8,
  54.    "Initial \!\(\*SubscriptBox[\"\[Theta]\", \"1\"]\)"}, 0, 2 Pi},
  55.  {{\[Theta]2, Pi/3,
  56.    "Initial \!\(\*SubscriptBox[\"\[Theta]\", \"2\"]\)"}, 0, 2 Pi},
  57.  {{d\[Theta]1, 2,
  58.    "Initial \!\(\*SubscriptBox[\"\[Theta]\", \"1\"]\)'"}, -5, 5},
  59.  {{d\[Theta]2, 1,
  60.    "Initial \!\(\*SubscriptBox[\"\[Theta]\", \"2\"]\)'"}, -5, 5},
  61.  {{m1, 1, "\!\(\*SubscriptBox[\"m\", \"1\"]\)"}, 0.1, 5},
  62.  {{m2, 2, "\!\(\*SubscriptBox[\"m\", \"2\"]\)"}, 0.1, 5},
  63.  {{l1, 1.5, "\!\(\*SubscriptBox[\"l\", \"1\"]\)"}, 0.1, 5},
  64.  {{l2, 1, "\!\(\*SubscriptBox[\"l\", \"2\"]\)"}, 0.1, 5},
  65.  {{g, 10, "gravity"}, 0, 20},
  66.  {{\[Alpha], 0.1, "friction"}, 0, 0.5},
  67.  {{ta, 0.001, "animate"}, 0.001, Infinity, 0.05,
  68.   ControlType -> Trigger}
  69.  ]
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top