20
$\begingroup$

What would be a theoretical explanation of an "ideal" 14:1 mass ratio in this experiment, also demonstrated in this video?

The experiment ties one nut to one end of the string and 14 nuts to the other, then holds the string like this and lets go:

enter image description here

The end with the single nut ends up wrapped round your finger and stops the nuts falling to the floor:

enter image description here

Why is a 14:1 mass ratio required for this to happen?

EDIT:

Here is the set of equations I'm trying myself for this problem:

$ l(t) = r^2(t) \alpha'(t) \\ T(t) = \text{max}(\mu g - k \alpha(t), 0) \\ l'(t) = g \cos{\alpha(t)} r(t) + T(t) r_0 \\ r''(t) = g \sin{\alpha(t)} - T(t) $

With $l(t)$ - angular momentum divided by smaller mass, $T(t)$ - string tension, $\mu$ - mass proportion, $k$ - friction coefficient, $r_0$ - pivot radius, $r(t)$ - string length, $\alpha(t)$ - string angle.

$\endgroup$
2
  • $\begingroup$ I can't justify it, but I bet it would also work with different (non-multiple of 14:1) ratios. I'm almost stopping working to try it myself... $\endgroup$ Mar 23, 2020 at 13:13
  • $\begingroup$ I tried to solve it analytically using Lagrangian mechanics but I can't solve the equations of motions. Does anyone know if a solution exist? $\endgroup$
    – PinkFloyd
    Oct 25, 2023 at 15:30

1 Answer 1

25
$\begingroup$

TL;DR: Mass ratio = 14 is not particularly special, but it is in a special region of mass ratios (about 11 to 14) that has optimal properties to wind the rope around the finger as much as possible.

If you want to know why read the answer. If you just want to look at pretty gifs check it out (hat tip to @Ruslan for the animation idea).

ed

One can actually learn a lot from these movies! Especially if one considers that friction kicks in after probably about 2 windings to stop the rope from slipping along the finger, one can identify which mass ratios should work in practice. Only the experiment can tell the full result since there are a lot more factors not considered in the model here (air resistance, non-ideal rope, finite finger thickness, finger movement...).

Code for animations if you want to run it yourself or adapt the equations of motion to someting fancy (such as including friction):

import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import cm

import numpy as np

# integrator for ordinary differential equations
from scipy.integrate import ode

def eoms_pendulum(t, y, params):
    """
    Equations of motion for the simple model.
    I was too dumb to do the geometry elegantly, so there are case distinctions...
    """
    # unpack #
    v1_x, v1_y, v2, x1, y1, y2 = y
    m1, m2, g, truncate_at_inversion = params
    if x1<=0 and y1<=0:
        # calc helpers #
        F1_g = m1*g
        F2_g = m2*g                        # _g for "gravity"
        L_swing = np.sqrt( x1**2 + y1**2 ) # distance of mass 1 to the pendulum pivot
        Theta = np.arctan(y1/x1)               # angle
        dt_Theta = ( v1_y/x1 - v1_x*y1/x1**2)/(1. + y1**2/x1**2) # derivative of arctan
        help_term = -F2_g/m2 - F1_g/m1 * np.sin(Theta) - v1_x*np.sin(Theta)*dt_Theta + v1_y*np.cos(Theta)*dt_Theta
        F_r = help_term / ( -1./m1 - 1./m2 ) # _r for "rope", this formula comes from requiring a constant length rope
        # calc derivatives
        dt_v1_x = ( F_r*np.cos(Theta) ) / m1
        dt_v1_y = ( -F1_g + F_r*np.sin(Theta) ) / m1
        dt_v2   = ( F_r - F2_g ) / m2
        dt_x1   = v1_x
        dt_y1   = v1_y
        dt_y2   = v2
    elif x1>=0 and y1<=0:
        # calc helpers #
        F1_g = m1*g
        F2_g = m2*g
        L_swing = np.sqrt( x1**2 + y1**2 )
        Theta = np.arctan(-x1/y1)
        dt_Theta = -( v1_x/y1 - v1_y*x1/y1**2)/(1. + x1**2/y1**2)
        help_term = -F2_g/m2 - F1_g/m1 * np.cos(Theta) - v1_x*np.cos(Theta)*dt_Theta - v1_y*np.sin(Theta)*dt_Theta
        F_r = help_term / ( -1./m1 - 1./m2 )
        # calc derivatives
        dt_v1_x = ( -F_r*np.sin(Theta) ) / m1
        dt_v1_y = ( -F1_g + F_r*np.cos(Theta) ) / m1
        dt_v2   = ( F_r - F2_g ) / m2
        dt_x1   = v1_x
        dt_y1   = v1_y
        dt_y2   = v2
    elif x1>=0 and y1>=0:
        # calc helpers #
        F1_g = m1*g
        F2_g = m2*g
        L_swing = np.sqrt( x1**2 + y1**2 )
        Theta = np.arctan(y1/x1)
        dt_Theta = ( v1_y/x1 - v1_x*y1/x1**2)/(1. + y1**2/x1**2)
        help_term = -F2_g/m2 + F1_g/m1 * np.sin(Theta) + v1_x*np.sin(Theta)*dt_Theta - v1_y*np.cos(Theta)*dt_Theta
        F_r = help_term / ( -1./m1 - 1./m2 )
        # calc derivatives
        dt_v1_x = ( -F_r*np.cos(Theta) ) / m1
        dt_v1_y = ( -F1_g - F_r*np.sin(Theta) ) / m1
        dt_v2   = ( F_r - F2_g ) / m2
        dt_x1   = v1_x
        dt_y1   = v1_y
        dt_y2   = v2
    elif x1<=0 and y1>=0:
        # calc helpers #
        F1_g = m1*g
        F2_g = m2*g
        L_swing = np.sqrt( x1**2 + y1**2 )
        Theta = np.arctan(-y1/x1)
        dt_Theta = -( v1_y/x1 - v1_x*y1/x1**2)/(1. + y1**2/x1**2)
        help_term = -F2_g/m2 + F1_g/m1 * np.sin(Theta) - v1_x*np.sin(Theta)*dt_Theta - v1_y*np.cos(Theta)*dt_Theta
        F_r = help_term / ( -1./m1 - 1./m2 )
        # calc derivatives
        dt_v1_x = ( F_r*np.cos(Theta) ) / m1
        dt_v1_y = ( -F1_g - F_r*np.sin(Theta) ) / m1
        dt_v2   = ( F_r - F2_g ) / m2
        dt_x1   = v1_x
        dt_y1   = v1_y
        dt_y2   = v2
    if truncate_at_inversion:
        if dt_y2 > 0.:
            return np.zeros_like(y)
    return [dt_v1_x, dt_v1_y, dt_v2, dt_x1, dt_y1, dt_y2]

def total_winding_angle(times, trajectory):
    """
    Calculates the total winding angle for a given trajectory
    """
    dt = times[1] - times[0]
    v1_x, v1_y, v2, x1, y1, y2 = [trajectory[:, i] for i in range(6)]
    dt_theta = ( x1*v1_y - y1*v1_x ) / np.sqrt(x1**2 + y1**2) # from cross-product
    theta_tot = np.cumsum(dt_theta) * dt
    return theta_tot

################################################################################
### setup ###
################################################################################
trajectories = []
m1 = 1
m2_list = np.arange(2, 20, 2)[0:9]
ntimes = 150
for m2 in m2_list:
    # params #       
    params = [
              m1,    # m1
              m2,  # m2
              9.81, # g
              False # If true, truncates the motion when m2 moves back upwards
             ]

    # initial conditions #
    Lrope = 1.0 # Length of the rope, initially positioned such that m1 is L from the pivot
    init_cond = [
                    0.0, # v1_x
                    0., # v1_y
                    0., # v2
                    -Lrope/2, # x1
                    0.0, # y1
                    -Lrope/2, # y2
                ]

    # integration time range #
    times = np.linspace(0, 1.0, ntimes)

    # trajectory array to store result #
    trajectory = np.empty((len(times), len(init_cond)), dtype=np.float64)

    # helper #
    show_prog = True

    # check eoms at starting position #
    #print(eoms_pendulum(0, init_cond, params))

    ################################################################################
    ### numerical integration ###
    ################################################################################

    r = ode(eoms_pendulum).set_integrator('zvode', method='adams',
                                               with_jacobian=False) # integrator and eoms
    r.set_initial_value(init_cond, times[0]).set_f_params(params)   # setup
    dt = times[1] - times[0] # time step

    # integration (loop time step) 
    for i, t_i in enumerate(times):
        trajectory[i,:] = r.integrate(r.t+dt) # integration

    trajectories.append(trajectory)

    # ### extract ###
    # x1 = trajectory[:, 3]
    # y1 = trajectory[:, 4]
    # x2 = np.zeros_like(trajectory[:, 5])
    # y2 = trajectory[:, 5]

    # L = np.sqrt(x1**2 + y1**2) # rope part connecting m1 and pivot
    # Ltot = -y2 + L             # total rope length


################################################################################
### Visualize trajectory ###
################################################################################


import numpy as np
from matplotlib import pyplot as plt
from matplotlib.animation import FuncAnimation
plt.style.use('seaborn-pastel')

n=3
m=3

axes = []
m1_ropes = []
m2_ropes = []
m1_markers = []
m2_markers = []

fig = plt.figure(figsize=(10,10))
for sp, m2_ in enumerate(m2_list):
    ax = fig.add_subplot(n, m, sp+1, xlim=(-0.75, 0.75), ylim=(-1, 0.5), xticks=[], yticks=[])
    m1_rope, = ax.plot([], [], lw=1, color='k')
    m2_rope, = ax.plot([], [], lw=1, color='k')
    m1_marker, = ax.plot([], [], marker='o', markersize=10, color='r', label=r'$m_1 = {}$'.format(m1))
    m2_marker, = ax.plot([], [], marker='o', markersize=10, color='b', label=r'$m_2 = {}$'.format(m2_))
    axes.append(ax)
    m1_ropes.append(m1_rope)
    m2_ropes.append(m2_rope)
    m1_markers.append(m1_marker)
    m2_markers.append(m2_marker)
    ax.set_aspect('equal', adjustable='box')
    ax.legend(loc='upper left', fontsize=12, ncol=2, handlelength=1, bbox_to_anchor=(0.1, 1.06))
    plt.tight_layout()

def init():
    for m1_rope, m2_rope, m1_marker, m2_marker in zip(m1_ropes, m2_ropes, m1_markers, m2_markers):
        m1_rope.set_data([], [])
        m2_rope.set_data([], [])
        m1_marker.set_data([], [])
        m2_marker.set_data([], [])
    return (*m1_ropes, *m2_ropes, *m1_markers, *m2_markers)
def animate(i):
    for sp, (m1_rope, m2_rope, m1_marker, m2_marker) in enumerate(zip(m1_ropes, m2_ropes, m1_markers, m2_markers)):
        x1 = trajectories[sp][:, 3]
        y1 = trajectories[sp][:, 4]
        x2 = np.zeros_like(trajectories[sp][:, 5])
        y2 = trajectories[sp][:, 5]
        m1_rope.set_data([x1[i], 0], [y1[i], 0])
        m2_rope.set_data([x2[i], 0], [y2[i], 0])
        m1_marker.set_data(x1[i], y1[i])
        m2_marker.set_data(x2[i], y2[i])
    return (*m1_ropes, *m2_ropes, *m1_markers, *m2_markers)

anim = FuncAnimation(fig, animate, init_func=init,
                               frames=len(trajectories[0][:, 0]), interval=500/ntimes, blit=True)


anim.save('PendulumAnim.gif', writer='imagemagick', dpi = 50)

plt.show()

Main argument

Winding angle behavior in the no friction, thin pivot case

My answer is based on a simple model for the system(no fricition, infinitely thin pivot, ideal rope, see also detailed description below), from which one can actually get some very nice insight on why the region around 14 is special.

As a quantity of interest, we define a winding angle as a function of time $\theta(t)$. It indicates, which total angle the small mass has travelled around the finger. $\theta(t)=2\pi$ corresponds to one full revolution, $\theta(t)=4\pi$ corresponds to two revolutions and so on.

One can then plot the winding angle as a function of time and mass ratio for the simple model:

enter image description here

The color axis shows the winding angle. We can clearly see that between mass ratio 12-14, the winding angle goes up continuously in time and reaches a high maximum. The first maxima in time for each mass ratio are indicated by the magenta crosses. Also note that the weird discontinuities are places where the swinging mass goes through zero/hits the finger, where the winding angle is not well defined.

To see the behaviour in a bit more detail, let us look at some slices of the 2D plot (2$\pi$ steps/full revolutions marked as horizontal lines):

enter image description here

We see that mass ratios 12, 13, 14 behave very similarly. 16 has a turning point after 4 revolutions, but I would expect this to still work in practice, since when the rope is wrapped 4 times around the finger, there should be enough friction to clip it.

For mass ratio 5, on the other hand, we do not even get 2 revolutions and the rope would probably slip.

If you want to reproduce these plots, here is my code. Feel free to make adaptions and post them as an answer. It would be interesting, for example, if one can include friction in a simple way to quantify the clipping effect at the end. I imagine this will be hard though, and one would need at least one extra parameter.

import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib import cm

import numpy as np

from scipy.signal import argrelextrema

# integrator for ordinary differential equations
from scipy.integrate import ode

def eoms_pendulum(t, y, params):
    """
    Equations of motion for the simple model.
    I was too dumb to do the geometry elegantly, so there are case distinctions...
    """
    # unpack #
    v1_x, v1_y, v2, x1, y1, y2 = y
    m1, m2, g, truncate_at_inversion = params
    if x1<=0 and y1<=0:
        # calc helpers #
        F1_g = m1*g
        F2_g = m2*g                        # _g for "gravity"
        L_swing = np.sqrt( x1**2 + y1**2 ) # distance of mass 1 to the pendulum pivot
        Theta = np.arctan(y1/x1)               # angle
        dt_Theta = ( v1_y/x1 - v1_x*y1/x1**2)/(1. + y1**2/x1**2) # derivative of arctan
        help_term = -F2_g/m2 - F1_g/m1 * np.sin(Theta) - v1_x*np.sin(Theta)*dt_Theta + v1_y*np.cos(Theta)*dt_Theta
        F_r = help_term / ( -1./m1 - 1./m2 ) # _r for "rope", this formula comes from requiring a constant length rope
        # calc derivatives
        dt_v1_x = ( F_r*np.cos(Theta) ) / m1
        dt_v1_y = ( -F1_g + F_r*np.sin(Theta) ) / m1
        dt_v2   = ( F_r - F2_g ) / m2
        dt_x1   = v1_x
        dt_y1   = v1_y
        dt_y2   = v2
    elif x1>=0 and y1<=0:
        # calc helpers #
        F1_g = m1*g
        F2_g = m2*g
        L_swing = np.sqrt( x1**2 + y1**2 )
        Theta = np.arctan(-x1/y1)
        dt_Theta = -( v1_x/y1 - v1_y*x1/y1**2)/(1. + x1**2/y1**2)
        help_term = -F2_g/m2 - F1_g/m1 * np.cos(Theta) - v1_x*np.cos(Theta)*dt_Theta - v1_y*np.sin(Theta)*dt_Theta
        F_r = help_term / ( -1./m1 - 1./m2 )
        # calc derivatives
        dt_v1_x = ( -F_r*np.sin(Theta) ) / m1
        dt_v1_y = ( -F1_g + F_r*np.cos(Theta) ) / m1
        dt_v2   = ( F_r - F2_g ) / m2
        dt_x1   = v1_x
        dt_y1   = v1_y
        dt_y2   = v2
    elif x1>=0 and y1>=0:
        # calc helpers #
        F1_g = m1*g
        F2_g = m2*g
        L_swing = np.sqrt( x1**2 + y1**2 )
        Theta = np.arctan(y1/x1)
        dt_Theta = ( v1_y/x1 - v1_x*y1/x1**2)/(1. + y1**2/x1**2)
        help_term = -F2_g/m2 + F1_g/m1 * np.sin(Theta) + v1_x*np.sin(Theta)*dt_Theta - v1_y*np.cos(Theta)*dt_Theta
        F_r = help_term / ( -1./m1 - 1./m2 )
        # calc derivatives
        dt_v1_x = ( -F_r*np.cos(Theta) ) / m1
        dt_v1_y = ( -F1_g - F_r*np.sin(Theta) ) / m1
        dt_v2   = ( F_r - F2_g ) / m2
        dt_x1   = v1_x
        dt_y1   = v1_y
        dt_y2   = v2
    elif x1<=0 and y1>=0:
        # calc helpers #
        F1_g = m1*g
        F2_g = m2*g
        L_swing = np.sqrt( x1**2 + y1**2 )
        Theta = np.arctan(-y1/x1)
        dt_Theta = -( v1_y/x1 - v1_x*y1/x1**2)/(1. + y1**2/x1**2)
        help_term = -F2_g/m2 + F1_g/m1 * np.sin(Theta) - v1_x*np.sin(Theta)*dt_Theta - v1_y*np.cos(Theta)*dt_Theta
        F_r = help_term / ( -1./m1 - 1./m2 )
        # calc derivatives
        dt_v1_x = ( F_r*np.cos(Theta) ) / m1
        dt_v1_y = ( -F1_g - F_r*np.sin(Theta) ) / m1
        dt_v2   = ( F_r - F2_g ) / m2
        dt_x1   = v1_x
        dt_y1   = v1_y
        dt_y2   = v2
    if truncate_at_inversion:
        if dt_y2 > 0.:
            return np.zeros_like(y)
    return [dt_v1_x, dt_v1_y, dt_v2, dt_x1, dt_y1, dt_y2]

def total_winding_angle(times, trajectory):
    """
    Calculates the total winding angle for a given trajectory
    """
    dt = times[1] - times[0]
    v1_x, v1_y, v2, x1, y1, y2 = [trajectory[:, i] for i in range(6)]
    dt_theta = ( x1*v1_y - y1*v1_x ) / (x1**2 + y1**2) # from cross-product
    theta_tot = np.cumsum(dt_theta) * dt
    return theta_tot

def find_nearest_idx(array, value):
    """
    Find the closest element in an array and return the corresponding index.
    """
    array = np.asarray(array)
    idx = (np.abs(array-value)).argmin()
    return idx

################################################################################
### setup ###
################################################################################
theta_tot_traj_list = []

# scan mass ratio
m2_list = np.linspace(5,17,200)
for m2_ in m2_list:
    # params #       
    params = [
              1,    # m1
              m2_,  # m2
              9.81, # g
              False # If true, truncates the motion when m2 moves back upwards
             ]

    # initial conditions #
    Lrope = 1.0 # Length of the rope, initially positioned such that m1 is L from the pivot
    init_cond = [
                    0.0, # v1_x
                    0., # v1_y
                    0., # v2
                    -Lrope/2, # x1
                    0.0, # y1
                    -Lrope/2, # y2
                ]

    # integration time range #
    times = np.linspace(0, 2.2, 400)

    # trajectory array to store result #
    trajectory = np.empty((len(times), len(init_cond)), dtype=np.float64)

    # helper #
    show_prog = True

    # check eoms at starting position #
    #print(eoms_pendulum(0, init_cond, params))

    ################################################################################
    ### numerical integration ###
    ################################################################################

    r = ode(eoms_pendulum).set_integrator('zvode', method='adams',
                                               with_jacobian=False) # integrator and eoms
    r.set_initial_value(init_cond, times[0]).set_f_params(params)   # setup
    dt = times[1] - times[0] # time step

    # integration (loop time step) 
    for i, t_i in enumerate(times):
        trajectory[i,:] = r.integrate(r.t+dt) # integration

    ### extract ###
    x1 = trajectory[:, 3]
    y1 = trajectory[:, 4]
    x2 = np.zeros_like(trajectory[:, 5])
    y2 = trajectory[:, 5]

    L = np.sqrt(x1**2 + y1**2) # rope part connecting m1 and pivot
    Ltot = -y2 + L             # total rope length

    theta_tot_traj = total_winding_angle(times, trajectory)

    theta_tot_traj_list.append(theta_tot_traj)

theta_tot_traj_list = np.asarray(theta_tot_traj_list)

#maxima_idxs = np.argmax(theta_tot_traj_list, axis=-1)
maxima_idxs = []
for i,m_ in enumerate(m2_list):
    maxima_idx = argrelextrema(theta_tot_traj_list[i,:], np.greater)[0]
    if maxima_idx.size == 0:
        maxima_idxs.append(-1)
    else:
        maxima_idxs.append(maxima_idx[0])
maxima_idxs = np.asarray(maxima_idxs)

### 2D plot ###
fig=plt.figure()

plt.axhline(14, color='r', linewidth=2, dashes=[1,1])
plt.imshow(theta_tot_traj_list, aspect='auto', origin='lower',
            extent = [times[0], times[-1], m2_list[0], m2_list[-1]])
plt.plot(times[maxima_idxs], m2_list,'mx')
plt.xlabel("Time")
plt.ylabel("Mass ratio")
plt.title("Winding angle")
plt.colorbar()

fig.savefig('winding_angle.png')

plt.show()

fig=plt.figure()
slice_list = [5, 12, 13, 14, 16]
for x_ in [0.,1.,2.,3.,4.,5.]:
    plt.axhline(x_*2.*np.pi, color='k', linewidth=1)
for i, slice_val in enumerate(slice_list):
    slice_idx = find_nearest_idx(m2_list, slice_val)
    plt.plot(times, theta_tot_traj_list[slice_idx, :], label='Mass ratio: {}'.format(slice_val))

plt.xlabel('Time')
plt.ylabel('Winding angle')
plt.legend()

fig.savefig('winding_angle2.png')
plt.show()

Details

The simple model

The simple model used above and (probably) the simplest way to model the system is to assume:

  • An ideal infinitely thin rope.
  • An infinitely thin pivot than the rope wraps around (the finger in the video).
  • No friction.

Especially the no friction assumption is clearly flawed, because the effect of stopping completely relies on friction. But as we saw above one can still get some insight into the initial dynamics anyway and then think about what friction will do to change this. If someone feels motivated I challenge you to include friction into the model and change my code!

Under these assumptions, one can set up a system of coupled differential equations using Newton's laws, which can easily be solved numerically. I won't go into detail on the geometry and derivation, I'll just give some code below for people to check and play with. Disclaimer: I am not sure my equations of motion are completely right. I did some checks and it looks reasonable, but feel free to fill in your own version and post an answer.

Geometry

The geometry assumed is like this:

enter image description here

From the picture, we can get the equations of motion as follows: $$ m_1 \dot{v}_{x,1} = F_\mathrm{rope} \cos(\theta) \,, \\ m_1 \dot{v}_{y,1} = -F_{g,1} + F_\mathrm{rope} \sin(\theta) \,, \\ m_2 \dot{v}_{y,2} = -F_{g,2} + F_\mathrm{rope} \,, \\ \dot{x}_1 = v_{x,1} \,, \\ \dot{y}_1 = v_{y,1} \,, \\ \dot{y}_2 = v_{y,2} \,. $$

This is just Newton's laws for the geometry wirtten as a set of first order coupled differential equations, which can easily be solved in scipy (see code).

The hard bit is to find the rope force $F_\textrm{rope}$. It is constraint by the ideal rope condition, that the total rope length does not change in time. Following this through I got

$$ F_\textrm{rope} = \frac{\frac{F_{g,2}}{m_2} + \frac{F_{g,1}}{m_1}\sin(\theta) + v_{x,1}\sin(\theta)\dot{\theta} - v_{y,1}\cos(\theta)\dot{\theta}}{\frac{1}{m_1} + \frac{1}{m_2}} \,. $$

Note that my way of writing the solution is not particularly elegant and as a result some of these formulas only apply in the lower left quadrant ($x_1<0$, $y_1<0$). The other quadrants are implemented in the code too.

As the initial position, we will consider $x_1 = -L/2$, $y_1 = -L/2$, similarly to the video. $y_1$ does not matter too much, it simply causes an overall displacement of mass 2. We set $L=1$ and $g=9.81$. Someone else can work out the units ;-)

Let's do it in python

I already gave some code snippets above. You need numpy and matplotlib to run it. Maybe python3 would be good. If you want to plot static trajectories you can used:

################################################################################
### setup ###
################################################################################

# params #       
params = [
          1,    # m1
          14.,  # m2
          9.81, # g
          False # If true, truncates the motion when m2 moves back upwards
         ]

# initial conditions #
Lrope = 1.0 # Length of the rope, initially positioned such that m1 is L from the pivot
init_cond = [
                0.0, # v1_x
                0., # v1_y
                0., # v2
                -Lrope/2, # x1
                0.0, # y1
                -Lrope/2, # y2
            ]

# integration time range #
times = np.linspace(0, 1.0, 400)

# trajectory array to store result #
trajectory = np.empty((len(times), len(init_cond)), dtype=np.float64)

# helper #
show_prog = True

# check eoms at starting position #
print(eoms_pendulum(0, init_cond, params))

################################################################################
### numerical integration ###
################################################################################

r = ode(eoms_pendulum).set_integrator('zvode', method='adams',
                                           with_jacobian=False) # integrator and eoms
r.set_initial_value(init_cond, times[0]).set_f_params(params)   # setup
dt = times[1] - times[0] # time step

# integration (loop time step) 
for i, t_i in enumerate(times):
    trajectory[i,:] = r.integrate(r.t+dt) # integration

### extract ###
x1 = trajectory[:, 3]
y1 = trajectory[:, 4]
x2 = np.zeros_like(trajectory[:, 5])
y2 = trajectory[:, 5]

L = np.sqrt(x1**2 + y1**2) # rope part connecting m1 and pivot
Ltot = -y2 + L             # total rope length


################################################################################
### Visualize trajectory ###
################################################################################
#
fig = plt.figure(figsize=(15,7))

plt.subplot(121)
titleStr = "m1: {}, m2: {}, g: {}, L: {}".format(params[0], params[1], params[2], Lrope)

fs = 8
plt.axvline(0, color='k', linewidth=1, dashes=[1,1])
plt.axhline(0, color='k', linewidth=1, dashes=[1,1])
plt.scatter(x1, y1, c=times, label="Mass 1")
plt.scatter(x2, y2, marker='x', c=times, label='Mass 2')
#plt.xlim(-1.5, 1.5)
#plt.ylim(-2, 1.)
plt.xlabel('x position', fontsize=fs)
plt.ylabel('y position', fontsize=fs)
plt.gca().set_aspect('equal', adjustable='box')
cbar = plt.colorbar()
cbar.ax.set_ylabel('Time', rotation=270, fontsize=fs)
plt.title(titleStr, fontsize=fs)
plt.legend()

plt.subplot(122)
plt.axhline(0., color='k', dashes=[1,1])
plt.plot(times, x1, '-', label="Mass 1, x pos")
plt.plot(times, y1, '-', label="Mass 1, y pos")
plt.plot(times, y2, '--', label="Mass 2, y pos")
plt.xlabel('Time')
plt.legend()

plt.tight_layout()
fig.savefig('{}-{}.pdf'.format(int(params[1]), int(params[0])))
plt.close()

# check that total length of the rope is constant #
plt.figure()

plt.axhline(0, color='k', linewidth=1, dashes=[1,1])
plt.axvline(0.4, color='k', linewidth=1, dashes=[1,1])
plt.plot(times, Ltot, label='total rope length')
plt.plot(times, L, label='rope from mass 1 to pivot')
plt.legend()
plt.tight_layout()

plt.close()

The dynamics for 1/14 mass ration

Here is what the dynamics of the pendulum look like for a mass ratio of 14 ($m_1 = 1$, $m_2=14$):

enter image description here

The left panel shows the trajectories of the two masses in the x-y plane, with time being indicated by the color axis. This is supposed to be a front view of the video performance. We see that mass 1 wraps around the pivot (at x=0, y=0) multiple times (see winding angle picture above). After a few revolutions, the model is probably not representative anymore. Instead, friction would start kicking in and clip the rope. In our simplified picture, the particle keeps going. What is interesting is that even without friction, the lower particle stops at some point and even comes back up, causing stable oscillation!!

What changes with the mass ratio?

We already saw what changes with the mass ratio in the winding angle picture. Just for visual intuiution, here is the corresponding picture for 1/5 mass ratio:enter image description here

For higher mass ratio (1/20):

enter image description here

$\endgroup$
6
  • 2
    $\begingroup$ Would be nice to see animations in addition to static trajectories. Something like this animation for the first trajectory $\endgroup$
    – Ruslan
    Mar 25, 2020 at 21:23
  • $\begingroup$ @Ruslan That is an excellent idea!!! Super nice animation too, how did you make it? Did you use the same equations of motion as in the answer? I'll figure something out regarding the animations, I'm sure there is a matplotlib solution. I think a whole panel with the movies for the different mass ratios running in parallel would be nice. Feel free to post an answer yourself though that I can upvote ;) $\endgroup$ Mar 25, 2020 at 21:58
  • 1
    $\begingroup$ Well, I just dumped the solution arrays from your code, reformatted them as CSV and imported into Wolfram Mathematica. I'm not very familiar with various Python packages, while Mathematica is what I normally use for such computations. There, Table, Graphics and Export to GIF were the functions I used. $\endgroup$
    – Ruslan
    Mar 25, 2020 at 22:02
  • 1
    $\begingroup$ @Ruslan check it out! Made in python using the matplotlib.animation.FuncAnimation tool :) $\endgroup$ Mar 27, 2020 at 1:08
  • 1
    $\begingroup$ Nice. Would be great if you posted the code you used to generate the animations. As for animations themselves, it might be better to supersample to avoid the jumpy points when they move by less than one pixel per frame. A simple way to do it is render in 5× larger resolution (both horizontally and vertically) and downsample before putting the frame into the GIF sequence. See these Q&A for more discussion of the problem and the solution (that's with Mathematica, but the idea is the same for any renderer). $\endgroup$
    – Ruslan
    Mar 27, 2020 at 5:38

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.