Three springs and a mass go to twitter
All good posts start with a tweet. How about this one?
I invented this problem at breakfast and it's remarkably difficult: a mass m inside an equilateral triangle is connected by springs. Find out how it will move.
— Maarten Mortier (@maartengm) May 17, 2022
I failed: could not eliminate z (not even with mentioned relationship), nor found a good coordinate system. Tbc.. pic.twitter.com/CmSgIOGVwP
Let’s solve this using Python. First, some drawings:
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
degree = np.pi / 180.
anchors = np.array([[0, 0],
[1, 0],
[np.cos(60 * degree), np.sin(60 * degree)]])
def draw_scene(x, y, ax):
for anchor in anchors:
ax.plot([anchor[0], x],
[anchor[1], y], "--k")
ax.plot(anchors[:, 0], anchors[:, 1], 'o')
ax.plot(x, y, 'o', ms=10)
fig, ax = plt.subplots()
draw_scene(0.5 * np.cos(30 * degree), 0.5 * np.sin(30 * degree), ax)
Let’s now introduce the equations that time-step the mass that is subject to the spring forces.
In [2]:
def spring_force(start, stop, k, l0):
AB = stop - start
l = np.linalg.norm(AB)
return -k * (l - l0) * AB / l
spring_force(np.array([0, 0]), np.array([1, 0]), 10, 0.5)
Out[2]:
array([-5., -0.])
In [3]:
def step(x, xdot, ks, l0s, anchors, dt):
force = np.zeros((2,))
for start, k, l0 in zip(anchors, ks, l0s):
force += spring_force(start, x, k, l0)
xdot += force * dt
x += xdot * dt
return x, xdot
x = np.array([0.5 * np.cos(30 * degree), 0.5 * np.sin(30 * degree)])
xdot = np.array([0, 0.])
step(x, xdot, ks=[1, 2, 3], l0s=[2, 4, 6], anchors=anchors, dt=0.1)
Out[3]:
(array([0.36669359, 0.12431156]), array([-0.66319116, -1.25688438]))
Let’s try some trajectories.
In [4]:
x = np.array([0.8, 0.8])
xdot = np.array([0.9, 0.6])
trajectory = []
for _ in range(700):
x, xdot = step(x, xdot, ks=[1, 2, 3], l0s=[0.6, 0.8, 0.5], anchors=anchors, dt=0.1)
trajectory.append(x.copy())
trajectory = np.array(trajectory)
fig, ax = plt.subplots()
ax.plot(trajectory[:, 0], trajectory[:, 1])
draw_scene(x[0], x[1], ax)
This looks good. Let’s make an animation.
In [5]:
from IPython.display import HTML
import matplotlib.animation as manim
fig, ax = plt.subplots()
def func(frame):
ax.cla()
ax.plot(trajectory[:int(frame), 0], trajectory[:int(frame), 1])
draw_scene(trajectory[int(frame), 0],
trajectory[int(frame), 1], ax)
animation = manim.FuncAnimation(fig, func, frames=np.linspace(0, 700, 100, dtype=int, endpoint=False))
plt.close(fig)
HTML(animation.to_jshtml())
Out[5]: