1
$\begingroup$

I am working on double physical pendulum with fixed freedom of one bar. Angle te2 (angle of lower bar) can not be bigger than angle te1 (angle of upper bar) How I can make WhenEvent to get that? Thank you for help.

Clear[poincarePend];
m1 = 1; m2 = 1.2; g = 9.81; L1 = 0.5; L2 = 0.6; L = 1.2;
te10 = -6 \[Pi]/12; dte10 = 1.0; te20 = -\[Pi]/2; dte20 = 2.0; tmax = 100; 
steps = 1000000; wd = \[Pi]/12;
poincarePend[m1_, m2_, g_, L1_, L2_, L_, te10_, dte10_, te20_, dte20_,
tmax_, steps_, wd_] :=  temp1 = NDSolve[{dte1[t] - te1'[t] == 0, dte2[t] - 
te2'[t] == 0, ((1/12)*2 L1*2 L1*m1 + m1*L1*L1 + m2*L*L)*dte1'[t] + 
   m2*L2*L*dte2'[t]*Cos[te2[t] - te1[t]] - 
   m2*L*L2*dte2[t]^2*Sin[te2[t] - te1[t]] + 
   g*(m1*L1 + m2*L)*Sin[te1[t]] == 
  0, ((1/12)*2 L2*2 L2*m2 + m2*L2*L2)*dte2'[t] + 
   m2*L1*L*dte1'[t]*Cos[te2[t] - te1[t]] + 
   m2*L*L2*dte1[t]^2*Sin[te2[t] - te1[t]] + g*m2*Sin[te2[t]] == 0,
  dte1[0] == dte10, te1[0] == te10, dte2[0] == dte20, 
 te2[0] == te20, 
 WhenEvent[te2[t] < te1[t], te2[t] -> te1[t]; 
  dte2[t] -> 0 ] }, {te1, dte1, te2, dte2}, {t, 0, tmax}, 
MaxSteps -> steps];
poincarePend[m1, m2, g, L1, L2, L, te10, dte10, te20, dte20, tmax, 
steps, wd];
temp1 = temp1[[1]];
Animate[Graphics[{{PointSize[.01], {Red, 
Point[{x1[t], y1[t]}]}, {Blue, Point[{x2[t], y2[t]}]}, 
Line[{{0, 0}, {x1[t], y1[t]}, {x2[t], y2[t]}}]} /. temp1, {Gray, 
Line[Map[Function[Evaluate[{x2[#], y2[#]} /. temp1]], 
 Range[0, t, 0.02]]]}}, PlotRange -> {{-1.2, 1.2}, {-1.5, .1}}, 
 Axes -> True, Ticks -> False, ImageSize -> 700], {t, 5, 15, .01}, 
 SaveDefinitions -> True]
$\endgroup$
  • 1
    $\begingroup$ Coordinates {x1[t], y1[t]}, {x2[t], y2[t]} not defined $\endgroup$ – Alex Trounev Nov 14 '18 at 18:04
2
$\begingroup$

Coordinates are not defined. If we remove the excess from the code, correct a couple of typos and add coordinates, then at the output we get

m1 = 1; m2 = 1.2; g = 9.81; L1 = 0.5; L2 = 0.6; L = 1.2;
te10 = -6 \[Pi]/12; dte10 = 1.0; te20 = -\[Pi]/
  2; dte20 = 2.0; tmax = 100;
steps = 1000000; wd = \[Pi]/12;
temp1 = NDSolve[{dte1[t] - te1'[t] == 0, 
    dte2[t] - te2'[t] == 
     0, ((1/12)*2 L1*2 L1*m1 + m1*L1*L1 + m2*L*L)*dte1'[t] + 
      m2*L2*L*dte2'[t]*Cos[te2[t] - te1[t]] - 
      m2*L*L2*dte2[t]^2*Sin[te2[t] - te1[t]] + 
      g*(m1*L1 + m2*L)*Sin[te1[t]] == 
     0, ((1/12)*2 L2*2 L2*m2 + m2*L2*L2)*dte2'[t] + 
      m2*L1*L*dte1'[t]*Cos[te2[t] - te1[t]] + 
      m2*L*L2*dte1[t]^2*Sin[te2[t] - te1[t]] + g*m2*Sin[te2[t]] == 0, 
    dte1[0] == dte10, te1[0] == te10, dte2[0] == dte20, 
    te2[0] == te20}, {te1, dte1, te2, dte2}, {t, 0, tmax}, 
   MaxSteps -> steps];
x1[t_] := L1*Sin[te1[t]] /. temp1
y1[t_] := Cos[te1[t]]*L1 /. temp1
x2[t_] := x1[t] + Sin[te2[t]]*L2 /. temp1
y2[t_] := y1[t] + Cos[te2[t]]*L2 /. temp1

{Plot[{x1[t], y1[t]}, {t, 0, 15}, PlotLegends -> "Expressions"], 
 Plot[{x2[t], y2[t]}, {t, 0, 15}, PlotLegends -> "Expressions"]}

fig1

ListAnimate[
  Table[Graphics[{{PointSize[.04], Red, 
      Point[-Flatten[{x1[t], y1[t]}]]}, {PointSize[.05], Blue, 
      Point[Flatten[-{x2[t], y2[t]}]]}, {Thick, Orange, 
      Line[{{0, 0}, Flatten[-{x1[t], y1[t]}], 
        Flatten[-{x2[t], y2[t]}]}]}}, 
    PlotRange -> {{-1.2, 1.2}, {-1.2, .6}}, Axes -> True, 
    Ticks -> False, ImageSize -> 300], {t, 0, 15, .05}]]

fig2

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.