Estimating The Period of an Oscillator
======================================

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:

1. 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.
    
2. 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.](https://en.wikipedia.org/wiki/Simpson's_rule).



```python
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](https://docs.scipy.org/doc/scipy-0.18.1/reference/integrate.html) 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.


```python
#
# 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)
```

Project 8, Version A
==========

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:

1. Analyzed the data collected and contributed some data to the full data set
2. Took at least one full cycle of low amplitude data (in the last few seconds of video) to get the right offset for the data you collected and to confirm the low amplitude graph.


Project 8, Version B
==========

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.




# Collecting data with a shaft encoder

Below are examples of pulling data from the digital shaft encoder.



```python
import pandas as pd
from scipy.optimize import curve_fit  # automated curve fitting

```


```python
df_alex = pd.read_csv('2020_PendData/alex1.csv')
df_alex.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>time</th>
      <th>angle</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0.000</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>0.001</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.002</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>3</th>
      <td>0.003</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.003</td>
      <td>0.0</td>
    </tr>
  </tbody>
</table>
</div>




```python
pl.plot(df_alex.time, df_alex.angle, 'b-')
```


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


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

# Repeat for other angles

You need to repeat this process for all the other angles and then produce a plot of $T/T_0$ vs $\theta_{max}$ for:

1. The small angle approximation (horizontal line)
2. The Taylor series expansion (a quadratic)
3. The exact integral (using simpson, or quad)
4. The data you collected

See the graph at the bottom of this notebook for an example of such a plot from our old video analysis data.


# Old Video Analysis Approach

The cells below are here for historical interest. We used to use video analysis to do this experiment. This was how we collected data then.


```python
#
# 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()

```


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

```


```python
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)
```

25a8d1d157cbab1d3f1c9523a8117071298eb0ba
ipynb
Jupyter Notebook
P08-LargeAmplitudePendulum.ipynb
sspickle/sci-comp-notebooks
c7f8afc087cc4013b5d13c698a96e5652a226dc3
MIT
2015-11-08T09:49:18.000Z
2021-08-09T16:07:51.000Z
P08-LargeAmplitudePendulum.ipynb
sspickle/sci-comp-notebooks
c7f8afc087cc4013b5d13c698a96e5652a226dc3
MIT
P08-LargeAmplitudePendulum.ipynb
sspickle/sci-comp-notebooks
c7f8afc087cc4013b5d13c698a96e5652a226dc3
MIT
2015-11-03T18:19:13.000Z
2022-01-26T20:15:30.000Z
Qwen/Qwen-72B
1. YES
2. YES