1
$\begingroup$

enter image description hereI want to make an animation of a double spring spherical pendulum. The first mass is connected to the ceiling with a spring, and the second mass is attached to the first one with a spring too. The generalised coordinates are the displacements of springs and four angles. Below is my code. It doesn't work in the sense that there somehow appears infinite expression 1/0.

m1 = 1; m2 = 1; l1 = 1; l2 = 1; k1 = 5; k2 = 5; g = 9.81;

x1[t_] = (l1 + s1[t]) Sin[\[Theta]1[t]] Sin[\[Phi]1[t]]
y1[t_] = (l1 + s1[t]) Sin[\[Theta]1[t]] Cos[\[Phi]1[t]]
z1[t_] = -(l1 + s1[t]) Cos[\[Theta]1[t]]
x2[t_] = (l1 + s1[t]) Sin[\[Theta]1[t]] Sin[\[Phi]1[t]] + (l2 + 
     s2[t]) Sin[\[Theta]2[t]] Sin[\[Phi]2[t]]
y2[t_] = (l1 + s1[t]) Sin[\[Theta]1[t]] Cos[\[Phi]1[t]] + (l2 + 
     s2[t]) Sin[\[Theta]2[t]] Cos[\[Phi]2[t]]
z2[t_] = -(l1 + s1[t]) Cos[\[Theta]1[t]] - (l2 + s2[t]) Cos[\[Theta]2[
     t]]

x1dot[t_] = D[x1[t], t]
y1dot[t_] = D[y1[t], t]
z1dot[t_] = D[z1[t], t]
x2dot[t_] = D[x2[t], t]
y2dot[t_] = D[y2[t], t]
z2dot[t_] = D[z2[t], t]

Kin[t_] := 
 m1 (x1dot[t]^2 + y1dot[t]^2 + z1dot[t]^2 + s1'[t]^2)/2 + 
  m2 (x2dot[t]^2 + y2dot[t]^2 + z2dot[t]^2 + s2'[t]^2)/2
Pot[t_] := m1*g*z1[t] + k1*s1[t]^2/2 + m2*g*z2[t] + k2*s2[t]^2/2

Lagrangian[t_] := Kin[t] - Pot[t]

eq0 := D[D[Lagrangian[t], \[Theta]1'[t]], t] - 
   D[Lagrangian[t], \[Theta]1[t]] == 0
eq1 := D[D[Lagrangian[t], \[Phi]1'[t]], t] - 
   D[Lagrangian[t], \[Phi]1[t]] == 0
eq2 := D[D[Lagrangian[t], s1'[t]], t] - D[Lagrangian[t], s1[t]] == 0
eq3 := D[D[Lagrangian[t], \[Theta]2'[t]], t] - 
   D[Lagrangian[t], \[Theta]2[t]] == 0
eq4 := D[D[Lagrangian[t], \[Phi]2'[t]], t] - 
   D[Lagrangian[t], \[Phi]2[t]] == 0
eq5 := D[D[Lagrangian[t], s2'[t]], t] - D[Lagrangian[t], s2[t]] == 0

sol := First[
  NDSolve[{eq0, eq1, eq2, eq3, eq4, 
    eq5, \[Theta]1[0] == Pi/3, \[Phi]1[0] == 0, 
    s1[0] == 1, \[Theta]2[0] == Pi/3, \[Phi]2[0] == 0, 
    s2[0] == 2, \[Theta]1'[0] == 0, \[Phi]1'[0] == 3, 
    s1'[0] == 0, \[Theta]2'[0] == 0, \[Phi]2'[0] == 0, 
    s2'[0] == 0}, {\[Theta]1, \[Phi]1, s1, \[Theta]2, \[Phi]2, 
    s2}, {t, 0, 100}]]

x1[t_] := 
 Evaluate[(l + s1[t]) Sin[\[Theta]1[t]] Sin[\[Phi]1[t]] /. sol]
y1[t_] := 
 Evaluate[(l + s1[t]) Sin[\[Theta]1[t]] Cos[\[Phi]1[t]] /. sol]
z1[t_] := Evaluate[-(l + s1[t]) Cos[\[Theta]1[t]] /. sol]
x2[t_] := 
 Evaluate[(l1 + s1[t]) Sin[\[Theta]1[t]] Sin[\[Phi]1[t]] + (l2 + 
       s2[t]) Sin[\[Theta]2[t]] Sin[\[Phi]2[t]] /. sol]
y2[t_] := 
 Evaluate[(l1 + s1[t]) Sin[\[Theta]1[t]] Cos[\[Phi]1[t]] + (l2 + 
       s2[t]) Sin[\[Theta]2[t]] Cos[\[Phi]2[t]] /. sol]
z2[t_] := 
 Evaluate[-(l1 + s1[t]) Cos[\[Theta]1[t]] - (l2 + 
       s2[t]) Cos[\[Theta]2[t]] /. sol]


frames = Table[
   Graphics3D[{Gray, Thick, 
     Line[{{0, 0, 0}, {x1[t], y1[t], z1[t]}, {x2[t], y2[t], z2[t]}}], 
     Darker[Blue], Sphere[{0, 0, 0}, .1], Darker[Red], 
     Sphere[{x1[t], y1[t], z1[t]}, .1], 
     Sphere[{x1[t], y1[t], z1[t]}, .2], 
     Sphere[{x2[t], y2[t], z2[t]}, .2]}, 
    PlotRange -> {{-10, 10}, {-10, 10}, {-5, 0}}], {t, 0, 20, .1}];
ListAnimate[frames]

Could someone check my code and tell me what's wrong? Also, above I added the picture with transition between the generalised coordinates and the Cartesian ones, as usual.

$\endgroup$
4

1 Answer 1

2
$\begingroup$

Parameters k1=k2=5 are too small, let take k1=k2=50 then we have

m1 = 1; m2 = 1; l1 = 1; l2 = 1; k1 = 50; k2 = 50; g = 9.81;

x1 = (l1 + s1[t]) Sin[\[Theta]1[t]] Sin[\[Phi]1[t]];
y1 = (l1 + s1[t]) Sin[\[Theta]1[t]] Cos[\[Phi]1[t]];
z1 = -(l1 + s1[t]) Cos[\[Theta]1[t]];
x2 = (l1 + s1[t]) Sin[\[Theta]1[t]] Sin[\[Phi]1[t]] + (l2 + 
      s2[t]) Sin[\[Theta]2[t]] Sin[\[Phi]2[t]];
y2 = (l1 + s1[t]) Sin[\[Theta]1[t]] Cos[\[Phi]1[t]] + (l2 + 
      s2[t]) Sin[\[Theta]2[t]] Cos[\[Phi]2[t]];
z2 = -(l1 + s1[t]) Cos[\[Theta]1[t]] - (l2 + s2[t]) Cos[\[Theta]2[t]];

x1dot = D[x1, t];
y1dot = D[y1, t];
z1dot = D[z1, t];
x2dot = D[x2, t];
y2dot = D[y2, t];
z2dot = D[z2, t];

Kin = m1 (x1dot^2 + y1dot^2 + z1dot^2)/2 + 
   m2 (x2dot^2 + y2dot^2 + z2dot^2)/2;
Pot = m1*g*z1 + k1*s1[t]^2/2 + m2*g*z2 + k2*s2[t]^2/2;

Lagrangian = Kin - Pot;

eq0 = D[D[Lagrangian, \[Theta]1'[t]], t] - 
    D[Lagrangian, \[Theta]1[t]] == 0;
eq1 = D[D[Lagrangian, \[Phi]1'[t]], t] - D[Lagrangian, \[Phi]1[t]] == 
   0;
eq2 = D[D[Lagrangian, s1'[t]], t] - D[Lagrangian, s1[t]] == 0;
eq3 = D[D[Lagrangian, \[Theta]2'[t]], t] - 
    D[Lagrangian, \[Theta]2[t]] == 0;
eq4 = D[D[Lagrangian, \[Phi]2'[t]], t] - D[Lagrangian, \[Phi]2[t]] == 
   0;
eq5 = D[D[Lagrangian, s2'[t]], t] - D[Lagrangian, s2[t]] == 0;


sol = First[
  NDSolve[{eq0, eq1, eq2, eq3, eq4, 
    eq5, \[Theta]1[0] == Pi/3, \[Phi]1[0] == 0, 
    s1[0] == 1, \[Theta]2[0] == Pi/3, \[Phi]2[0] == 0, 
    s2[0] == 2, \[Theta]1'[0] == 0, \[Phi]1'[0] == 3, 
    s1'[0] == 0, \[Theta]2'[0] == 0, \[Phi]2'[0] == 0, 
    s2'[0] == 0}, {\[Theta]1, \[Phi]1, s1, \[Theta]2, \[Phi]2, 
    s2}, {t, 0, 100}, 
   Method -> {"EquationSimplification" -> "Residual"}]]


X1[t_] := 
 Evaluate[(l1 + s1[t]) Sin[\[Theta]1[t]] Sin[\[Phi]1[t]] /. sol];
Y1[t_] := 
 Evaluate[(l1 + s1[t]) Sin[\[Theta]1[t]] Cos[\[Phi]1[t]] /. sol];
Z1[t_] := Evaluate[-(l1 + s1[t]) Cos[\[Theta]1[t]] /. sol];
X2[t_] := 
 Evaluate[(l1 + s1[t]) Sin[\[Theta]1[t]] Sin[\[Phi]1[t]] + (l2 + 
       s2[t]) Sin[\[Theta]2[t]] Sin[\[Phi]2[t]] /. sol];
Y2[t_] := 
 Evaluate[(l1 + s1[t]) Sin[\[Theta]1[t]] Cos[\[Phi]1[t]] + (l2 + 
       s2[t]) Sin[\[Theta]2[t]] Cos[\[Phi]2[t]] /. sol];
Z2[t_] := 
 Evaluate[-(l1 + s1[t]) Cos[\[Theta]1[t]] - (l2 + 
       s2[t]) Cos[\[Theta]2[t]] /. sol];

Visualization

R = Sqrt[Max[Table[X2[t]^2 + Y2[t]^2 + Z2[t]^2, {t, 0, 5, .01}]]];

frames = 
  Table[Graphics3D[{Green, Opacity[.1], Sphere[{0, 0, 0}, R], 
     Opacity[.9], Gray, Thick, 
     Line[{{0, 0, 0}, {X1[t], Y1[t], Z1[t]}, {X2[t], Y2[t], Z2[t]}}], 
     Darker[Blue], Sphere[{0, 0, 0}, .1], Darker[Red], 
     Sphere[{X1[t], Y1[t], Z1[t]}, .1], 
     Sphere[{X2[t], Y2[t], Z2[t]}, .2]}, 
    PlotRange -> {{-R, R}, {-R, R}, {-R, R}}, Boxed -> False, 
    Axes -> False, ImageSize -> 200], {t, 0, 4.9, .025}];

ListAnimate[frames]

Figure 1

$\endgroup$
4
  • $\begingroup$ Thank you very much, Alex! $\endgroup$ Nov 6, 2023 at 7:15
  • $\begingroup$ don't you know how to trace path of the movement onto Oxy, Oyz, Oxz planes? $\endgroup$ Nov 6, 2023 at 8:49
  • $\begingroup$ Use functions X1[t],Y1[t],X2[t],Y2[t] to visualize on Oxy, etc. $\endgroup$ Nov 6, 2023 at 10:23
  • $\begingroup$ Sure thing, thanks :^) $\endgroup$ Nov 6, 2023 at 11:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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