Say we have the simplest pendulum imaginable -- one in which we ignore friction and air resistance, and we assume the rod connecting our mass to its pivot point is rigid and massless. Say the rod has some length $l$, which is connected to a bob of mass $m$, and we displace it from the horizontal by some angle $\theta$ before letting it go. We are curious about the motion of the pendulum. To investigate what happens, we will derive the system's equations of motion, then numerically solve the equations using Python.
We begin by drawing and investigating a free body diagram, seen below. Note that we use polar coordinates instead of cartesian coordinates to simplify the problem greatly.
Free Body Diagram of simple pendulum
To derive the equations of motion, we start by setting up Newton's second law in the radial and anagular directions. First, let us consider the radial direction. Since we are assuming the rod of our pendulum is perfectly rigid, this means its length $l$ is a constant, and our mass will have no acceleration in the this direction. This means that the "tension" force $T_r$ exerted on the mass by the rod is always equal and opposite to the radial component of our weight $F_r$
$$
\sum F_r = m\frac{d^2r}{dt^2} = 0
$$
Now we consider the angular direction. If we let α be our angular acceleration, setting up Newton's second law yields
$$
\sum F_\theta = ma_\theta = ml\alpha = ml\frac{d^2\theta}{dt^2} = -mg\sin\theta
$$
\begin{equation}
\Rightarrow \frac{d^2\theta}{dt^2}=-\frac{g}{l}\sin\theta = -\omega_0^2\sin\theta
\end{equation}
Where g=9.81m/s2 is acceleration due to gravity, and we have defined the angular frequency (which interestingly is only a function of the pendulum length, since we might consider g a constant on the surface of the Earth)to be
$$
\omega_0 = \sqrt{\frac{g}{l}}
$$
Our second order differential equation can easily be solved numerically, but has a rather disgusting analytic solution involving Jacobian elliptic functions. Before numerically solving the exact equation, we will investigate what happens when we use the small angle approximation. This way we might compare the two methods of treating this situation, and see how well the small angle approximations holds.
For the uninitiated, the small angle approximation says that for small angles (in radians),
$$
\sin\theta\approx\theta
$$
This is by no means a rigorous proof, but to convice yourself this is a decent approximation, consider the following table:
θ (degrees) |
θ (radians) |
sinθ
0
0
0
1
0.0175
0.0175
5
0.0873
0.0872
10
0.1745
0.1736
15
0.2618
0.2588
20
0.3491
0.3420
25
0.4363
0.4226
30
0.5236
0.5000
35
0.6109
0.5736
40
0.6981
0.6428
45
0.7854
0.7071
50
0.8727
0.7660
90
1.5708
1.0000
120
2.0944
0.8660
It is clear that the approximation is quite good for small angles, and gets worse the larger your angle is (shocking result, given its name). This means that simplifying our differential equation to
$$
\frac{d^2\theta}{dt^2}=\omega_0^2\theta
$$
will act as a good approximation, as long as we keep the initial angular displacement of our pendulum fairly small. With this simplification, our differential equation is quite easy to solve analytically (guess and check method is sufficient). This yields the result
$$
\theta(t) = \theta_{max}\sin(\omega_0 t)
$$
where θmax is our initial angular displacement.
We can make an estimation of what we should expect the period T of our pendulum to be. Recalling the relationship between frequency f and angular frequency,
$$
f = \frac{1}{2 \pi}\omega_0 = \frac{1}{2\pi}\sqrt{\frac{g}{l}}
$$
$$
T = \frac{1}{f} = 2\pi\sqrt{\frac{l}{g}}
$$
So if we let the length of our pendulum be l=0.25m, we should expect our period to be almost exactly equal to 1 second. Below, find the code I wrote to simulate our pendulum when dropped from an initial angle of 10 degrees. Below that, find a plot of our angle vs. time, and lets see if this prediction for the period is accurate.
# Use the small angle approximation to plot theta vs time for a simple pendulum# Import packages neededimportmatplotlib.animationasanimationfromnumpyimportsqrt,sin,cos,arange,pi,appendfrompylabimportplot,xlabel,ylabel,title,show,axhline,savefig,subplots_adjust,\
figure,xlim,rcParams,rc,rc_context,subplot,tight_layout,axvline# define a function to take in an initial angle and a time, then outputs the angle after a timestepdefsimpPendApprox(theta_max,t):# calculate our angle at time t, don't forget to add our phase!theta_new=theta_max*sin(omega0*t+pi/2)# calculate our angular speed at time tomega_new=omega0*theta_max*cos(omega0*t+pi/2)# return the angle and angular speedreturn(theta_new,omega_new)# set up the parameters of our situationl=0.25# length of pendulumg=9.81# acceleration due to gravityomega0=sqrt(g/l)# resonant frequencytheta_initial_deg=10# initial angle in degreestheta_initial=theta_initial_deg*pi/180# convert initial angle into degrees# set up a domain (time interval of interest)a=0.0# interval startb=16# interval enddt=0.005# timestept_points=arange(a,b,dt)# array of times# create a list of angles and angular speeds to fill then plotangle_points=[]omega_points=[]fortint_points:next_angle,next_omega=simpPendApprox(theta_initial,t)angle_points.append(next_angle)omega_points.append(next_omega)# convert angles back to degrees, and put angular speed in deg/sangle_points_deg=[]omega_points_deg=[]foriinrange(len(t_points)):angle_points_deg.append(angle_points[i]*180/pi)omega_points_deg.append(omega_points[i]*180/pi)##### Now to plot/animate ###### begin with some styling options# fontsize of titles and labelsrcParams.update({'font.size':14})# make our axes a little thickerrc('axes',linewidth=1.5)# update the colors of... pretty much everything!withrc_context({'axes.edgecolor':'white','xtick.color':'white', \
'ytick.color':'white','figure.facecolor':'darkslategrey',\
'axes.facecolor':'darkslategrey','axes.labelcolor':'white',\
'axes.titlecolor':'white'}):# create a figure with an appropriate sizefig1=figure(figsize=(10,9))# make our figures less cramped (ironic name for such a function)fig1.tight_layout()subplots_adjust(wspace=0.3,hspace=0.3)# create subplot for our pendulum animation# specify that it is a 2x2 grid, in top leftax_pend=subplot(2,2,1,aspect='equal',adjustable='datalim')title("Pendulum (approx)")# make axis ticks same color as backgroundax_pend.tick_params(axis='both',colors="darkslategrey")# create a suplot for our phase diagram, in top rightax_phase=subplot(2,2,2,aspect='equal',adjustable='datalim')title("Angular Speed vs. Angle")xlabel("Angle (deg)")ylabel("Angular Speed (deg/s)")axhline(color="k",lw=1)axvline(color="k",lw=1)# create a subplot for our angle vs. time graph# specify that we want it to be two rows longax_ang=subplot(2,2,(3,4))title("Angle Vs. Time")xlabel("Time (s)")ylabel("Angle (degrees)")axhline(color="k",lw=1)##### Now we animate ###### create an empty list to fill with pictures of our plot each step# set an index to keep track of how many steps we've takenims=[]index=0whileindex<len(t_points):### our pendulum ### # the rigid rod of length lln,=ax_pend.plot([0,l*sin(angle_points[index])],[0,-l*cos(angle_points[index])],\
color='k',lw=3)# our mass on the endbob,=ax_pend.plot(l*sin(angle_points[index]),-l*cos(angle_points[index]),'o',\
markersize=22,color="magenta")### phase diagram ###phase_curve,=ax_phase.plot(angle_points_deg[:(index+1)],omega_points_deg[:(index+1)],\
color="coral")phase_dot,=ax_phase.plot(angle_points_deg[index:(index+1)],\
omega_points_deg[index:(index+1)],\
color="magenta",marker="o",markersize=14)### angle vs. time ###ang,=ax_ang.plot(t_points[:(index+1)],angle_points_deg[:(index+1)],\
color="magenta",lw=2.5)# add our pictures to the ims listims.append([ln,bob,phase_curve,phase_dot,ang])# only show every 2 time stepsindex+=2# save animationani=animation.ArtistAnimation(fig1,ims,interval=100)writervideo=animation.FFMpegWriter(fps=60)ani.save('../simplePendulum/10_simplePend0.mp4',writer=writervideo)
Looking at our angle vs. time graph, we can see that one period does indeed appear to last almost exactly one second. We also take our first look into the world of phase diagrams (top right), how exciting! A quick synopsis of phase diagrams is that they are a pictoral representation of all possible states a system can exist in. Here it shows every combination of speed (momentum) and position. Notice that the phase diagram is an ellipse. Moving forward on this page, we will use a pendulum length of 9.81, to make our phase diagram for our small-angle-approximated pendulum a circle! Note that this geometry is due to us treating the angular position (proportional to a sine term) and the angular speed (proportional to a cosine term) as a set of parametric equations. Letting the length of our pendulum equal 9.81 makes our ω0 term equal 1, so the coefficient in front of each trig term is just θmax. Thus our phase diagram must yield a circle with radius θmax. How exhilerating! See below how this change alters our pendulum system; specifically take a look at our new period.
This is great, but we want to push this approximation to its limits to find at what point it starts to be a poor decision to apply the small angle approximation. So, let us now solve the exact equation (without using the small angle approximation). We have
$$
\frac{d^2\theta}{dt^2} = \omega_0^2\sin\theta
$$
In order to solve this second order differential equation, we use a cute trick to turn it into two coupled first order differential equations. The trick is to let
$$
\frac{d\theta}{dt} = \omega
$$
such that,
$$
\frac{d\omega}{dt} = \omega_0^2\sin\theta = f(\theta,\omega,t) = f(\vec{\eta})
$$
Where we are arbitrarily defining a vector η=(θ,ω,t). We will be using the fourth order Runge-Kutta method so solve this system of equations, but you should check out my page on solving the equations of motion of different pendulum systems using both the Runge-Kutta method, and the Euler-Cromer method with different step sizes.
See the code block below to follow the coding process. Notice that the first chunk of code is almost identical to the code above, as we still want to gather information on the approximated pendulum.
A good coding exercise would be to turn the code below into a user defined function that takes the initial angle and initial angular velocity of the pendulums in as parameters. In this way, you can produce animations for different initial conditions merely by calling your function.
# solve the simple pendulum using fourth order Runge-Kutta method# Import packages neededimportmatplotlib.animationasanimationfromnumpyimportsqrt,sin,cos,arange,pi,append,arrayfrompylabimportplot,xlabel,ylabel,title,show,axhline,savefig,subplots_adjust,\
figure,xlim,rcParams,rc,rc_context,subplot,tight_layout,axvline,legend##### First we solve the simple case ###### define a function to take in an initial angle and a time, then outputs the angle after a timestepdefsimpPendApprox(theta_max,t):# calculate our angle at time t, don't forget to add our phase!theta_new=theta_max*sin(omega0*t+pi/2)# calculate our angular speed at time tomega_new=omega0*theta_max*cos(omega0*t+pi/2)# return the angle and angular speedreturn(theta_new,omega_new)# an integer to scale our intial angle, and the length of our pendulum# for animatingn=1# set up the parameters of our situationl=9.81# length of pendulumg=9.81# acceleration due to gravityomega0=sqrt(g/l)# resonant frequencytheta_initial_deg=n*10# initial angle in degreestheta_initial=theta_initial_deg*pi/180# convert initial angle into degrees# set up a domain (time interval of interest)a=0.0# interval startb=38# interval enddt=0.005# timestept_points=arange(a,b,dt)# array of times# create a list of angles and angular speeds to fill then plottheta_points_app=[]omega_points_app=[]fortint_points:next_angle,next_omega=simpPendApprox(theta_initial,t)theta_points_app.append(next_angle)omega_points_app.append(next_omega)# convert angles back to degrees, and put angular speed in deg/stheta_points_app_deg=[]omega_points_app_deg=[]foriinrange(len(t_points)):theta_points_app_deg.append(theta_points_app[i]*180/pi)omega_points_app_deg.append(omega_points_app[i]*180/pi)################################################################################### Now we solve the exact case ###### define a function to calculate our slopes at a given positiondeff(eta):# what is our angle equal to?theta=eta[0]# what is our angular speed equal to?omega=eta[1]# slope of our anglef_theta=omega# slope of our angular speedf_omega=-omega0**2*sin(theta)# return an array of our slopesreturn(array([f_theta,f_omega],float))# initial conditionseta=array([theta_initial,0.0],float)# initial angular speed of zero# create empty sets to update with values of interest, then invoke Runge-Kuttatheta_points_ex=[]omega_points_ex=[]fortint_points:# add angle of current steptheta_points_ex.append(eta[0])# add angular speed of current stepomega_points_ex.append(eta[1])# calculate where we think we are going from slopes of current positionk1=dt*f(eta)# calculated where we think we are going from slopes at the midpoint of# where we are and where we just calculated we would end upk2=dt*f(eta+0.5*k1)# another mid point calculationk3=dt*f(eta+0.5*k2)# final time through we use the slope of where we think we are goingk4=dt*f(eta+k3)# take a weighted average of all 4 estimations to update our conditionseta+=(k1+2*k2+2*k3+k4)/6# convert angles back to degrees, and put angular speed in deg/stheta_points_ex_deg=[]omega_points_ex_deg=[]foriinrange(len(t_points)):theta_points_ex_deg.append(theta_points_ex[i]*180/pi)omega_points_ex_deg.append(omega_points_ex[i]*180/pi)##### Now to plot/animate ###### begin with some styling options# fontsize of titles and labelsrcParams.update({'font.size':14})# make our axes a little thickerrc('axes',linewidth=1.5)# update the colors of... pretty much everything!withrc_context({'axes.edgecolor':'white','xtick.color':'white', \
'ytick.color':'white','figure.facecolor':'darkslategrey',\
'axes.facecolor':'darkslategrey','axes.labelcolor':'white',\
'axes.titlecolor':'white'}):##### the animations function is not set up to handle a subplot with 3 rows##### so we will make a 2x2 grid with phase plots on the top row with the##### pendulum animations embedded within, then the angle vs time graph below# create figure with an appropriate sizefig=figure(figsize=(10,9))# make our figures less cramped (ironic name for such a function)fig.tight_layout()subplots_adjust(wspace=0.3,hspace=0.3)# create a suplots for our phase diagrams and pendulums## approx in top left ##ax_phase_app=subplot(2,2,1,aspect='equal',adjustable='datalim')title("Angular Speed vs. Angle (Approx)")xlabel("Angle (deg)")ylabel("Angular Speed (deg/s)")axhline(color="k",lw=1)axvline(color="k",lw=1)## exact in top right ##ax_phase_ex=subplot(2,2,2,aspect='equal',adjustable='datalim')title("Angular Speed vs. Angle (Exact)")xlabel("Angle (deg)")ylabel("Angular Speed (deg/s)")axhline(color="k",lw=1)axvline(color="k",lw=1)# create a subplot for our angle vs. time graphs## approx and exact on same plot, two grids wideax_ang=subplot(2,2,(3,4))title("Angle Vs. Time")xlabel("Time (s)")ylabel("Angle (degrees)")axhline(color="k",lw=1)##### Now we animate ###### create an empty list to fill with pictures of our plot each step# set an index to keep track of how many steps we've takenims=[]index=0whileindex<len(t_points):### our pendulums ##### approx ### the rigid rod of length lln_app,=ax_phase_app.plot([0,n*l*sin(theta_points_app[index])],\
[0,-n*l*cos(theta_points_app[index])],color='white',lw=3)# our mass on the endbob_app,=ax_phase_app.plot(n*l*sin(theta_points_app[index]),\
-n*l*cos(theta_points_app[index]),'o',\
markersize=22,color="magenta",zorder=100)## exact ### the rigid rod of length lln_ex,=ax_phase_ex.plot([0,n*l*sin(theta_points_ex[index])],\
[0,-n*l*cos(theta_points_ex[index])],color='white',lw=3)# our mass on the endbob_ex,=ax_phase_ex.plot(n*l*sin(theta_points_ex[index]),\
-n*l*cos(theta_points_ex[index]),'o',\
markersize=22,color="cyan",zorder=100)### phase diagrams ##### approx ##phase_curve_app,=ax_phase_app.plot(theta_points_app_deg[:(index+1)],\
omega_points_app_deg[:(index+1)],\
color="coral",lw=2.5)phase_dot_app,=ax_phase_app.plot(theta_points_app_deg[index:(index+1)],\
omega_points_app_deg[index:(index+1)],\
color="magenta",marker="o",markersize=14)## exact ##phase_curve_ex,=ax_phase_ex.plot(theta_points_ex_deg[:(index+1)],\
omega_points_ex_deg[:(index+1)],\
color="coral",lw=2.5)phase_dot_ex,=ax_phase_ex.plot(theta_points_ex_deg[index:(index+1)],\
omega_points_ex_deg[index:(index+1)],\
color="cyan",marker="o",markersize=14)### angle vs. time ###ang_app,=ax_ang.plot(t_points[:(index+1)],theta_points_app_deg[:(index+1)],\
color="magenta",lw=3)ang_ex,=ax_ang.plot(t_points[:(index+1)],theta_points_ex_deg[:(index+1)],\
color="cyan",lw=3)# add our pictures to the ims1 list for phase diagrams and pendulumsims.append([phase_dot_app,phase_curve_app,phase_dot_ex,phase_curve_ex,\
ln_app,bob_app,ln_ex,bob_ex,ang_app,ang_ex])# only show every 14 time steps (produces smaller gif and speeds things up)index+=14# save animationani=animation.ArtistAnimation(fig,ims,interval=100)writervideo=animation.FFMpegWriter(fps=60)ani.save('../simplePendulum/10_simplePendz.mp4',writer=writervideo)
Simple Pendulum [$\theta_i$ = 10 degrees]
Notice that you can not even tell a difference between the two pendulums, (barring their difference in color). If you look closely on the angle vs time graph, toward the end you can just start to tell that the curves are becoming out of sync. This is fantastic news! It means that our use of the small angle approximation is good for small angles -- at least up to 10 degrees.
I will stop giving you words to read at this point, and leave you to inspect the following graphs and think about what you are seeing. Some things to consider: Pendulums are anharmonic oscillators, which means that as the amplitude increases, so does the period. At small angles, the pendulum does indeed exhibit properties of a harmonic oscillator, meaning its period does not change with increasing amplitude. Pay close attention to how the periods of both the approximated case, and the exact case change with an increase of initial displacement. See how the phase diagrams change. Watch how the pendulums start getting out of sync faster each time you increase the initial displacement. Ask yourself, when does the small angle approximation begin to stop being useful?
