CoCalc Public Filespython-labs / coupled_oscillator / CO-SingleOscillator.ipynbOpen in with one click!
Authors: Eric Howard, Adam Joyce, Alessandro Maini, Jason Twamley
Views : 20
License: GNU General Public License v3.0

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

Read function

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 f.readline() 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 TT as a function of LL 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 t, v1 = readCSV(r'YOURFILENAMEHERE.csv') #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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend - this can stay empty plt.legend() #add title 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 t, v1 = readCSV(r'YOURFILENAMEHERE.csv') #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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend - this can stay empty plt.legend() #add title 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 t, v1 = readCSV(r'YOURFILENAMEHERE.csv') #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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend - this can stay empty plt.legend() #add title 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 t, v1 = readCSV(r'YOURFILENAMEHERE.csv') #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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend - this can stay empty plt.legend() #add title 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 t, v1 = readCSV(r'YOURFILENAMEHERE.csv') #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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend - this can stay empty plt.legend() #add title 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 t, v1 = readCSV(r'YOURFILENAMEHERE.csv') #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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend - this can stay empty plt.legend() #add title 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 t, v1 = readCSV(r'YOURFILENAMEHERE.csv') #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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend - this can stay empty plt.legend() #add title 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 t, v1 = readCSV(r'YOURFILENAMEHERE.csv') #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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend - this can stay empty plt.legend() #add title 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 T2T^2 being proportional to LL. 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 T2T^2 as a function of LL. 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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend plt.legend() #add title 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") #add x label plt.xlabel("X LABEL HERE") #add y label plt.ylabel("Y LABEL HERE") #add legend plt.legend() #add title 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.