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:

simple-pendulum-with-time-no-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:

simple-pendulum-with-time-and-air-resistance-2

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:

spherical-pendulum-no-air-resistance

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:

spherical-pendulum-no-air-resistance-1-air-resistance-comparison

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

spherical-pendulum-with-0-02-air-resistance-long

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):

spherical-pendulum-no-air-resistance-5-air-resistance-comparison

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

spherical-pendulum-with-0-05-air-resistance-long

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

def main():
 print("Starting spherical pendulum")
 t0 = time.time()
 [course,constants,timeSteps] = inputInitialConditions()
 course = completeInitialTimeStep(course, constants, timeSteps)
 course = calculatePendulumCourse(course, constants, timeSteps)
 #course.to_csv('C:\Users\david\Desktop\spherical_pendulum_0.05_air_resistance.csv')
 print (time.time()-t0)
 plotPendulumCourse(course)
 print("Completed spherical pendulum")


main()
</pre>
<pre>

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')
plt.axis('equal')
plt.show()

Leave a Reply