Spherical Pendulums

Finding the governing equations, and forming a strategy for creating a really basic solver given some initial conditions:

Download (PDF, 122KB)

Then using this code to draw some pendulums and their movements through time:

For a simple pendulum, plotting its x and y coordinated with time, without any air resistance:


Then including some air resistance, see pdf for detail, shows a roughly exponential decay of the amplitude of the oscillations over time, as would be expected:


The same approach is then taken with the more complicated spherical pendulum, which, with no air resistance considered, has a path that looks like this:


Adding in a small air resistance of alpha=0.01 (see pdf for expanation), produces a slightly different path, as would be expected. The ‘no air resistance path’ is in blue, the ‘small air resistance path’ is in green:


Following the green path for a longer period of time, we see the following pattern (now plotted in blue):


Setting alpha=0.05 leads to the following difference from the no air resistance path (no air resistance in blue, with air resistance in green):


Following this green path for a longer period of time we see the following pattern (now plotted in blue rather than green):


Perhaps most nicely, the later, gentle, spiraling behaviour is clear here.


The code used to generate these plots was as follows, for detail on how it was derived see the pdf at the top of the post:

The simple pendulum:

import numpy
import matplotlib
import pandas
import matplotlib.pyplot as plt
import time
from mpl_toolkits.mplot3d import Axes3D

def inputInitialConditions():
 print("installising arrays and inputting initial conditions")
 t_initial = 0
 theta_initial = numpy.pi/5
 thetadot_initial = 0
 psi_initial = 0
 psidot_initial = 4
 numberOfTimeSteps = 20000
 lengthOfTimeStep = 0.001
 m = 1
 l = 1
 g = 9.81
 alpha = 0.04
 timeSteps = numpy.arange(numberOfTimeSteps)
 columnTitles = ['t', 'theta', 'thetadot', 'thetadotdot', 'psi', 'psidot', 'psidotdot', 'x', 'y', 'z']
 constants = {'m':m, 'l':l, 'g':g, 'lengthOfTimeStep':lengthOfTimeStep, 'numberOfTimeSteps':numberOfTimeSteps, 't_initial':t_initial, 'alpha':alpha}
 course = pandas.DataFrame(numpy.zeros([numberOfTimeSteps,10],dtype=long), index = timeSteps, columns = columnTitles)
 course.loc[0] = [t_initial, theta_initial,thetadot_initial,0,psi_initial,psidot_initial,0,0,0,0]
 return [course, constants, timeSteps]

def completeInitialTimeStep(course, constants, timeSteps):
 g = constants['g']
 l = constants['l']
 m = constants['m']
 alpha = constants['alpha']
 t_initial = constants['t_initial']
 theta = course.loc[0,'theta']
 psi = course.loc[0,'psi']
 if numpy.absolute(theta) < 0.000001:
 theta = 0.000001
 print 'Important Warning: Very small theta on timestep 0 generating errors, use simple pendulum program'
 if numpy.absolute(theta) > numpy.pi/2:
 theta = numpy.pi/2 
 print 'Important Warning: Very large theta, potential errors or unphysical behaviour'
 sintheta = numpy.sin(theta)
 costheta = numpy.cos(theta)
 tantheta = sintheta / costheta
 sinpsi = numpy.sin(psi)
 cospsi = numpy.cos(psi)
 psidot = course.loc[0, 'psidot']
 thetadot = course.loc[0,'thetadot']

 course.loc[0,'thetadotdot'] = ((psidot**2)*costheta-g/l)*sintheta-alpha*l/m*sintheta*numpy.absolute(thetadot)*thetadot
 course.loc[0,'psidotdot'] = -2*psidot*thetadot/tantheta-alpha*l/m*numpy.absolute(psidot)*psidot
 course.loc[0,'t'] = t_initial
 course.loc[0,'x'] = l*sintheta*cospsi
 course.loc[0,'y'] = l*sintheta*sinpsi
 course.loc[0,'z'] = l*costheta
 return course

def calculatePendulumCourse(course, constants, timeSteps):
 g = constants['g']
 l = constants['l']
 m = constants['m']
 alpha = constants['alpha']
 lengthOfTimeStep = constants['lengthOfTimeStep']
 previoustheta = course.loc[0,'theta']
 previousthetadot = course.loc[0, 'thetadot']
 previousthetadotdot = course.loc[0, 'thetadotdot']
 previouspsi = course.loc[0,'psi']
 previouspsidot = course.loc[0,'psidot']
 previouspsidotdot = course.loc[0,'psidotdot']
 for timeStep in timeSteps[1:]:
 theta = previoustheta + (previousthetadot * lengthOfTimeStep) + ((previousthetadotdot * lengthOfTimeStep**2)/2)
 course.loc[timeStep, 'theta'] = theta
 sintheta = numpy.sin(theta)
 costheta = numpy.cos(theta)
 tantheta = numpy.tan(theta)
 if numpy.absolute(theta) < 0.000001:
 sintheta = 0.000001
 tantheta = sintheta
 print 'Important Warning: Very small theta on timestep' + str(timeStep)
 psi = course.loc[timeStep, 'psi'] = previouspsi + (previouspsidot * lengthOfTimeStep) + ((previouspsidotdot * lengthOfTimeStep**2)/2)
 course.loc[timeStep, 'psi'] = psi
 sinpsi = numpy.sin(psi)
 cospsi = numpy.cos(psi)
 if numpy.absolute(psi) < 0.000001:
 sinpsi = 0.000001
 print 'Important Warning: Very small psi on timestep' + str(timeStep)
 tanpsi = sinpsi / cospsi 
 thetadot = previousthetadot + previousthetadotdot * lengthOfTimeStep
 course.loc[timeStep, 'thetadot'] = thetadot
 psidot = previouspsidot + previouspsidotdot * lengthOfTimeStep
 course.loc[timeStep, 'psidot'] = psidot 
 thetadotdot = ((psidot**2)*costheta-g/l)*sintheta-alpha*l/m*numpy.absolute(thetadot)*thetadot
 course.loc[timeStep,'thetadotdot'] = thetadotdot
 psidotdot = -2*psidot*thetadot/tantheta-alpha*l/m*numpy.absolute(psidot)*psidot*sintheta
 course.loc[timeStep,'psidotdot'] = psidotdot
 x = l*sintheta*cospsi
 y = l*sintheta*sinpsi
 z = l*costheta
 course.loc[timeStep,'x'] = x
 course.loc[timeStep,'y'] = y
 course.loc[timeStep,'z'] = z
 previoustheta = theta
 previousthetadot = thetadot
 previousthetadotdot = thetadotdot
 previouspsi = psi
 previouspsidot = psidot
 previouspsidotdot = psidotdot
 return course
def plotPendulumCourse(course):
 xs = course['x']
 ys = course['y']
 zs = course['z']
 fig = plt.figure()
 ax = fig.gca(projection='3d')
 ax.plot(xs, ys, zs)

def main():
 print("Starting spherical pendulum")
 t0 = time.time()
 [course,constants,timeSteps] = inputInitialConditions()
 course = completeInitialTimeStep(course, constants, timeSteps)
 course = calculatePendulumCourse(course, constants, timeSteps)
 print (time.time()-t0)
 print("Completed spherical pendulum")


The comparison plot:

import numpy
import matplotlib
import pandas
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

course0 = pandas.read_csv('C:\Users\david\Desktop\spherical_pendulum_no_air_resistance.csv')
course1 = pandas.read_csv('C:\Users\david\Desktop\spherical_pendulum_0.01_air_resistance.csv')
course5 = pandas.read_csv('C:\Users\david\Desktop\spherical_pendulum_0.05_air_resistance.csv')
course10 = pandas.read_csv('C:\Users\david\Desktop\spherical_pendulum_0.1_air_resistance.csv')

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(course0['x'], course0['y'], course0['z'], label = 'alpha = 0')
ax.plot(course1['x'], course1['y'], course1['z'], label = 'alpha = 0')
ax.plot(course5['x'], course5['y'], course5['z'], label = 'alpha = 0')
ax.plot(course10['x'], course10['y'], course10['z'], label = 'alpha = 0')
ax.set_xlabel('x value')
ax.set_ylabel('y value')
ax.set_zlabel('z value')

Leave a Reply