The Wilberforce Pendulum


The Wilberforce Pendulum is a mass suspended on a helical spring which is free to move up and down as well as to rotate as the spring coils and uncoils with its extension. You can see one in action below.

An analysis of the motion of a Wilberforce pendulum was carried out by Berg and Marshall ["Wilberforce pendulum oscillations and normal modes", Am. J. Phys. 59 (1991)] who give the following Lagrangian equations of motion for the extension, $z$, and torsion, $\theta$ of the spring:

\begin{align*} m\frac{\mathrm{d}^2z}{\mathrm{d}t^2} + kz + \frac{1}{2}\epsilon\theta & = 0,\\ I\frac{\mathrm{d}^2\theta}{\mathrm{d}t^2} + \delta\theta + \frac{1}{2}\epsilon z & = 0, \end{align*}

where $k$ and $\delta$ are the longitudinal and torsional spring constants respectively and $m$ and $I$ are the mass and moment of inertia of the suspended bob. $\epsilon$ is the coupling constant between the linear and twisting motions.

These equations may be put into canonical form as \begin{align*} \frac{\mathrm{d}\dot{z}}{\mathrm{d}t} = -\omega^2z - \frac{\epsilon}{2m}\theta & \quad \frac{\mathrm{d}z}{\mathrm{d}t} = \dot{z},\\ \frac{\mathrm{d}\dot{\theta}}{\mathrm{d}t} = -\omega^2\theta - \frac{\epsilon}{2I}z & \quad \frac{\mathrm{d}\theta}{\mathrm{d}t} = \dot{\theta}, \end{align*} where $\omega = \omega_z = \omega_\theta$ for resonance between the two modes with frequencies $\omega_z = \sqrt{k/m}$ and $\omega_\theta = \sqrt{\delta/I}$.

To simulate this system using SciPy, we use scipy.integrate.odeint in the code below, which generates plots of $z$ and $\theta$ against time as well as cartesian and polar plots of $z$ against $\theta$.

This code is also available on my github page.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Parameters for the system
omega = 2.314       # rad.s-1
epsilon = 9.27e-3   # N
m = 0.4905          # kg
I = 1.39e-4         # kg.m2

def deriv(y, t, omega, epsilon, m, I):
    """Return the first derivatives of y = z, zdot, theta, thetadot."""
    z, zdot, theta, thetadot = y
    dzdt = zdot
    dzdotdt = -omega**2 * z - epsilon / 2 / m * theta
    dthetadt = thetadot
    dthetadotdt = -omega**2 * theta - epsilon / 2 / I * z
    return dzdt, dzdotdt, dthetadt, dthetadotdt

# The time grid in s
t = np.linspace(0,40,50000)
# Initial conditions: theta=2pi, z=zdot=thetadot=0
y0 = [0, 0, 2*np.pi, 0]

# Do the numerical integration of the equations of motion
y = odeint(deriv, y0, t, args=(omega, epsilon, m, I))
# Unpack z and theta as a function of time
z, theta = y[:,0], y[:,2]

# Plot z vs. t and theta vs. t on axes which share a time (x) axis
fig, ax_z = plt.subplots(1,1)
l_z, = ax_z.plot(t, z, 'g', label=r'$z$')
ax_z.set_xlabel('time /s')
ax_z.set_ylabel(r'$z /\mathrm{m}$')
ax_theta = ax_z.twinx()
l_theta, = ax_theta.plot(t, theta, 'orange', label=r'$\theta$')
ax_theta.set_ylabel(r'$\theta /\mathrm{rad}$')

# Add a single legend for the lines of both twinned axes
lines = (l_z, l_theta)
labels = [line.get_label() for line in lines]
plt.legend(lines, labels)

# Plot theta vs. z on a cartesian plot
fig, ax1 = plt.subplots(1,1)
ax1.plot(z, theta, 'r', alpha=0.4)
ax1.set_xlabel(r'$z /\mathrm{m}$')
ax1.set_ylabel(r'$\theta /\mathrm{rad}$')

# Plot z vs. theta on a polar plot
fig, ax2 = plt.subplots(1,1, subplot_kw={'projection': 'polar'})
ax2.plot(theta, z, 'b', alpha=0.4)

The plots generated are shown below.

The Wilberforce Pendulum: z and theta vs. time The Wilberforce Pendulum: theta vs z (cartesian plot) The Wilberforce Pendulum: z vs theta (polar plot)

Current rating: 5


Comments are pre-moderated. Please be patient and your comment will appear soon.

Jose Luis Castro 3 years, 1 month ago

Es un gran trabajo
y computacional !!!

Link | Reply
Currently unrated

christian 3 years, 1 month ago

Muchas gracias – I'm glad you found it interesting.

Link | Reply
Currently unrated

New Comment


required (not published)