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:

\begin{align*}
x &= l\sin\theta & \dot{x} &= \dot{l}\sin\theta + l\dot{\theta}\cos\theta\\
y &= -l\cos\theta & \dot{y} &= -\dot{l}\cos\theta + l\dot{\theta}\sin\theta
\end{align*}

The potential and kinetic energies are therefore:

\begin{align*}
V &= mgy + \tfrac{1}{2}k(l-l_0)^2 = -mgl\cos\theta + \tfrac{1}{2}k(l-l_0)^2,\\
T &= \tfrac{1}{2}mv^2 = \tfrac{1}{2}m(\dot{x}^2 + \dot{y}^2) = \tfrac{1}{2}m(\dot{l}^2 + l^2\dot{\theta}^2),
\end{align*}

and so the Lagrangian is

\begin{align*}
\mathcal{L} &= \tfrac{1}{2}m(\dot{l}^2 + l^2\dot{\theta}^2) - \tfrac{1}{2}k(l-l_0)^2 + mgl\cos\theta.
\end{align*}

Applying the Euler-Lagrange relations for $\theta$ and $l$ as before gives

\begin{align*}
2\dot{l}\dot{\theta} + l\ddot{\theta} + g\sin\theta & = 0\\
m\ddot{l} - ml\dot{\theta}^2 + k(l-l_0) - mg\cos\theta &= 0.
\end{align*}

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

\begin{align*}
\dot{z}_1 &= -\frac{1}{l}[g\sin\theta + 2z_1z_2],\\
\dot{z}_2 &= \frac{1}{m}[mlz_1^2 - k(l-l_0) + mg\cos\theta].
\end{align*}

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 1 year, 8 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 | Reply## christian 1 year, 8 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 1 year, 4 months ago

hello, does this code use the Euler Lagrangian method or runge kutta? thank you also for the code :)

Link | Reply## christian 1 year, 4 months 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 | Reply## New Comment