# Spherical Pendulums

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

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
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)
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']

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']
previouspsi = course.loc[0,'psi']
previouspsidot = course.loc[0,'psidot']
previouspsidotdot = course.loc[0,'psidotdot']

for timeStep in timeSteps[1:]:
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

psidot = previouspsidot + previouspsidotdot * lengthOfTimeStep
course.loc[timeStep, 'psidot'] = psidot

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

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

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

```