text: 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) 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.
  Share on Twitter
  Share on Facebook
  Comments are pre-moderated. Please be patient and your comment will appear soon.
  Jose Luis Castro 3 years, 1 month ago
  Link | Reply
  Es un gran trabajo
  experimental
  teorico
  y computacional !!!
  christian 3 years, 1 month ago
  Muchas gracias – I'm glad you found it interesting.Link | Reply
  New Comment