import numpy as np
from scipy import optimize
%matplotlib inline
import matplotlib.pyplot as plt
class Chain:
def __init__(self, nlinks, hdist):
if nlinks < hdist:
raise ValueError('a hanging chain with {} links of '
'length 1 cannot have endpoints '
'separated by a distance {}'.format(
nlinks, hdist))
self.nlinks = nlinks
self.hdist = hdist
def equilibrium(self):
phi0 = np.arcsin(self.hdist/self.nlinks)
result = optimize.minimize(
self.f_energy,
np.linspace(-phi0, phi0, self.nlinks),
method='SLSQP',
constraints=[{'type': 'eq', 'fun': self.x_constraint},
{'type': 'eq', 'fun': self.y_constraint}])
return result.x
def x_constraint(self, phi_vals):
"""ensure the correct horizontal distance
phi_vals: angles of links with respect to the horizontal
"""
return np.sum(np.cos(phi_vals))-self.hdist
def y_constraint(self, phi_vals):
"""ensure that endpoints are at the same height
phi_vals: angles of links with respect to the horizontal
"""
return np.sum(np.sin(phi_vals))
def f_energy(self, phi_vals):
"""potential energy of all links
"""
return np.sum(np.arange(self.nlinks, 0, -1)*np.sin(phi_vals))
def angles_to_coords(phi_vals):
"""convert angles to coordinates of link endpoints
phi_vals: angles of links with respect to the horizontal
x, y: coordinates of link endpoints
"""
dim = phi_vals.shape[0]+1
x = np.zeros(dim)
x[1:] = np.cumsum(np.cos(phi_vals))
y = np.zeros(dim)
y[1:] = np.cumsum(np.sin(phi_vals))
return x, y
nlinks = 20
hdist = 16
c = Chain(nlinks, hdist)
plt.axes().set_aspect('equal')
_ = plt.plot(*angles_to_coords(c.equilibrium()), 'o-')
W. Tomaszewski, P. Pieranski, J.-C. Geminard
The motion of a free falling chain tip
American Journal of Physics 74, 776 (2006)
with $$a_i = n-i+\frac{1}{2}$$ and $$m_{i,j} = \begin{cases} n-i+\frac{1}{3} & i=j \\ n-\max(i,j)+\frac{1}{2} & i\neq j \end{cases}$$
parameters:
$m$ mass of chain element
$\ell$ length of chain element
$g$ acceleration due to gravity
$r$ damping constant
Introduce a new variable $$p_\varphi = \dot\varphi$$ to replace the second-order differential equation $$\ddot\varphi = f(\dot\varphi, \varphi)$$ by two first-order differential equations $$ \begin{aligned} \dot\varphi &= p_\varphi\\ \dot p_\varphi &= f(p_\varphi, \varphi) \end{aligned} $$
import numpy.linalg as LA
from scipy import integrate
class FallingChain(Chain):
def __init__(self, nlinks, hdist, damping):
super(FallingChain, self).__init__(nlinks, hdist)
self.set_matrix_m()
self.set_vector_a()
self.set_matrix_damping(damping)
def set_matrix_m(self):
m = np.fromfunction(lambda i, j: self.nlinks-np.maximum(i, j)-0.5,
(self.nlinks, self.nlinks))
self.m = m-np.identity(self.nlinks)/6
def set_vector_a(self):
self.a = np.arange(self.nlinks, 0, -1)-0.5
def set_matrix_damping(self, damping):
self.damping = (-2*np.identity(self.nlinks, dtype=np.float64)
+ np.eye(self.nlinks, k=1)
+ np.eye(self.nlinks, k=-1))
self.damping[0, 0] = -1
self.damping[self.nlinks-1, self.nlinks-1] = -1
self.damping = damping*self.damping
def solve_eq_of_motion(self, time_i, time_f, nt):
y0 = np.zeros(2*self.nlinks, dtype=np.float64)
y0[self.nlinks:] = self.equilibrium()
self.solution = integrate.solve_ivp(
self.diff, (time_i, time_f), y0,
method='RK45',
t_eval=np.linspace(time_i, time_f, nt))
def diff(self, t, y):
momenta = y[:self.nlinks]
angles = y[self.nlinks:]
d_angles = momenta
ci = np.cos(angles)
cij = np.cos(angles[:, np.newaxis]-angles)
sij = np.sin(angles[:, np.newaxis]-angles)
mcinv = LA.inv(self.m*cij)
d_momenta = -np.dot(self.m*sij, momenta*momenta)
d_momenta = d_momenta + np.dot(self.damping, momenta)
d_momenta = d_momenta - self.a*ci
d_momenta = np.dot(mcinv, d_momenta)
d = np.empty_like(y)
d[:self.nlinks] = d_momenta
d[self.nlinks:] = d_angles
return d
from matplotlib import rc
import matplotlib.animation as animation
rc('animation', html='jshtml')
The following cell should always be executed before running the animation cell below. In this way, an undesired extra static image will be avoided.
%%capture
fig = plt.figure()
def animate(i):
x, y = angles_to_coords(c.solution.y[:, i][c.nlinks:])
line.set_data(x, y)
return line,
def init():
line.set_data([], [])
return line,
nlinks = 50
hdist = 40
damping = 1
ti = 0
tf = 200
tsteps = 1000
c = FallingChain(nlinks, hdist, damping)
c.solve_eq_of_motion(ti, tf, tsteps)
ax = fig.add_subplot(111, autoscale_on=False,
xlim=(-nlinks, nlinks), ylim=(-nlinks, 0.3*nlinks))
ax.set_aspect('equal')
line, = ax.plot([], [])
anim = animation.FuncAnimation(fig, animate, tsteps,
interval=20, blit=True,
init_func=init)
anim