# 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
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
dthetadotdt = -omega**2 * theta - epsilon / 2 / I * z

# The time grid in s
t = np.linspace(0,40,50000)
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)
plt.show()

# 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}$')
plt.show()

# 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)
plt.show()


The plots generated are shown below.   Current rating: 5 #### Jose Luis Castro 3 years, 1 month ago

Es un gran trabajo
experimental
teorico
y computacional !!!

Currently unrated #### christian 3 years, 1 month ago

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