4
$\begingroup$

I would like to animate a spring pendulum. Below is my piece of code. The Lagrangian and equations are correct, in which I am sure, but the animation looks terrible and non-realistic.

Could someone help me with plotting the solution of my system of ODEs? Did I make a mistake in replacing the generalised coordinates to the Cartesian ones?

Clear[s, ϕ, t];
sol = First[
   NDSolve[{4 s''[t] + s[t]/2 + 9.81*Sin[ϕ[t]] - 
       2 (1 + s[t]) ϕ'[t]^2 == 0, 
     2 ϕ''[t] (1 + s[t])^2 + 4 ϕ'[t] (1 + s[t]) s'[t] + 
       9.81 (1 + s[t]) Cos[ϕ[t]] == 0, 
     s[0] == 1/5, ϕ[0] == 1, 
     s'[0] == 0, ϕ'[0] == 0}, {s, ϕ}, {t, 10}]];

x[t_] := Evaluate[(1 + s[t]) Cos[-ϕ[t]] /. sol]
y[t_] := Evaluate[-(1 + s[t]) Sin[-ϕ[t]] /. sol]

frames = Table[
   Graphics[{Gray, Thick, Line[{{0, 0}, {x[t], y[t]}, {x[t], y[t]}}], 
     Darker[Blue], Disk[{0, 0}, .1], Darker[Red], 
     Disk[{x[t], y[t]}, .1], Disk[{x[t], y[t]}, .1]}, 
    PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}], {t, 0, 10, .1}];
ListAnimate[frames]
Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ I am not sure that your equations are correct. Please show us the Laplacian and the constants used like mass and force constant. $\endgroup$
    – Daniel Huber
    Commented Nov 2, 2023 at 10:27
  • $\begingroup$ @DanielHuber my Lagrangian is as follows: $$m(\dot{s}^2+(l+s)^2\dot{\varphi}^2)/2+mg(l+s)\cos\varphi-ks^2/2$$ and I've taken the mass equal to 2 and $l=1, k=1/2$ $\endgroup$
    – Denis D. Bavrin
    Commented Nov 2, 2023 at 10:53
  • $\begingroup$ @DanielHuber and after re-checking I've got these ODE's: $$m\ddot{\varphi}(l+s)^2+2m\dot{\varphi}(l+s)\dot{s}+mg(l+s)\sin\varphi=0$$ $$m\ddot{s}-m(l+s)\dot{\varphi}-mg\cos\varphi+ks=0$$ But id doesn't work either $\endgroup$
    – Denis D. Bavrin
    Commented Nov 2, 2023 at 11:16

1 Answer 1

Reset to default
7
$\begingroup$

Here is the corrected code:

m = 2; l = 1; k = 1/2; g = 9.81;
la = m (s'[t]^2 + (l + s[t])^2 ph'[t]^2)/2 + 
  m g (l - s[t]) Cos[ph[t]] - k s[t]^2/2
eq1 = D[D[la, s'[t]], t] - D[la, s[t]] == 0;
eq2 = D[D[la, ph'[t]], t] - D[la, ph[t]] == 0;
sol = NDSolve[{eq1, eq2, s[0] == 1/5, ph[0] == 1, s'[0] == 0, 
    ph'[0] == 0}, {s, ph}, {t, 0, 10}][[1]]

;

The length and phi:

Plot[Evaluate[{ s[t], ph[t]} /. sol], {t, 0, 10}, PlotRange -> All]

enter image description here

x and y:

x[t_] := (1 + s[t]) Sin[ph[t]] /. sol
y[t_] := -(1 + s[t]) Cos[ph[t]] /. sol
Plot[{x[t], y[t]}, {t, 0, 10}]

enter image description here

And a manipulate to play with:

d = 1.1;
Manipulate[
 Graphics[{Line[{{0, 0}, {x[t], y[t]}}], Disk[{x[t], y[t]}, 0.05]}, 
  PlotRange -> {{-d, d}, {-d, 0}}]
 , {t, 0, 10}]

enter image description here

Addendum

There is actually still an error that I did not notice. In the potential energy m g (l - s[t]) Cos[ph[t]] the sign of s[t] should be positive. But then it turns out tat the force constant is so small, that the pendulum get stretched very much.

To get a more realistic result, where you see interaction between angular and length motions,, I increase the force constant by 50:

m = 2; l = 1; k = 50 1/2; g = 9.81; tmax = 10;
la = m (s'[t]^2 + (l + s[t])^2 ph'[t]^2)/2 + 
  m g (l + s[t]) Cos[ph[t]] - k s[t]^2/2
eq1 = D[D[la, s'[t]], t] - D[la, s[t]] == 0;
eq2 = D[D[la, ph'[t]], t] - D[la, ph[t]] == 0;
sol = NDSolve[{eq1, eq2, s[0] == 1/5, ph[0] == 1, s'[0] == 0, 
     ph'[0] == 0}, {s, ph}, {t, 0, tmax}][[1]];
Plot[Evaluate[{s[t], ph[t]} /. sol], {t, 0, tmax}, PlotRange -> All]
x[t_] := (l + s[t]) Sin[ph[t]] /. sol
y[t_] := -(l + s[t]) Cos[ph[t]] /. sol
Plot[{x[t], y[t]}, {t, 0, tmax}]
d = 3;
Manipulate[
 Graphics[{Line[{{0, 0}, {x[t], y[t]}}], Disk[{x[t], y[t]}, 0.05]}, 
  PlotRange -> {{-2, 2}, {-d, 0}}], {t, 0, tmax}]

![enter image description here

Improve this answer
$\endgroup$
8
  • $\begingroup$ thank you very much, I appreciate your labour and cooperation! $\endgroup$
    – Denis D. Bavrin
    Commented Nov 2, 2023 at 15:22
  • $\begingroup$ In this case you should upvote (clicking the upper triangle) and/or accept (clicking the checkmark) the answer. $\endgroup$
    – eldo
    Commented Nov 2, 2023 at 16:19
  • $\begingroup$ @DanialHuber Nice answer. Shouldn't the potential energy be -m g(l+s)Cos[ph]+1/2 k s^2 ? $\endgroup$
    – Ulrich Neumann
    Commented Nov 2, 2023 at 22:12
  • $\begingroup$ It depends on the direction of the y axis. $\endgroup$
    – Daniel Huber
    Commented Nov 3, 2023 at 8:03
  • $\begingroup$ @Ulrich Neumann Due to your comment I had a second look at the code and noticed an additional bug in the original code that I overlocked. Look at the addendum. $\endgroup$
    – Daniel Huber
    Commented Nov 3, 2023 at 9:10

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.