Following on from the previous post on the double pendulum, here is a similar Python script for plotting the behaviour of the "spring pendulum": a bob of mass $m$ suspended from a fixed anchor by a massless spring.
As before, the code places frames of the animation in a directory frames/
which can be turned into a gif with, for example, ImageMagick's convert
utility. A derivation of the Lagrangian equations used is given under the code.
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
# Pendulum equilibrium spring length (m), spring constant (N.m)
L0, k = 1, 40
m = 1
# The gravitational acceleration (m.s-2).
g = 9.81
def deriv(y, t, L0, k, m):
"""Return the first derivatives of y = theta, z1, L, z2."""
theta, z1, L, z2 = y
thetadot = z1
z1dot = (-g*np.sin(theta) - 2*z1*z2) / L
Ldot = z2
z2dot = (m*L*z1**2 - k*(L-L0) + m*g*np.cos(theta)) / m
return thetadot, z1dot, Ldot, z2dot
# Maximum time, time point spacings and the time grid (all in s).
tmax, dt = 20, 0.01
t = np.arange(0, tmax+dt, dt)
# Initial conditions: theta, dtheta/dt, L, dL/dt
y0 = [3*np.pi/4, 0, L0, 0]
# Do the numerical integration of the equations of motion
y = odeint(deriv, y0, t, args=(L0, k, m))
# Unpack z and theta as a function of time
theta, L = y[:,0], y[:,2]
# Convert to Cartesian coordinates of the two bob positions.
x = L * np.sin(theta)
y = -L * np.cos(theta)
# Plotted bob circle radius
r = 0.05
# Plot a trail of the m2 bob's position for the last trail_secs seconds.
trail_secs = 1
# This corresponds to max_trail time points.
max_trail = int(trail_secs / dt)
def plot_spring(x, y, theta, L):
"""Plot the spring from (0,0) to (x,y) as the projection of a helix."""
# Spring turn radius, number of turns
rs, ns = 0.05, 25
# Number of data points for the helix
Ns = 1000
# We don't draw coils all the way to the end of the pendulum:
# pad a bit from the anchor and from the bob by these number of points
ipad1, ipad2 = 100, 150
w = np.linspace(0, L, Ns)
# Set up the helix along the x-axis ...
xp = np.zeros(Ns)
xp[ipad1:-ipad2] = rs * np.sin(2*np.pi * ns * w[ipad1:-ipad2] / L)
# ... then rotate it to align with the pendulum and plot.
R = np.array([[np.cos(theta), -np.sin(theta)],
[np.sin(theta), np.cos(theta)]])
xs, ys = - R @ np.vstack((xp, w))
ax.plot(xs, ys, c='k', lw=2)
def make_plot(i):
"""
Plot and save an image of the spring pendulum configuration for time
point i.
"""
plot_spring(x[i], y[i], theta[i], L[i])
# Circles representing the anchor point of rod 1 and the bobs
c0 = Circle((0, 0), r/2, fc='k', zorder=10)
c1 = Circle((x[i], y[i]), r, fc='r', ec='r', zorder=10)
ax.add_patch(c0)
ax.add_patch(c1)
# The trail will be divided into ns segments and plotted as a fading line.
ns = 20
s = max_trail // ns
for j in range(ns):
imin = i - (ns-j)*s
if imin < 0:
continue
imax = imin + s + 1
# The fading looks better if we square the fractional length along the
# trail.
alpha = (j/ns)**2
ax.plot(x[imin:imax], y[imin:imax], c='r', solid_capstyle='butt',
lw=2, alpha=alpha)
# Centre the image on the fixed anchor point, and ensure the axes are equal
ax.set_xlim(-np.max(L)-r, np.max(L)+r)
ax.set_ylim(-np.max(L)-r, np.max(L)+r)
ax.set_aspect('equal', adjustable='box')
plt.axis('off')
plt.savefig('frames/_img{:04d}.png'.format(i//di), dpi=72)
# Clear the Axes ready for the next image.
plt.cla()
# Make an image every di time points, corresponding to a frame rate of fps
# frames per second.
# Frame rate, s-1
fps = 10
di = int(1/fps/dt)
# This figure size (inches) and dpi give an image of 600x450 pixels.
fig = plt.figure(figsize=(8.33333333, 6.25), dpi=72)
ax = fig.add_subplot(111)
for i in range(0, t.size, di):
print(i // di, '/', t.size // di)
make_plot(i)
There are two degrees of freedom in this problem, which are taken to be the angle of the pendulum from the vertical and the total length of the spring. Let the origin be at the anchor with the $y$-axis pointing up and the $x$ axis pointing to the right. The spring is assumed to have a force constant $k$ and an equilibrium length $l_0$ such that its potential energy when extended or compressed to a length $l$ is $\frac{1}{2}k(l-l_0)^2$. The components of the bob's position and velocity are:
The potential and kinetic energies are therefore:
and so the Lagrangian is
Applying the Euler-Lagrange relations for $\theta$ and $l$ as before gives
To cast this pair of second-order differential equations into four first-order differential equations, let $z_1 = \dot{\theta}$ and $z_2 = \dot{l}$ to give
which are the relations that appear in the function deriv
.
Comments
Comments are pre-moderated. Please be patient and your comment will appear soon.
Miguel 5 years, 4 months ago
I try to make the simulation of the pendulum, but it does not give me a shame, please I need you to help me
Link | Replychristian 5 years, 4 months ago
Hi Miguel,
Link | ReplyCould you be a little more specific about what the problem is, exactly?
I'd like to help,
Christian
james 5 years ago
hello, does this code use the Euler Lagrangian method or runge kutta? thank you also for the code :)
Link | Replychristian 5 years ago
It solves the Euler-Lagrange equations using scipy.integrate.odeint, which is a Python wrapper to the LSODA algorithm (which is not a Runge-Kutta method).
Link | Replysue 2 years, 6 months ago
I copied your code in colab
Link | Replybut it didn't work
error: FileNotFoundError: [Errno 2] No such file or directory: 'frames/_img0000.png'
christian 2 years, 6 months ago
Dear Sue, You need to make sure that the directory frames/ exists in the same location as the script – the code doesn't create it for you.
Link | ReplyI hope that helps,
Christian
Diego Carrera 6 months, 1 week ago
hey can you help me? I'm new in python and I don't know what to do with the error error: FileNotFoundError: [Errno 2] No such file or directory: 'frames/_img0000.png'. I understand why this happen, but I don't know how to solve it. Sorry my bad English haha :p
Link | Replychristian 6 months, 1 week ago
Do you have a directory called frames in the directory of your script? That's where the animation frames get written to.
Link | ReplyNew Comment