Animation of swinging pendulum on spring

: home
: PDF (letter size)
: PDF (legal size)
: Mathematica notebook

Animation of swinging pendulum on spring

February 22, 2025 Compiled on February 22, 2025 at 5:59pm

Let original length of pendulum be $L$ and extension at instance shown be $x$ . The position vector of bob is

\bar{r} = (L + x) {\bar{b}}_{1}

Hence the velocity vector is

\begin{aligned} {(\frac{d}{d t} \bar{r})}_{A} & = {(\frac{d}{d t} \bar{r})}_{B} + \bar{ω} \times \bar{r} \\ {\bar{v}}_{A} & = {(\frac{d}{d t} (L + x) {\bar{b}}_{1})}_{B} + \bar{ω} \times ((L + x) {\bar{b}}_{1}) \\ = \dot{x} {\bar{b}}_{1} + \dot{θ} {\bar{b}}_{3} \times ((L + x) {\bar{b}}_{1}) \\ = \dot{x} {\bar{b}}_{1} + \dot{θ} (L + x) {\bar{b}}_{2} \end{aligned}

And the acceleration is

{(\frac{d}{d t} \bar{v})}_{A} = {(\frac{d}{d t} \bar{v})}_{B} + \bar{ω} \times \bar{v}

But

\begin{aligned} {(\frac{d}{d t} \bar{v})}_{B} & = \frac{d}{d t} {(\dot{x} {\bar{b}}_{1} + \dot{θ} (L + x) {\bar{b}}_{2})}_{B} \\ (2) & = \ddot{x} {\bar{b}}_{1} + (\ddot{θ} (L + x) + \dot{θ} \dot{x}) {\bar{b}}_{2} \end{aligned}

And

\begin{aligned} \bar{ω} \times \bar{v} & = \bar{ω} \times (\dot{x} {\bar{b}}_{1} + \dot{θ} (L + x) {\bar{b}}_{2}) \\ = \dot{θ} {\bar{b}}_{3} \times (\dot{x} {\bar{b}}_{1} + \dot{θ} (L + x) {\bar{b}}_{2}) \\ (3) & = \dot{θ} \dot{x} {\bar{b}}_{2} + \dot{θ} \dot{x} {\bar{b}}_{2} - {\dot{θ}}^{2} (L + x) {\bar{b}}_{1} \end{aligned}

Substituting (2,3) into (1) gives the acceleration of the bob in frame $A$ as

\begin{aligned} {\bar{a}}_{A} & = \ddot{x} {\bar{b}}_{1} + \ddot{θ} (L + x) {\bar{b}}_{2} + \dot{θ} \dot{x} {\bar{b}}_{2} + \dot{θ} \dot{x} {\bar{b}}_{2} - {\dot{θ}}^{2} (L + x) {\bar{b}}_{1} \\ = \ddot{x} {\bar{b}}_{1} + \ddot{θ} (L + x) {\bar{b}}_{2} + 2 \dot{θ} \dot{x} {\bar{b}}_{2} - {\dot{θ}}^{2} (L + x) {\bar{b}}_{1} \end{aligned}

To obtain $\bar{F} = m \bar{a}$ we now just need to ﬁnd $\bar{F}$ . The force on bob is given by just the weight and the spring force acting on it. The spring force is proportional to spring coeﬃcient $k$ times the extension. $F_{s} = - k x$ . Hence the force vector is

\bar{F} = m g \cos θ {\bar{b}}_{1} - m g \sin θ {\bar{b}}_{2} - k x {\bar{b}}_{1}

Therefore the equation of motion is

\begin{aligned} \bar{F} & = m \bar{a} \\ m g \cos θ {\bar{b}}_{1} - m g \sin θ {\bar{b}}_{2} - k x {\bar{b}}_{1} & = m (\ddot{x} {\bar{b}}_{1} + \ddot{θ} (L + x) {\bar{b}}_{2} + 2 \dot{θ} \dot{x} {\bar{b}}_{2} - {\dot{θ}}^{2} (L + x) {\bar{b}}_{1}) \\ g \cos θ {\bar{b}}_{1} - g \sin θ {\bar{b}}_{2} - \frac{k}{m} x {\bar{b}}_{1} & = \ddot{x} {\bar{b}}_{1} + \ddot{θ} (L + x) {\bar{b}}_{2} + 2 \dot{θ} \dot{x} {\bar{b}}_{2} - {\dot{θ}}^{2} (L + x) {\bar{b}}_{1} \end{aligned}

Or, by equating each vector component

\begin{aligned} \ddot{x} - {\dot{θ}}^{2} (L + x) + \frac{k}{m} x & = g \cos θ \\ \ddot{θ} (L + x) + 2 \dot{θ} \dot{x} & = - g \sin θ \end{aligned}

\begin{aligned} \ddot{x} & = {\dot{θ}}^{2} (L + x) - \frac{k}{m} x + g \cos θ \\ \ddot{θ} & = - \frac{2 \dot{θ} \dot{x}}{L + x} - \frac{g}{L + x} \sin θ \end{aligned}

Let initial conditions be $L = 1, x (0) = 0.1, \dot{x} (0) = 0, θ (0) = 20^{0}, \dot{θ} (0) = .1$ (rad/sec).

We can now solve for $x, θ$ as function of time and make the animations. We can assume values for $m, k$ and $g = 9.81$ .

The solution will give $x (t), θ (t)$ . To plot the path, we need to express $x (t)$ in frame $A$ coordinates of course which is the ﬁxed frame. But this is easy, since the pendulum frame $B$ is ﬁxed at the base on the frame A origin.

The following is plot of the solution using $k = 0.99, m = 0.09$

The following is a quick animation of the above

The following is the code used

(*Nasser M. Abbasi, Jan 25, 2025*) 
 
SetDirectory[NotebookDirectory[]] 
 
L = 1; 
k = .99; 
m = .09; 
g = 9.81; 
ode1 = x''[t] == z'[t]^2*(L + x[t]) + g*Cos[z[t]] - k/m*x[t]; 
ode2 = z''[t] == -2 z'[t]*x'[t]/(L + x[t]) - g*Sin[z[t]]/(L + x[t]); 
IC = {x[0] == .1, x'[0] == 0, z[0] == 20 Degree, z'[0] == .1}; 
sol = NDSolve[{ode1, ode2, IC}, {x, z}, {t, 0, 100}] 
 
Manipulate[ 
 Module[{currentX, currentY}, 
  currentX = (L + Evaluate[x[t0] /. sol])*Sin[Evaluate[z[t0] /. sol]]; 
  currentY = (L + Evaluate[x[t0] /. sol])*Cos[Evaluate[z[t0] /. sol]]; 
  Grid[{ 
    {Graphics[{ 
       Line[{{0, 0}, {First@currentX, -First@currentY}}], 
       {Blue, Disk[{First@currentX, -First@currentY}, .1]} 
       }, 
      PlotRange -> {{-2, 2}, {-4, 1}}, GridLines -> Automatic, 
      GridLinesStyle -> LightGray 
      ]}}] 
  ], 
 {{t0, 0, "time"}, 0, 10, .01}, 
 TrackedSymbols :> {t0} 
 ]