This week we'll be investigating the period of a pendulum. This is a problem that has been studied for a very long time (e.g., one of earliest workers in the pendulum field was Galeleo himself!). However, most folks (e.g., physics students) just investigate the low amplitude behavior (when the angle between the pendulum and the vertical is pretty small, say 5 or 10 degrees). We are not going to satisfied with such a simple analysis! No, on the contrary, we'll be looking at how the period of the pendulum varies with amplitude. We could use the momentum principle as our startign point for this analsysis, as is usually done, but we're also going to be somewhat unconventional and work instead with the energy principle. Note that this same technique works with any non-linear oscilator (think resonance circuits in electronics, or lasers, or molecules). So, lots to learn! This project has two primary components:
Integrating a function, that is not integrable, numerically so that it's integral can be used as part of the analysis of a problem. As we work out the exact expression for the period of a pendulum we'll encounter such a function and we'll integrate it numerically.
Applying our techniques from last week of fitting functions to data to extract very precise times from a video that we'll collect of the pendulum in motion. This will allow us to compare our analytical approximation and the numerical integral with actual data gathered from a real pendulum.
Let's say we have a system with a potential energy given by some, possibly complex expression $U(x)$. For any system with position $x$ and a potential energy $U$ we have:
\begin{equation} E = {1 \over 2}m\left(\frac{dx}{dt}\right)^2 + U(x). \end{equation}If we solve this for $dx/dt$ we get the result we want: \begin{equation} \frac{dx}{dt} = \sqrt{2(E - U(x))/m}. \end{equation}
or, solving for $dt$ we get:
\begin{equation} dt = \sqrt{m \over 2 (E-U(x)) } dx. \end{equation}If the system is bound we can get the time for one oscillation by integrating the time for one complete period! If the potential is symmetric with respect to $x$ then we only have to do 1/4 of a period and multiply by 4:
\begin{equation} T = 4 \int_0^{x_{\rm max}} \sqrt{m \over 2 (E-U(x)) } dx. \end{equation}Let's apply this idea to the simple pendulum. Use this potential (and use polar coordinates!) to study the simple pendulum:
\begin{equation} U(\theta) = -mgl {\rm cos}(\theta) \end{equation}Use this idea to show that for small values of $\theta_m$ the result is the same as we got in General Physics ($T = 2 \pi \sqrt{l/g}$). You can watch the slides for this lesson see how to carry out the expansion to the next order and find the second order correction to this result. To do the integral it helps to first use: ${\rm cos}\,\theta = 1 - 2 {\rm sin}^2\,(\theta/2)$ and a similar substitution for ${\rm cos}\,\theta_m$. Then change variables to ${\rm sin}\,\phi = {\rm sin}\,(\theta/2)/ {\rm sin}\,(\theta_m/2)$.Then expand in powers of $\theta_m$ At the last step before expansion you should have:
\begin{equation} 4 \sqrt {l/g} \int_0^{\pi/2} {d\phi \over {\sqrt{1 - \sin^2\phi \sin^2{\theta_m \over 2}}}} \end{equation}The integral here is called the "Complete Elliptic Integral of the First Kind". It's famous. It's impossible to do in closed form.
However we can approximate it in several ways. First, we can evalutate the integral numerically (see Simpon's Rule below) or we can estimate the value of the integral by using a Taylor Expansion WRT a small parameter (in our case the amplitude of the pedulum) to get a series that gives us a correction to the simple theory. We'll do both of these things and comapare the results.
Below is a simple python implementation of Simpson's Rule..
In [1]:
import numpy as np
import matplotlib.pyplot as pl
def simpson(f, a, b, N, args=None):
"""
Use Simpson's Rule to estimate an integral
f: function to be integrated
a, b: lower, upper limits
N: "half the number of samples" to take.
args: list of additional arguments (besides 'x') for the function.
"""
if args is None:
args=[]
x=np.linspace(a,b,2*N+1) # sample the function
h=(b-a)/(2*N)
y=f(*([x] +list(args)))
evens = y[2:-2:2]
odds = y[1:-1:2]
return (y[0] + y[-1] + 4*odds.sum() + 2*evens.sum())*h/3.0
def K_int(phi, m):
"""
Integrand for K(m,N) below. This defines the
"Complete Elliptic Integral of the First Kind".
"""
return 1.0/np.sqrt(1-m*np.sin(phi)**2)
def K(m, N=1000):
"""
K(m,N) computes the:
"Complete Elliptic Integral of the First Kind"
with 2*N samples using the simple Simpson's Rule integrator above.
"""
return simpson(K_int, 0, np.pi/2, N=N, args=(m,))
The scipy.integrate package has a number of numerical integrators that you can study. For the purposes of this example we'll use quad which is a simple general purpose integrator. Let's compare that with the Simpson's Rule integrator above.
In [2]:
#
# scipy has an easy "built in" integration routine: integrate.quad
# We'll use this as well for comparison.
#
from scipy.integrate import quad # import a numerical integrator
siList=[] # scipy integrate.quad list
saList=[] # small angle list
srList=[] # Simpson's Rule, N=1000 list
thetArr = np.linspace(0.0, 70, 150) # range of theta values
for theta in thetArr:
thetaRad=theta*np.pi/180.0
m=np.sin(thetaRad/2.0)**2
simpresult = K(m)
quadResult = quad(K_int, 0, np.pi/2, args=(m,))
smallAngle = (np.pi/2)*(1.0+m/4)
srList.append(2*simpresult/np.pi)
siList.append(2*quadResult[0]/np.pi)
saList.append(2*smallAngle/np.pi)
pl.title("Compare integration results (normalized)")
pl.xlabel("$\\theta$")
pl.ylabel("K(m)/K(0)")
pl.plot(thetArr, srList, 'b-', label="Simpson's Rule (N=1000)")
pl.plot(thetArr, siList, 'r-', label="scipy integrate.quad")
pl.plot(thetArr, saList, 'g-', label="Small Angle")
pl.legend(loc=2)
Out[2]:
For project 8 we'll take data from a simple pendulum. This year I plan to collect new data. We have enough students that we can divide and conquer to get a significant data set. If you wish to participate in Version A please let me know. From that data we'll get experimental values for period. I'd like you to compare those experimental results with a) The small angle analytical results and b) the exact theory using Simpson's Rule. You can find below example code that illustrates the idea. In your project you should show how you:
Find any application you like of numerical integration. Use either the Simpson's Rule integrator above, or one of the built-in scipy integrators (e.g., quad) to numerically evaluate the integral. Find some way to validate your code. Either take a simple case where you can compare to an analytical integration, or use "real" data to show that your code is producing something reasonable.
In [3]:
import pandas as pd
from scipy.optimize import curve_fit # automated curve fitting
In [4]:
df_alex = pd.read_csv('2020_PendData/alex1.csv')
df_alex.head()
Out[4]:
In [5]:
pl.plot(df_alex.time, df_alex.angle, 'b-')
Out[5]:
In [6]:
good_time = (df_alex.time > 8.5) & (df_alex.time<16)
time = df_alex[good_time].time
angle = df_alex[good_time].angle
pl.plot(time, angle ,'b.')
pl.grid()
In [7]:
T0 = 2
A0 = 45
def model(t, A, omega, phi, off):
return A*np.sin(omega*t + phi) + off
par, cov = curve_fit(model, time, angle, p0 = (A0, 2*np.pi/T0, 0, 0))
A, omega, phi, off = par
dA, dOm, dPhi, dOff = np.sqrt(np.diag(cov))
print("A = {0:3.4f}+/-{1:3.4f} degrees, omega = {2:3.4f} +/- {3:3.4f} rad/sec".format(A, dA, omega, dOm))
pl.plot(time, angle ,'b.', label='data')
pl.plot(time, model(time, A, omega, phi, off), 'r-', label="fit")
pl.grid()
pl.title("Fit at around 50 degrees")
pl.xlabel('time (sec)')
pl.ylabel('angle (deg)')
pl.legend()
Out[7]:
You need to repeat this process for all the other angles and then produce a plot of $T/T_0$ vs $\theta_{max}$ for:
See the graph at the bottom of this notebook for an example of such a plot from our old video analysis data.
In [15]:
#
# This is an example of getting data out of a video with tracker.
# Note that we can get the angle and the time with the same fit.
# You could use a similar technique with shaft encoder data, or
# you could simply fit large amplitude data with a cos/sin model
# to get the period. Even though the overall fit will not be
# good, the uncertainty in the period would be quite low.
#
import os
import csv
from scipy.optimize import curve_fit # automated curve fitting
def fQuad(t, theta0, alpha, t0):
"""
Use a quadratic to fit the angle vs. time
"""
return theta0 + alpha*(t-t0)**2
fname = 'track270.csv'
dataFile = open(os.path.join('p8-data',fname))
reader = csv.reader(dataFile, delimiter=',')
lineCnt=0
tlist=[]
xlist=[]
ylist=[]
thetaList=[]
for line in reader:
if lineCnt in [0,1]:
pass
else:
try:
tlist.append(float(line[0]))
xlist.append(float(line[1])) # tracker has x-y backwards for some reason.
ylist.append(float(line[2]))
thetaList.append(np.arctan2(ylist[-1],xlist[-1]))
except:
print("Ack!", line)
lineCnt+=1
print("Number of data points:", len(tlist))
tarr=np.array(tlist)
thArr=np.array(thetaList)
sigTheta=0.18/58.0*np.ones(len(tarr)) # ~0.18 cm (1 pixel) error, r~= 58cm
thDeg=thArr*180/np.pi # get theta array in degress for display
theta0e = thArr.sum()/len(thArr)
dt=tarr[1]-tarr[0]
om0 = (thArr[1]-thArr[0])/dt
omf = (thArr[-1]-thArr[-2])/dt
alpha0e = (omf-om0)/(tarr[-1]-tarr[0])
t0e = tarr.sum()/len(tarr)
popt, pcov = curve_fit(fQuad, tarr, thArr, p0=(theta0e, alpha0e, t0e),sigma=sigTheta)
theta0=popt[0]
th0Deg=theta0*180.0/np.pi # theta0 in degrees for display
thErr=np.sqrt(pcov[0,0])*180.0/np.pi # theta0 error in degrees
alpha=popt[1]
t0=popt[2]
t0Err=np.sqrt(pcov[2,2])
print("fit parameters:"+fname)
print("-----------------")
print("theta0:", th0Deg,"+/-",thErr)
print("alpha:", alpha)
print("t0", t0,"+/-",t0Err)
print("K(sin^2(theta0/2)", K(np.sin(theta0/2.0)**2)) # estimate elliptic integral
tstar = np.linspace(tarr[0],tarr[-1],50) # range of times for quadratic model
ystar = fQuad(tstar, theta0, alpha, t0) # value of quadratic model
trange = tarr[-1]-tarr[0] # get range of times
thRange = max(thDeg)-min(thDeg) # range of displayed angles
pl.title("Data: "+fname)
pl.xlabel("time(s)")
pl.ylabel("angle(deg)")
pl.errorbar(tarr, thArr*180/np.pi,fmt='b.',yerr=sigTheta*180.0/np.pi, label="theta")
pl.plot(tstar, ystar*180/np.pi, 'r-', label="quad fit")
pl.plot([tarr[0],tarr[-1]],[th0Deg,th0Deg],'k-') # draw horiz. cross hair
pl.plot([t0,t0],[min(thDeg)-0.3*thRange,max(thDeg)],'k-') # draw vertical cross hair
pl.legend()
Out[15]:
In [25]:
def fPeriodic(t, A, omega, phi, offset):
"""
Use a quadratic to fit the angle vs. time
"""
return A*np.cos(omega*t + phi) + offset
fname = 'track_full_lowamp.csv'
dataFile = open(os.path.join('p8-data',fname))
reader = csv.reader(dataFile, delimiter=',')
lineCnt=0
tlist=[]
xlist=[]
ylist=[]
thetaList=[]
for line in reader:
if lineCnt in [0,1]:
pass
else:
try:
tlist.append(float(line[0]))
xlist.append(float(line[1])) # tracker has x-y backwards for some reason.
ylist.append(float(line[2]))
thetaList.append(np.arctan2(ylist[-1],xlist[-1]))
except:
print("Ack!", line)
lineCnt+=1
print(len(tlist))
tarr=np.array(tlist)
thArr=np.array(thetaList)
sigTheta=0.18/58.0*np.ones(len(tarr)) # ~0.18 cm (1 pixel) error, r~= 58cm
omega0 = 2*np.pi/(120.1-118.6) # estimate omega (from the data)
A0 = (6-(-8))*np.pi/360.0 # estimate Amplitude ((max-min)/2)*(pi/180)
phi0 = 0 # estimate phi (just a guess)
offset0 = 0.0 # estimate offset
print("initial estimates")
print("-----------------")
print("omega0", omega0)
print("A0", A0)
popt, pcov = curve_fit(fPeriodic, tarr, thArr, p0=(A0, omega0, phi0, offset0),sigma=sigTheta)
A=popt[0]
omega=popt[1]
phi=popt[2]
offset=popt[3]
T0 = (2*np.pi/omega)/(1+(A**2)/16) # T0 corrected for amplitude effect
print("fit parameters")
print("-----------------")
print("A:", A*180.0/np.pi)
print("omega:", omega, "corrected period:", T0, "s")
print("phi:", phi*180.0/np.pi)
print("offset:", offset*180.0/np.pi)
print()
tstar = np.linspace(tarr[0],tarr[-1],50)
ystar = fPeriodic(tstar, A, omega, phi, offset)
pl.title("Data: "+fname)
pl.xlabel("time(s)")
pl.ylabel("angle(deg)")
pl.errorbar(tarr, thArr*180/np.pi, fmt='b.',yerr=sigTheta*180.0/np.pi, label="theta")
pl.plot(tstar, ystar*180/np.pi, 'r-', label="quad fit")
pl.legend()
Out[25]:
In [26]:
fname = 'raw-data.csv'
dataFile = open(os.path.join('p8-data',fname))
reader = csv.reader(dataFile, delimiter=',')
lineCnt=0
t1_list=[]
t2_list=[]
theta1_list=[]
theta2_list=[]
for line in reader:
if lineCnt in [0]:
pass
else:
try:
t1_list.append(float(line[1]))
theta1_list.append(float(line[2])*np.pi/180.0) # convert to radians
t2_list.append(float(line[3]))
theta2_list.append(float(line[4])*np.pi/180.0) # convert to radians
except:
print("Ack!", line)
lineCnt+=1
N=len(t1_list)
t1arr=np.array(t1_list)
t2arr=np.array(t2_list)
th1arr=np.array(theta1_list)+offset # theta 1 corrected for offset
th2arr=np.array(theta2_list)+offset # theta 2 corrected for offset
ratio = (t2arr-t1arr)/T0 # array of experimental periods devided by zero amplitude period
thAvg = (th1arr+th2arr)/2.0
vsaa = 1.0*np.ones(N)
saa = 1.0+thAvg*thAvg/16
sraList=[]
for theta in thAvg:
sraList.append(2*K(np.sin(theta/2.0)**2)/np.pi)
thDeg=thAvg*180/np.pi
# array of average amplitudes (in radians)
pl.title("Amplitude of a simple pendulum")
pl.xlabel("$\\theta_{\\rm max}$ (degrees)")
pl.ylabel("T/To")
pl.axis([30,70,0.99,1.1])
pl.plot(thDeg, ratio,'b.', label="experiment")
pl.plot(thDeg, saa, 'r-',label="small angle approx.")
pl.plot(thDeg, vsaa, 'g-', label="very small angle approx.")
pl.plot(thDeg, sraList, 'c-', label="exact theory")
pl.legend(loc=2)
Out[26]: