1
$\begingroup$

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}.$$

Thanks in advance.the graphic

$\endgroup$
16
  • $\begingroup$ Check your constraints in NDSolve , they must be defined with "==" not "=" $\endgroup$ Nov 3, 2023 at 8:12
  • $\begingroup$ 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. $\endgroup$ Nov 3, 2023 at 11:29
  • $\begingroup$ @UlrichNeumann I tried putting extra = to my code, however, it didn't help. $\endgroup$ Nov 3, 2023 at 11:35
  • $\begingroup$ @JyrkiLahtonen do you mean a bunch of False things in NDSolve? $\endgroup$ Nov 3, 2023 at 11:37
  • $\begingroup$ 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. $\endgroup$ Nov 3, 2023 at 11:50

1 Answer 1

2
$\begingroup$

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}] 

enter image description here

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}]

enter image description here

Hope it helps!

$\endgroup$

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.