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