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