Project 8

Description

In this project we were learning how to use simpson law to perform an estimation and we compaed that with the quad integrator in the scipy integrate package. In this project, we learnt about using numerical techniques to integrate functions that are not integrable. We will use this idea to investigate the period of a simple pendulum. To do so, we collected data from a real pendulum and we will compare them to the exact theory, obtained using simpson's integration technique, to the small angle approximation, and to the very small angle approximation.

Algorithm and Discussion

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:

$$ E = {1 \over 2}m\left(\frac{dx}{dt}\right)^2 + U(x$$

If we solve this for $dx/dt$ we get the result we want: $$ \frac{dx}{dt} = \sqrt{2(E - U(x))/m} $$

or, solving for $dt$ we get:

$$ dt = \sqrt{m \over 2 (E-U(x)) } dx $$

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:

$$ T = 4 \int0^{x{\rm max}} \sqrt{m \over 2 (E-U(x)) } dx $$

Let's apply this idea to the simple pendulum. Use this potential (and use polar coordinates!) to study the simple pendulum:

$$ U(\theta) = -mgl {\rm cos}(\theta) $$

Use this idea to show that for small values of $\theta_m$ the result is the same when ($T = 2 \pi \sqrt{l/g}$)

$$ 4 \sqrt {l/g} \int_0^{\pi/2} {d\phi \over {\sqrt{1 - \sin^2\phi \sin^2{\theta_m \over 2}}}} $$

Code and Implementation

import numpy as np import matplotlib.pyplot as pl import pandas as pd from scipy.optimize import curve_fit # automated curve fitting

df_alex = pd.read_csv('alex1.csv') df_alex.head()

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

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

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

df_dimi1 = pd.read_csv('dimitris1.csv') df_dimi1.head()

pl.plot(df_dimi1.time, df_dimi1.angle, 'b-')

good_time = (df_dimi1.time > 8.5) & (df_dimi1.time<16) time = df_dimi1[good_time].time angle = df_dimi1[good_time].angle pl.plot(time, angle ,'b.') pl.grid()

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.xlabel('time (sec)') pl.ylabel('angle (deg)') pl.legend()

df_dimi2 = pd.read_csv('dimitris2.csv') df_dimi2.head()

pl.plot(df_dimi2.time, df_dimi2.angle, 'b-')

good_time = (df_dimi2.time > 8.5) & (df_dimi2.time<16) time = df_dimi2[good_time].time angle = df_dimi2[good_time].angle pl.plot(time, angle ,'b.') pl.grid()

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

df_mitch1 = pd.read_csv('mitchell1.csv') df_mitch1.head()

pl.plot(df_mitch1.time, df_mitch1.angle, 'b-')

good_time = (df_mitch1.time > 8.5) & (df_mitch1.time<16) time = df_mitch1[good_time].time angle = df_mitch1[good_time].angle pl.plot(time, angle ,'b.') pl.grid()

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

df_soren1 = pd.read_csv('soren1.csv') df_soren1.head()

pl.plot(df_soren1.time, df_soren1.angle, 'b-')

good_time = (df_soren1.time > 8.5) & (df_soren1.time<16) time = df_soren1[good_time].time angle = df_soren1[good_time].angle pl.plot(time, angle ,'b.') pl.grid()

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

df_soren2 = pd.read_csv('soren2.csv') df_soren2.head()

pl.plot(df_soren2.time, df_soren2.angle, 'b-')

bad_time = (df_soren2.time > 8.5) & (df_soren2.time<16) time = df_soren2[bad_time].time angle = df_soren2[bad_time].angle pl.plot(time, angle ,'b.') pl.grid()

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

import os import csv fname = 'rawdata.csv' dataFile = open(os.path.join('rawdata',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)+off # theta 1 corrected for offset th2arr=np.array(theta2_list)+off # 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 def simpson(f, a, b, N, args=None): if args is None: args=[] x=np.linspace(a,b,2*N) # sample the function h=(b-a)/(2*N-1) y=f(*([x] +list(args))) evens = y[2:-2:2] odds = y[1:-1:2] return (y[0] + y[-1] + 2*odds.sum() + 4*evens.sum())*h/3.0 def K_int(phi, m): return 1.0/np.sqrt(1-m*np.sin(phi)**2) def K(m, N=1000): return simpson(K_int, 0, np.pi/2, N=N, args=(m,)) 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)

Results and Conclusion

We fit the angle vs time to a sinusoidal model with amplitude and freqency. We were able to visualize the good time and bad time.

Omega : Angular velocity
Amplitude (A): Measure of the change over a single period.
Offset: It depends on the postion of the coordinate system. It may vary, but it should not affect the outcome of the experiment.

alex1.csv:

Good time
A = 49.6569+/-0.0853 degrees, omega = 3.4957 +/- 0.0008 rad/sec

dimitris1.csv:

Good time
A = 23.9288+/-0.0583 degrees, omega = 3.6381 +/- 0.0011 rad/sec

dimitris2.csv:

Good time
A = 36.5456+/-1.6959 degrees, omega = 3.4562 +/- 0.0209 rad/sec

mitchell1.csv:

Good time
A = 31.3807+/-0.0897 degrees, omega = 3.6081 +/- 0.0021 rad/sec

soren1.csv:

Good time
A = 9.2986+/-0.0468 degrees, omega = 3.6744 +/- 0.0024 rad/sec

soren2.csv:

Bad time
A = 3.3849+/-0.7718 degrees, omega = 3.3386 +/- 0.1057 rad/sec

We have just demonstrated that for large angles, the period of oscillations of a simple pendulum $2\pi(\frac{l}{g})^.5$ fails. In order to study how the pendulum behaves over time, at different angles. Being unable to solve such an equation, we used a numerical integration technique,