CoCalc Public Filespython-labs / coupled_oscillator / CO-SingleOscillator.ipynb
Views : 20

# PHYS2010: Coupled Oscillators

## Single Pendulum Analysis

NB: Read the instructions given before each grey code cell, edit what you need to, and run it. Make sure you run the cells in order so there's a higher chance they'll work!

When you're done, save/screenshot/print any graphs you've made and any fitted parameters that you find, and include them in your lab notebook/report.

Step 1: We need scientific packages to do the calculations and fitting that we want to do, so let's import Numpy, Pyplot from Matplotlib, and curve_fit from the relevant python packages.

In [ ]:
#you can just run this cell
import numpy as np
import matplotlib.pyplot as plt

from scipy.optimize import curve_fit


Step 2: To be able to import your experimental data from the file from Picoscope, you will need to define the function to read the pico file format (a .csv file with some header information). Run the cell below to define it.

In [ ]:
def readCSV(fileName):

with  open(fileName,'r') as f:
for k in range(3): #The first 3 lines of the file are all header information that we don't need, so we read them in and ignore them

data = []

for line in f:
data.append([float(x) for x in line.split(",")]) #The data are (or should be) floats, and they're separated by commas

dArray = np.array(data)

t = dArray[:,0] # this is the time
v1 = dArray[:,1] # this is the voltage of the first oscillator

return t, v1



## Fit function

Fitting a function to some data is a delicate art form. You need to understand what you're trying to do and what the algorithm you're using does. In this case we are using the curve_fit function from scipy.optimize. Take a look at the documentation of the function by googling curve_fit python - googling functions and reading their documentation is a vital skill for a programmer.

In a nutshell, any fitting algorithm minimises the difference between the data points and a given function. The function you use must have a set of free parameters that the algorithm changes in order to minimise the differences.

You typically have to provide the algorithm with a "starting point" for the parameters, so that the algorithm knows roughly where the expected optimum parameters are. Minimising nonlinear functions in multidimensional spaces is very tricky and you can easily fall in a "local minimum", that is a set of parameters which are not optimal, but that are better than any parameters in their neighbourhood.

Step 3: Run the next cell to define a function that fits a curve/function (generically called 'myFunc') to a data set of x's and y's, given a set of starting parameter values.

You can just run this.

In [ ]:
def fitCurve(myFunc, x, y, startPars):
optimisedPars, covarPars = curve_fit(myFunc, x, y, startPars)
yFit = myFunc(x, *optimisedPars)
uncertPars=np.sqrt(np.diag(covarPars))
return yFit, optimisedPars, uncertPars
#This function returns the optimised parameters (optimisedPars) and the 'data' resulting from using the
#chosen function with the optimal parameters (yFit) - this is the fitted curve you want to plot.
#Also returns an estimate of the uncertainties of the optimised parameters (uncertPars)


Step 4: Define a particular function based loosely on Equation (2a) in the lab notes that describes a single pendulum's displacement as a function of time. The function should take as arguments: an array of time values, t, amplitude A, angular frequency omega, phase shift phi and a non-zero offset C.

In [ ]:
#complete the function
def harmonicOscillator(t, A, omega, phi, C):
return #PUT A FUNCTION HERE


You should have gathered displacement data for an oscillating, single pendulum for at least 8 different lengths. For each length, there is a different .csv file containing the data.

You are going to read in the pico file for a given length and fit the harmonicOscillator function to the experimental data, using your picoscope measurements of amplitude, phase and period as starting estimates for the fit. The fitted angular frequency parameter will a more accurate estimate of the period of the pendulum at that particular length.

Once you have repeated this for all the lengths of the single pendulum, you can plot a graph of $T$ as a function of $L$ and compare it to the relationship expected for an ideal, simple pendulum (Equation (1) in the notes).

Step 5: Load in your first .csv file, set the starting parameters for the curve fit from your experimental results, perform the curve fit and obtain a better estimate of the angular frequency and its uncertainty.

In [ ]:
# enter the correct directory and filename for a single pendulum result file - easiest if it's in the same folder as this file

#plot experiment results to make sure everything's hunky-dory
plt.plot(t, v1)

#set starting parameters for curve fit from your picoscope measurements: A, omega, phi, C
A = 0.5            #change this
T1exp = 1        #change this
deltaT1exp = 0.005 #change this
phi = 3.14         #change this
omega = 2*np.pi/T1exp
C = v1.mean()   #starting value for your offset can just be calculated from your data like this

L1 =  #Put your length in metres here
deltaL1 =  #Put your length uncertainty in metres here

startingPars = (A, omega, phi, C)

#fit the function from Equation (2a) to your data using fitCurve
v1Fit, fitPars, deltaPars = fitCurve(harmonicOscillator, t, v1, startingPars)



Step 6: Next plot the data and the optimal curve found by the algorithm together to see how good the fit is. If you're satisfied with the fit continue to Step 7!

In [ ]:
#plot experimental data and label it
plt.plot(ADD X AND Y VALUES HERE, label="expt")

#plot fitted data and label it
plt.plot(ADD X AND Y VALUES HERE, label="fit")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

#add legend - this can stay empty
plt.legend()
plt.suptitle("THIS NEEDS A TITLE")

#add an appropriate file name for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")


Step 7: Save your new fitted estimate of period and its uncertainty so you can use them later.

In [ ]:
#complete the following
T1fit = 2*np.pi/fitPars[1]
deltaT1fit = deltaPars[1]*2*np.pi/fitPars[1]**2



Step 8: Repeat Steps 5-7 for every length of the single pendulum you tested until you have a complete set of fitted periods, their associated uncertainties and all the corresponding pendulum lengths and their uncertainties.

## Length 2

In [ ]:
# enter the correct directory and filename for a single pendulum result file - easiest if it's in the same folder as this file

#plot experiment results to make sure everything's hunky-dory
plt.plot(t, v1)

#set starting parameters for curve fit from your picoscope measurements: A, omega, phi, C
A = 0.5            #change this
T2exp = 1        #change this
deltaT2exp = 0.005 #change this
phi = 3.14         #change this
omega = 2*np.pi/T2exp
C = v1.mean()   #starting value for your offset can just be calculated from your data like this

L2 =  #Put your length in metres here
deltaL2 =  #Put your length uncertainty in metres here

startingPars = (A, omega, phi, C)

#fit the function from Equation (2a) to your data using fitCurve
v1Fit, fitPars, deltaPars = fitCurve(harmonicOscillator, t, v1, startingPars)

#plot experimental data and label it
plt.plot(ADD X AND Y VALUES HERE, label="expt")

#plot fitted data and label it
plt.plot(ADD X AND Y VALUES HERE, label="fit")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

#add legend - this can stay empty
plt.legend()
plt.suptitle("THIS NEEDS A TITLE")

#add appropriate filename for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

#complete the following
T2fit = 2*np.pi/fitPars[1]
deltaT2fit = deltaPars[1]*2*np.pi/fitPars[1]**2



## Length 3

In [ ]:
# enter the correct directory and filename for a single pendulum result file - easiest if it's in the same folder as this file

#plot experiment results to make sure everything's hunky-dory
plt.plot(t, v1)

#set starting parameters for curve fit from your picoscope measurements: A, omega, phi, C
A = 0.5            #change this
T3exp = 1        #change this
deltaT3exp = 0.005 #change this
phi = 3.14         #change this
omega = 2*np.pi/T3exp
C = v1.mean()   #starting value for your offset can just be calculated from your data like this

L3 =  #Put your length in metres here
deltaL3 =  #Put your length uncertainty in metres here

startingPars = (A, omega, phi, C)

#fit the function from Equation (2a) to your data using fitCurve
v1Fit, fitPars, deltaPars = fitCurve(harmonicOscillator, t, v1, startingPars)

#plot experimental data and label it
plt.plot(ADD X AND Y VALUES HERE, label="expt")

#plot fitted data and label it
plt.plot(ADD X AND Y VALUES HERE, label="fit")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

#add legend - this can stay empty
plt.legend()
plt.suptitle("THIS NEEDS A TITLE")

#add appropriate filename for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

#complete the following
T3fit = 2*np.pi/fitPars[1]
deltaT3fit = deltaPars[1]*2*np.pi/fitPars[1]**2



## Length 4

In [ ]:
# enter the correct directory and filename for a single pendulum result file - easiest if it's in the same folder as this file

#plot experiment results to make sure everything's hunky-dory
plt.plot(t, v1)

#set starting parameters for curve fit from your picoscope measurements: A, omega, phi, C
A = 0.5            #change this
T4exp = 1        #change this
deltaT4exp = 0.005 #change this
phi = 3.14         #change this
omega = 2*np.pi/T4exp
C = v1.mean()   #starting value for your offset can just be calculated from your data like this

L4 =  #Put your length in metres here
deltaL4 =  #Put your length uncertainty in metres here

startingPars = (A, omega, phi, C)

#fit the function from Equation (2a) to your data using fitCurve
v1Fit, fitPars, deltaPars = fitCurve(harmonicOscillator, t, v1, startingPars)

#plot experimental data and label it
plt.plot(ADD X AND Y VALUES HERE, label="expt")

#plot fitted data and label it
plt.plot(ADD X AND Y VALUES HERE, label="fit")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

#add legend - this can stay empty
plt.legend()
plt.suptitle("THIS NEEDS A TITLE")

#add appropriate filename for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

#complete the following
T4fit = 2*np.pi/fitPars[1]
deltaT4fit = deltaPars[1]*2*np.pi/fitPars[1]**2



## Length 5

In [ ]:
# enter the correct directory and filename for a single pendulum result file - easiest if it's in the same folder as this file

#plot experiment results to make sure everything's hunky-dory
plt.plot(t, v1)

#set starting parameters for curve fit from your picoscope measurements: A, omega, phi, C
A = 0.5            #change this
T5exp = 1        #change this
deltaT5exp = 0.005 #change this
phi = 3.14         #change this
omega = 2*np.pi/T5exp
C = v1.mean()   #starting value for your offset can just be calculated from your data like this

L5 =  #Put your length in metres here
deltaL5 =  #Put your length uncertainty in metres here

startingPars = (A, omega, phi, C)

#fit the function from Equation (2a) to your data using fitCurve
v1Fit, fitPars, deltaPars = fitCurve(harmonicOscillator, t, v1, startingPars)

#plot experimental data and label it
plt.plot(ADD X AND Y VALUES HERE, label="expt")

#plot fitted data and label it
plt.plot(ADD X AND Y VALUES HERE, label="fit")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

#add legend - this can stay empty
plt.legend()
plt.suptitle("THIS NEEDS A TITLE")

#add appropriate filename for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

#complete the following
T5fit = 2*np.pi/fitPars[1]
deltaT5fit = deltaPars[1]*2*np.pi/fitPars[1]**2



## Length 6

In [ ]:
# enter the correct directory and filename for a single pendulum result file - easiest if it's in the same folder as this file

#plot experiment results to make sure everything's hunky-dory
plt.plot(t, v1)

#set starting parameters for curve fit from your picoscope measurements: A, omega, phi, C
A = 0.5            #change this
T6exp = 1        #change this
deltaT6exp = 0.005 #change this
phi = 3.14         #change this
omega = 2*np.pi/T6exp
C = v1.mean()   #starting value for your offset can just be calculated from your data like this

L6 =  #Put your length in metres here
deltaL6 =  #Put your length uncertainty in metres here

startingPars = (A, omega, phi, C)

#fit the function from Equation (2a) to your data using fitCurve
v1Fit, fitPars, deltaPars = fitCurve(harmonicOscillator, t, v1, startingPars)

#plot experimental data and label it
plt.plot(ADD X AND Y VALUES HERE, label="expt")

#plot fitted data and label it
plt.plot(ADD X AND Y VALUES HERE, label="fit")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

#add legend - this can stay empty
plt.legend()
plt.suptitle("THIS NEEDS A TITLE")

#add appropriate filename for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

#complete the following
T6fit = 2*np.pi/fitPars[1]
deltaT6fit = deltaPars[1]*2*np.pi/fitPars[1]**2



## Length 7

In [ ]:
# enter the correct directory and filename for a single pendulum result file - easiest if it's in the same folder as this file

#plot experiment results to make sure everything's hunky-dory
plt.plot(t, v1)

#set starting parameters for curve fit from your picoscope measurements: A, omega, phi, C
A = 0.5            #change this
T7exp = 1        #change this
deltaT7exp = 0.005 #change this
phi = 3.14         #change this
omega = 2*np.pi/T7exp
C = v1.mean()   #starting value for your offset can just be calculated from your data like this

L7 =  #Put your length in metres here
deltaL7 =  #Put your length uncertainty in metres here

startingPars = (A, omega, phi, C)

#fit the function from Equation (2a) to your data using fitCurve
v1Fit, fitPars, deltaPars = fitCurve(harmonicOscillator, t, v1, startingPars)

#plot experimental data and label it
plt.plot(ADD X AND Y VALUES HERE, label="expt")

#plot fitted data and label it
plt.plot(ADD X AND Y VALUES HERE, label="fit")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

#add legend - this can stay empty
plt.legend()
plt.suptitle("THIS NEEDS A TITLE")

#add appropriate filename for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

#complete the following
T7fit = 2*np.pi/fitPars[1]
deltaT7fit = deltaPars[1]*2*np.pi/fitPars[1]**2



## Length 8

In [ ]:
# enter the correct directory and filename for a single pendulum result file - easiest if it's in the same folder as this file

#plot experiment results to make sure everything's hunky-dory
plt.plot(t, v1)

#set starting parameters for curve fit from your picoscope measurements: A, omega, phi, C
A = 0.5            #change this
T8exp = 1        #change this
deltaT8exp = 0.005 #change this
phi = 3.14         #change this
omega = 2*np.pi/T8exp
C = v1.mean()   #starting value for your offset can just be calculated from your data like this

L8 =  #Put your length in metres here
deltaL8 =  #Put your length uncertainty in metres here

startingPars = (A, omega, phi, C)

#fit the function from Equation (2a) to your data using fitCurve
v1Fit, fitPars, deltaPars = fitCurve(harmonicOscillator, t, v1, startingPars)

#plot experimental data and label it
plt.plot(ADD X AND Y VALUES HERE, label="expt")

#plot fitted data and label it
plt.plot(ADD X AND Y VALUES HERE, label="fit")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

#add legend - this can stay empty
plt.legend()
plt.suptitle("THIS NEEDS A TITLE")

#add appropriate filename for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

#complete the following
T8fit = 2*np.pi/fitPars[1]
deltaT8fit = deltaPars[1]*2*np.pi/fitPars[1]**2



## Testing the relationship between period and length

Now that you've plotted and fitted all those curves, let's group together all the periods (both experimental and fitted) and the lengths, and see if we can find a relationship between them.

In [ ]:
#so long as you've entered your lengths and periods as they are above, this should just work. If you've added any more points, you'll need to change the names of those variables and add them in below

expPeriods = np.array([T1exp, T2exp, T3exp, T4exp, T5exp, T6exp, T7exp, T8exp])
expDeltaPeriods = np.array([deltaT1exp, deltaT2exp, deltaT3exp, deltaT4exp, deltaT5exp, deltaT6exp, deltaT7exp, deltaT8exp])

fitPeriods = np.array([T1fit, T2fit, T3fit, T4fit, T5fit, T6fit, T7fit, T8fit])
fitDeltaPeriods = np.array([deltaT1fit, deltaT2fit, deltaT3fit, deltaT4fit, deltaT5fit, deltaT6fit, deltaT7fit, deltaT8fit])

lengths = np.array([L1, L2, L3, L4, L5, L6, L7, L8])
deltaLengths = np.array([deltaL1, deltaL2, deltaL3, deltaL4, deltaL5, deltaL6, deltaL7, deltaL8])



Theoretically, from Equation (1), the period of oscillation of a simple pendulum is proportional to the square root of the pendulum's length. To fit this using polyfit and polyval, it's easier to change it to $T^2$ being proportional to $L$. According to Equation (1), this should be a linear relationship.

Step 9: Use Numpy's polyfit and polyval to find the best fitting parameters to your T and L experimental data, given Equation (1). Do the same for the fitted period values of $T^2$ as a function of $L$. Make sure to estimate the uncertainty in the fitted parameters by setting the cov flag to "True". Plot the resulting fit. Is a linear trendline required for a good fit, or some other polynomial?

In [ ]:
#polyfit generates the co-efficients (and covariance) for a polynomial fit to the input. The output coefficients list from highest power to lowest. Try for different degrees (starting with 1, linear):
expFitPars, expCovarPars = np.polyfit(lengths, expPeriods**2, 1, cov=True)
fittedFitPars, fittedCovarPars = np.polyfit(lengths, fitPeriods**2, 1, cov=True)

#polyval then uses these co-efficients to generate a fitted curve:
expPeriodSquaredFit = np.polyval(expFitPars, lengths)
fittedPeriodSquaredFit = np.polyval(fittedFitPars, lengths)

#we can work out the uncertainties from the covariance spat out by polyfit:
expPeriodSquaredDelta = np.sqrt(np.diag(expCovarPars))
fittedPeriodSquaredDelta = np.sqrt(np.diag(fittedCovarPars))


In [ ]:
#let's plot the fit generated from the EXPERIMENTAL period values
plt.plot(lengths, expPeriodSquaredFit, label="fit")
plt.plot(lengths, expPeriods**2, label="expt")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

plt.legend()
plt.suptitle("NEED A TITLE HERE")

#add appropriate filename for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

In [ ]:
#let's plot the fit generated from the FITTED period values
plt.plot(lengths, fittedPeriodSquaredFit, label="fit")
plt.plot(lengths, fitPeriods**2, label="expt")

plt.xlabel("X LABEL HERE")
plt.ylabel("Y LABEL HERE")

plt.legend()
plt.suptitle("NEED A TITLE HERE")

#add appropriate filename for pdf of graph - this will save a pdf version in the same folder as this notebook
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")


Step 10: Print the fitted parameters and their uncertainties and record them in your lab book.

In [ ]:
print("For the experimental values, the fit is y = "+str(expFitPars[0])+"x + "+str(expFitPars[1])+".", "The uncertainty in the gradient is "+str(expPeriodSquaredDelta[0])+" and the uncertainty in the constant is "+str(expPeriodSquaredDelta[1]))

print("For the fitted values, the fit is y = "+str(fittedFitPars[0])+"x + "+str(fittedFitPars[1])+".", "The uncertainty in the gradient is "+str(fittedPeriodSquaredDelta[0])+" and the uncertainty in the constant is "+str(fittedPeriodSquaredDelta[1]))


Is this what you expected? Continue with the questions in your lab notes and discuss. If you have time, you could try modifying the above for a different degree polynomial fit (i.e. quadratic rather than linear) if the linear fit doesn't look that great.