A moving reference frame with conservation of angular momentum

Hide code cell content

 
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
plt.style.use('fivethirtyeight')

A moving reference frame with conservation of angular momentum#

rolling table with a spring-mass system attached to its center

A mass is connected to a spring on top of a table. The system is at rest with the spring unstretched oriented towards the top of the table, \(\theta(t=0) = 90^o\).

Then, the table is pushed so it has a constant velocity of \(v_0\). The mass was not pushed, so its initial velocity realtive to ground is \(\bar{0}\)

\(v_{P/O} = \bar{0} = v_O\hat{i} +v_{P/O'} \rightarrow v_{P/O'} = -v_0 \hat{i}\)

The table is moving at a constant speed and does not rotate. It is an inertial reference frame. You can apply Newton’s second law to the global coordinate system or the moving table coordinate system. This notebook predicts the motion of the point \(P\) after the table is pushed for

  • \(k = 16~N/m\) spring constant

  • \(m = 0.1~kg\) mass of object at point P

  • \(L_0 = 0.1\) unstretched spring length

  • \(v_0 = 1~m/s\) speed of the moving table

  • \(r_0 = L_0\) initial position of mass is \(\mathbf{r} = r_0\hat{e}_r\) where \(\hat{e}_r = 1\hat{j}\)

  • \(\dot{\theta}_0 = \frac{v_0}{r_0}\) initial angular velocity

Conservation of angular momentum in a moving system#

The spring does not create any moment on the system, so the angular momentum is constant.

\(\sum M_{O'} = 0 = \frac{d}{dt}\bar{h}_{O'}\)

\(h_{O'} = \mathbf{r}_{P/O'}\times m \mathbf{v}_{P/O'}\)

\(h_{O'} = mr^2\dot{\theta} = constant\)

Conservation of energy in moving system#

Ignoring friction and drag, you have a system with kinetic and potential energy, \(T~and~V\), respectively. There are no nonconservative forces acting on the system.

  • \(T = \frac{1}{2}mv^2 = \frac{m}{2}(\dot{r}^2 + r^2\dot{\theta}^2)\)

  • \(V = \frac{1}{2}k(r-l_0)^2\)

\(T_1 + V_1 = T(t) + V(t)\)

\(\frac{m}{2}(\dot{r_0}^2 + r^2\dot{\theta}_0^2) + \frac{1}{2}k(r_0-l_0)^2 = \frac{m}{2}(\dot{r}^2 + r^2\dot{\theta}^2) + \frac{1}{2}k(r-l_0)^2\)

Create equation of motion and substitute \(\theta(t)\) for \(f(r)\)#

The conservation of angular momentum creates a relation between \(\dot{\theta}~and~r\) as such,

\(\dot{\theta} = \frac{r_0^2\dot{\theta}_0}{r^2}\)

You can plug this into the conservation of energy equation as such,

\(\frac{m}{2}(\dot{r_0}^2 + r^2\dot{\theta}_0^2) + \frac{1}{2}k(r_0-l_0)^2 = \frac{m}{2}\left(\dot{r}^2 + r^2\left(\frac{r_0^2\dot{\theta}_0}{r^2}\right)^2\right) + \frac{1}{2}k(r-l_0)^2\)

but, you’re still left with two unknowns on the right-hand-side \(r~and~\dot{r}\). What you have now is a relationship between distance from origin, \(r\), and its speed, \(\dot{r}\). If you want to know maximum or minimum distances, \(r_{max}~or~r_{min}\), set \(\dot{r}=0\) and solve the quadratic equation.

Here, you want to know the motion of the system, so you create the equations of motion from Newton’s second law,

  • \(\sum F_r = -k(r - l_0) = m(\ddot{r} - r\dot{\theta}^2)\)

  • \(\sum F_\theta = 0 = m(r\ddot{\theta} + 2 \dot{r}\dot{\theta})\)

Plugging in \(\dot{\theta}^2 = \left(\frac{r_0^2\dot{\theta}_0}{r^2}\right)^2\) to the \(\sum F_r\), you are left with one second order differential equation that describes \(\ddot{r} = f(r)\)

\(\ddot{r} = -\frac{k}{m}(r-l_0) + \frac{(r_0^2\dot{\theta}_0)^2}{r^3}\)

Solving for the position of the point \(P\)#

Now, you have one second order differential equation and one function for \(\dot{\theta} = f(r)\). The position of the point \(P\), \(\mathbf{r}_{P/O'} = r\hat{e}_r = r(\cos\theta \hat{i} + \sin\theta \hat{j})\). So you can solve for \(r~and~theta\) in 3 steps

  1. integrate the second order ODE to find \(r(t)\)

  2. plug in \(r(t)\) into \(\dot{\theta} = f(r)\) to find \(\dot{\theta}\)

  3. sum the \(\dot{\theta}dt\) values to find \(\theta(t)\) and animate the path of the point \(P\) as \(r\cos\theta-vs-r\sin\theta\)

1. integrate the second order ODE#

Here, you define the ODE as a function that accepts the time, t, and state, y \(=[r,~\dot{r}]\), and returns the derivatives, dy \(=\dot{y} = [\dot{r},~\ddot{r}]\).

def table_ode(t, y, dtheta0 = 10, r0 = 0.2, k = 50, m = 0.1, L0 = 0.1):
    dy = np.zeros(np.shape(y))
    dy[0] = y[1]
    dy[1] = -k/m*(y[0]-L0) + r0**4*dtheta0**2/y[0]**3
    return dy

Use the solve_ivp to integrate the ODE for 2 seconds, then plot \(r(t)\).

Hide code cell source

 
tmax = 2.5

L0 = 0.1
r0 = L0
v0 = 1
dtheta0 = v0/r0
sol = solve_ivp(lambda t, y: table_ode(t, y, r0 = r0, dtheta0 = dtheta0, L0 = L0), 
                [0, tmax], [0.2, 0], t_eval = np.linspace(0, tmax, 500))

r = sol.y[0] # radial position is r = y[0]

plt.plot(sol.t, r)
plt.title('Result from ODE integration r(t)')
plt.xlabel('time (s)')
plt.ylabel('r (m)');
../_images/114253e231c985e8aab5ccac4b32624973374c40c2cae40d63638dfbddf3940c.png

2. plug in \(r(t)\) into \(\dot{\theta} = f(r)\) to find \(\dot{\theta}\)#

The angular momentum is conserved, so if \(r\) decreases, \(\dot{\theta}\) increases. Here, you plot the angular speed of the point \(P\), \(\dot{\theta}\).

Hide code cell source

 
dtheta = dtheta0*r0**2/sol.y[0]**2
plt.plot(sol.t, dtheta*180/np.pi)
plt.title(r'$\dot{\theta} = h_0/(mr^2)$')
plt.xlabel('time (s)')
plt.ylabel(r'$\dot{\theta}$ (deg)');
../_images/d297d5bb4bd5ad6a3e80dd9c6ed9addb6865408993de02533cdb618fbc333e16.png

3. sum the \(\dot{\theta}dt\) values to find \(\theta(t)\) and animate the path of the point \(P\)#

The position is \(r_{P/O'} = r(\cos\theta \hat{i} + \sin \theta \hat{j})\). The angle, \(\theta(t) = \int_0^t \dot{\theta}dt = \sum_0^i \dot{\theta}(t_i)dt\). Here, you use the np.cumsum to calculate the integral of the angular velocity and find angle.

Hide code cell source

 
theta = np.pi/2 + np.cumsum(dtheta*sol.t[1])
plt.plot(sol.t, theta*180/np.pi)
plt.xlabel('time (s)')
plt.ylabel(r'$\theta$ (deg)');
../_images/23735155064b9fab3327f961b31fceb912729acad04fa0d7231a7187580e3451.png

Now, you have \(r(t)\) and \(\theta(t)\), so you can plot the path of the point \(P\).

Hide code cell source

 
xp = r*np.cos(theta)
yp = r*np.sin(theta)

plt.plot(xp, yp, 'o-')
ax_lims = np.array(plt.axis('equal'))
ax_lims

array([-0.21907951,  0.21764961, -0.2161024 ,  0.21979797])
../_images/8d6056cb5adbf2d1fc59c141f7b9232f4efaa25f1a0c4cbdbf56706fbe93c574.png

Now, you can animate the results and watch the the object move on the table.

Hide code cell source

 
from matplotlib import animation
from IPython.display import HTML

fig, ax = plt.subplots()

ax.axis(ax_lims*1.3)
ax.set_aspect('equal')

ax.set_xlabel('x-position (m)')
ax.set_ylabel('y-position (m)')

line, = ax.plot([], [], '--')
marker, = ax.plot([], [], 'o', markersize=10)

# 2. Create an initializing (`init`) function that clears the previous line and marker

def init():
    line.set_data([], [])
    marker.set_data([], [])
    return (line,marker,)

# 3. Create an animating (`animate`) function that updates the line

def animate(i):
    '''function that updates the line and marker data
    arguments:
    ----------
    i: index of timestep
    outputs:
    --------
    line: the line object plotted in the above ax.plot(...)
    marker: the marker for the end of the 2-bar linkage plotted above with ax.plot('...','o')'''
    line.set_data(xp[:i], yp[:i])
    marker.set_data([xp[i-1], yp[i-1]])
    return (line, marker, )

# 4. Create an animation (`anim`) variable using the `animation.FuncAnimation`

anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=range(0,len(sol.t)), interval=30, 
                               blit=True)

HTML(anim.to_html5_video())

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[7], line 43
     37 # 4. Create an animation (`anim`) variable using the `animation.FuncAnimation`
     39 anim = animation.FuncAnimation(fig, animate, init_func=init,
     40                                frames=range(0,len(sol.t)), interval=30, 
     41                                blit=True)
---> 43 HTML(anim.to_html5_video())

File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:1306, in Animation.to_html5_video(self, embed_limit)
   1302 Writer = writers[mpl.rcParams['animation.writer']]
   1303 writer = Writer(codec='h264',
   1304                 bitrate=mpl.rcParams['animation.bitrate'],
   1305                 fps=1000. / self._interval)
-> 1306 self.save(str(path), writer=writer)
   1307 # Now open and base64 encode.
   1308 vid64 = base64.encodebytes(path.read_bytes())

File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:1122, in Animation.save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs, progress_callback)
   1119 for data in zip(*[a.new_saved_frame_seq() for a in all_anim]):
   1120     for anim, d in zip(all_anim, data):
   1121         # TODO: See if turning off blit is really necessary
-> 1122         anim._draw_next_frame(d, blit=False)
   1123         if progress_callback is not None:
   1124             progress_callback(frame_number, total_frames)

File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:1157, in Animation._draw_next_frame(self, framedata, blit)
   1153 def _draw_next_frame(self, framedata, blit):
   1154     # Breaks down the drawing of the next frame into steps of pre- and
   1155     # post- draw, as well as the drawing of the frame itself.
   1156     self._pre_draw(framedata, blit)
-> 1157     self._draw_frame(framedata)
   1158     self._post_draw(framedata, blit)

File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:1789, in FuncAnimation._draw_frame(self, framedata)
   1785     self._save_seq = self._save_seq[-self._save_count:]
   1787 # Call the func with framedata and args. If blitting is desired,
   1788 # func needs to return a sequence of any artists that were modified.
-> 1789 self._drawn_artists = self._func(framedata, *self._args)
   1791 if self._blit:
   1793     err = RuntimeError('The animation function must return a sequence '
   1794                        'of Artist objects.')

Cell In[7], line 34, in animate(i)
     25 '''function that updates the line and marker data
     26 arguments:
     27 ----------
   (...)     31 line: the line object plotted in the above ax.plot(...)
     32 marker: the marker for the end of the 2-bar linkage plotted above with ax.plot('...','o')'''
     33 line.set_data(xp[:i], yp[:i])
---> 34 marker.set_data([xp[i-1], yp[i-1]])
     35 return (line, marker, )

File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/lines.py:680, in Line2D.set_data(self, *args)
    677 else:
    678     x, y = args
--> 680 self.set_xdata(x)
    681 self.set_ydata(y)

File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/lines.py:1304, in Line2D.set_xdata(self, x)
   1291 """
   1292 Set the data array for x.
   1293 
   (...)   1301 set_ydata
   1302 """
   1303 if not np.iterable(x):
-> 1304     raise RuntimeError('x must be a sequence')
   1305 self._xorig = copy.copy(x)
   1306 self._invalidx = True

RuntimeError: x must be a sequence
../_images/04e29510ab9e99f80c3cef6f311221761ff198b1487e380fb8942ef5141397ae.png

Wrapping up#

In this notebook, you used conservation of angular momentum and Newton’s second law to create an equation of motion for the radius of a spring-mass system on a moving table. Then, you plotted the results and watched the path of the object as the table slides along the floor.

Next steps:

  • What happens if you change the parameters of the system, \(k,~L0,~etc.\)?

  • What happens if you change the speed of the table?

  • How would you incorporate friction into the analysis?

Contents