Animating a double spring spherical pendulum

I 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.

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]


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