# Animation of double pendulum with a spring (3DOF)

I've got a problem with animating a double pendulum where the first mass is connected with a ceiling with a bar and the second mass is connected with the first one with a spring. Below is my code which doesn't work and the reason is "Coordinate {ReplaceAll[Sin[ph[0]], NDSolve[{False, False, False, False, False, 0, False, False, 0}, {s, ph, th}, {t, 0, 100}]], ReplaceAll[-Cos[ph[0]], NDSolve[{False, False, False, False should be a pair of numbers, or a Scaled or Offset form." What did I make wrong?

m1 = 1; m2 = 1; l1 = 1; l2 = 1; k = 1/10; g = 9.81;
lagrangian :=
m1*l1^2*ph'[t]^2/2 +
m2*((l1 Cos[ph[t]] ph'[t] + Sin[th[t]] s'[t] +
Cos[th[t]] (l2 + s[t]) th'[t])^2 + (-l1 Sin[ph[t]] ph'[t] +
Sin[th[t]] s'[t] + Cos[th[t]] (l2 + s[t]) th'[t])^2 +
s'[t]^2/2)/2 + m1*9.81*l1*Cos[ph[t]] +
m2*9.81*(l1*Cos[ph[t]] + (l2 + s[t])*Sin[th[t]]) - k*s[t]^2/2;
eq1 := D[D[lagrangian, s'[t]], t] - D[lagrangian, s[t]] == 0;
eq2 := D[D[lagrangian, ph'[t]], t] - D[lagrangian, ph[t]] == 0;
eq3 := D[D[lagrangian, th'[t]], t] - D[lagrangian, th[t]] == 0;
sol := NDSolve[{eq1, eq2, eq3, s[0] == 1/2, ph[0] == Pi/2, th[0] = 0,
s'[0] == 0, ph'[0] == 0, th'[0] = 0}, {s, ph, th}, {t, 0, 100}]

x1[t_] := l1*Sin[ph[t]] /. sol
y1[t_] := -l1*Cos[ph[t]] /. sol
x2[t_] := l1*Sin[ph[t]] + (l2 + s[t])*Sin[th[t]] /. sol
y2[t_] := -l1*Cos[ph[t]] - (l2 + s[t])*Sin[th[t]] /. sol

d = 5;
Manipulate[
Graphics[{Line[{{x1[t], y1[t]}, {x2[t], y2[t]}}],
Disk[{x2[t], y2[t]}, 0.05]}, PlotRange -> {{-d, d}, {-d, d}}], {t,
0, 10}]


My Lagrangian is as follows:

$$\frac{m_1l_1^2\dot{\varphi}^2}{2}+\frac{m_2}{2}(l_1^2\dot{\varphi}^2+2l_2^2\dot{\theta}^2\cos^2\theta+2\dot{s}^2\sin^2\theta+4\dot{s}s\dot{\theta}\cos\theta+2s^2\theta^2\cos^2\theta+\dot{s}^2)+m_1gl_1\cos\varphi+m_2g(l_1\cos\varphi+(l_2+s)\sin\theta)-\frac{ks^2}{2}$$

One can check the correctness of the KE within the Lagrangian by checking the sums of (squared) tangent velocities which are:

$$\dot{x}_1=l_1\dot{\varphi}\cos\varphi,$$ $$\dot{y}_1=l_1\dot{\varphi}\sin\varphi,$$ $$\dot{x}_2=(l_1\dot{\varphi}\cos\varphi+l_2\dot{\theta}\cos\theta)+(\dot{s}\sin\theta+s\dot{\theta}\cos\theta),$$ $$\dot{y}_2=(l_1\dot{\varphi}\sin\varphi-l_2\dot{\theta}\cos\theta)-(\dot{s}\sin\theta+s\dot{\theta}\cos\theta),$$

and the normal velocity of the spring which is $$\dot{s}$$. The PE of the Lagrangian is: $$-m_1gl_1\cos\varphi-m_2g(l_1\cos\varphi+(l_2+s)\sin\theta)+\frac{ks^2}{2}.$$

• Check your constraints in NDSolve , they must be defined with "==" not "=" Nov 3, 2023 at 8:12
• Running your code one command at a time leads to the following. When trying to use sol (causing it to be evaluated) I get a bunch of division by zero errors. Nov 3, 2023 at 11:29
• @UlrichNeumann I tried putting extra = to my code, however, it didn't help. Nov 3, 2023 at 11:35
• @JyrkiLahtonen do you mean a bunch of False things in NDSolve? Nov 3, 2023 at 11:37
• Further testing (with a narrower interval for $t$, and modified initial values) revealed the following problem: x1[t] is not a number, rather it is a list of length one. I haven't used NDSolve but may be it is like all the other variants of solving equations, and outputs a list of solutions. If I edit your commands to read like x1[t_] := l1*Sin[ph[t]] /. sol[[1]] at least that problem goes away. Nov 3, 2023 at 11:50

Here my solution

parameters (changed k !)

m1 = 1; m2 = 1; l1 = 1; l2 = 1; k =100(* 1/10*); g = 9.81;


geometry

{e1, e2} = IdentityMatrix[2];
e[phi_] := e1 Sin[phi] + e2 Cos[phi]
p1 = l1 e[\[CurlyPhi]1[t]];  (* position m1 *)
p2 = l1 e[\[CurlyPhi]1[t]] + (l2 + s[t]) e[\[CurlyPhi]2[t]]; (*position m2 *)


Energies

T = 1/2 m1 D[p1, t] . D[p1, t] + 1/2 m2 D[p2, t] . D[p2, t] //Expand // FullSimplify
U = m1 g p1 . (-e2) + m2 g p2 . (-e2) + 1/2 k s[t]^2


Lagrange

lagrangian = T - U;
var = {s[t], \[CurlyPhi]1[t], \[CurlyPhi]2[t]};
eqn = D[D[lagrangian, {D[var, t]}] , t] - D[lagrangian, {var  }] ;


solution (changed s[0]==m2 g/k end simualtion time )

sol = NDSolve[{Map[# == 0 &, eqn],
s[0] ==  m2 g/k(*1/2*), \[CurlyPhi]1[0] == Pi/2, \[CurlyPhi]2[0] ==0,s'[0] == 0, \[CurlyPhi]1'[0] == 0, \[CurlyPhi]2'[0] ==0}
,{s, \[CurlyPhi]1, \[CurlyPhi]2}, {t, 0, 10}][[1]];

Plot[Evaluate[{s[t], \[CurlyPhi]1[t], \[CurlyPhi]2[t]} /. sol] , {t, 0, 10}]


animation

P1 = Function[t, Evaluate[DiagonalMatrix[{1, -1}] . (p1 /. sol)]] ;
P2 = Function[t, Evaluate[DiagonalMatrix[{1, -1}] . (p2 /. sol)]] ;
Manipulate[
Graphics[{Line[{{0, 0}, P1[t], P1[t] + P2[t]} ],Point[{{0, 0}, P1[t], P1[t] + P2[t]}]} ,
PlotRange ->{3 {-1, 1}, 4 {-1, 0}}], {t, 0, 10}]


Hope it helps!