Kernel: Python 3 (Anaconda 2019)

PHYS2010: Coupled Oscillators

Coupled Oscillator 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 [0]:

#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 file format (a .csv file with some header information). Run the cell below to define it.

NB: This 'readCSV' function differs from the one defined in CO-SingleOscillator.ipynb because it expects two voltage signals, one from each oscillator, to be found in the file.

In [0]:

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
    v2 = dArray[:,2] # this is the voltage of the second oscillator; if you get an error, check that it's in your .csv file. If it's not, you're gonna need to redo it.

    return t, v1, v2

Fit function

Again you will need to perform curve fitting in this Python notebook using the curve_fit function from scipy.optimize. Your initial experimental measurements of amplitude, phase and period or angular frequency will provide the fitting algorithm good "starting points" for determining the optimum parameters. Some you'll have to guess from the values you've recorded.

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

In [0]:

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)

The differential equations describing the behaviour of the coupled oscillators can be 'decoupled' by taking the sum and difference of the two oscillators amplitude functions (see Equation (8) and (9) in the notes). The solutions to Equations (8) and (9) are given in Equation (10) and (11) - both being of the same form.

Step 4: Define a particular function of the form described in Equations (10) and (11) from the notes. The function should take as arguments two amplitudes, A and B, an angular frequency component omega and a constant offset C.

In [0]:

def summedOscillators(t, A, B, omega, C):
    return #PUT A FUNCTION HERE

Step 5: Define another function based on Equation (17) in the notes. To make a function that's good for fitting, it's better to have more levers than less - add in some flexibility around the phase of each cos curve (with phi_1 and phi_2) and allow for a constant offset of the baseline with C.

In [0]:

def beatingOscillators(t, phi0, omega, omega_delta, phi_1, phi_2, C):
    return #PUT A FUNCTION HERE TOO

In-phase coupled oscillators

Step 6: Read in the data file containing your in-phase coupled oscillator trace.

We need to change from Equations (6) and (7) to Equations (8) and (9), i.e. from individual oscillators to sums and differences of oscillators. Do that below by taking the sum and difference of your data, and assigning it to phi_plus and phi_minus respectively.

Plot the sum and the difference against time.

Save this to pdf to print out and stick in your lab book.

In [0]:

#change the filename to read your in-phase oscillator file and complete the rest of the cell
t, v1, v2 = readCSV(r'YOURFILENAMEHERE.csv')

#take the sum and difference
phi_plus =
phi_minus =

In [0]:

#make a plot of phi_plus against time, and phi_minus against time
plt.plot(ADD X AND Y VALUES HERE, label="phi_plus")
plt.plot(ADD X AND Y VALUES HERE, label="phi_minus")

#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
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

Theoretically, perfectly in-phase motion should have $\varphi_-=0$ . Your experimental $\varphi_-$ probably won't be a nice flat 0, but hopefully it's close. Since it's meant to be zero, let's focus on the part that's meant to do something, according to theory.

Step 7: Fit the appropriate function to your $\varphi_+$ using fitCurve and the appropriate function you defined above.

Make starting estimates of the parameters first based on your initial picoscope measurements. What could be good guesses for A_plus and B_plus?

Check your fitted curve against your experimental data when complete to make sure it's a good fit. If it wasn't, try different starting estimates of parameter values, or try a different function to fit it to, and re-run it.

When you are happy with the fit, record your fitted parameters and their uncertainty estimates in your lab book. Calculate a better estimate of $T_+$ from your fitted parameter $\omega_+$ , and its uncertainty.

In [0]:

#set up your parameters
T_plus =      #add your period here
omega_plus = 2*np.pi/T_plus
A_plus =      #estimate A_plus
B_plus =      #estimate B_plus
C = phi_plus.mean()
startPars = (A_plus, B_plus, omega_plus, C)

phi_plusFit, fitPars, deltaPars = fitCurve(summedOscillators, t, phi_plus, startPars)

In [0]:

plt.plot(ADD X AND Y VALUES HERE, label="expt")
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
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

In [0]:

T_plusFit = 2*np.pi/fitPars[2]
deltaT_plusFit = 2*np.pi*deltaPars[2]/fitPars[2]**2

print(T_plusFit, deltaT_plusFit)

Out-of-phase coupled oscillators.

Step 8: Read in the data file containing your out-of-phase coupled oscillator trace.

Like above, we need to change from Equations (6) and (7) to Equations (8) and (9), i.e. from individual oscillators to sums and differences of oscillators. Do that below by taking the sum and difference of your data again, and assigning it to phi_plus and phi_minus respectively.

Plot the sum and the difference against time.

Save this to pdf to print out and stick in your lab book.

In [0]:

#change the filename to read your out-of-phase oscillator file and complete the rest of the cell
t, v1, v2 = readCSV(r'YOURFILENAMEHERE.csv')

#take the sum and difference
phi_plus =
phi_minus =

In [0]:

#make a plot of phi_plus against time, and phi_minus against time
plt.plot(ADD X AND Y VALUES HERE, label="phi_plus")
plt.plot(ADD X AND Y VALUES HERE, label="phi_minus")

#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
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

Theoretically, perfectly out-of-phase motion has $\varphi_+=0$ - that's the opposite to above, just as you'd expect. Again, it probably won't be perfectly zero - good to discuss, but won't help us get more analysis done. Let's focus on $\varphi_-$ .

Step 9: Fit the appropriate function to your $\varphi_-$ using the fitCurve function and the function you defined above. Make starting estimates of the parameters first based on your initial picoscope measurements as you did above.

When you are happy with the fit, record your fitted parameters and their uncertainty estimates in your lab book.

In [0]:

#set up your parameters
T_minus =    #add your period here
omega_minus = 2*np.pi/T_minus
A_minus =    #estimate A_minus
B_minus =    #estimate B_minus
C = phi_minus.mean()
startPars = (A_minus, B_minus, omega_minus, C)

phi_minusFit, fitPars, deltaPars = fitCurve(summedOscillators, t, phi_minus, startPars)

In [0]:

plt.plot(ADD X AND Y VALUES HERE, label="expt")
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
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

In [0]:

T_minusFit = 2*np.pi/fitPars[2]
deltaT_minusFit = 2*np.pi*deltaPars[2]/fitPars[2]**2

print(T_minusFit, deltaT_minusFit)

Maximum beat amplitude

Step 10: Read in the data file containing your maximum beat coupled oscillator trace. Plot the trace of both oscillators against time. Save this to pdf to print out and stick in your lab book.

In [0]:

#change the filename to read your maximum beat amplitude oscillator file and complete the rest of the cell
t, v1, v2 = readCSV(r'YOURFILENAMEHERE.csv')

#make a plot of phi_plus against time, and phi_minus against time
plt.plot(ADD X AND Y VALUES HERE, label="oscillator 1")
plt.plot(ADD X AND Y VALUES HERE, label="oscillator 2")

#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
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

Theoretically, one of the oscillators has the same trace as the other oscillator but shifted by $\pi$ , as indicated by Equation (17) and (18) in the notes.

Step 11: Fit the appropriate function to both oscillators to get fitted estimates of the parameters $\omega_\Delta$ , $\omega$ and therefore the beat period $T_\Delta$ and period $T$ . Use your picoscope measurements as starting values for the fitted parameters. Plot your fitted results to make sure they're sensible. Finally record your fitted parameter values and their uncertainties in your lab book.

Oscillator 1

In [0]:

#initial estimates from picoscope measurements for oscillator 1
T_Delta = 
omega_Delta = 2*np.pi/T_Delta
T = 
omega= 2*np.pi/T
phi0 = 
phi1 = 1
phi2 = 1
C = v1.mean()
startPars=(phi0, omega, omega_Delta, phi1, phi2, C)

#do the fitting thing
v1Fit, v1FitPars, v1DeltaPars = fitCurve(beatingOscillators, t, v1, startPars)

In [0]:

plt.plot(ADD X AND Y VALUES HERE, label="expt")
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
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

In [0]:

#calculate a fitted value of T_Delta and T and their uncertainties from the fitted parameters
T_DeltaFit1 = 2*np.pi/v1FitPars[2]
print("T_delta is", T_DeltaFit1)
deltaT_DeltaFit1 = 2*np.pi*v1DeltaPars[2]/v1FitPars[2]**2
print("its uncertainty is", deltaT_DeltaFit1)

TFit1= 2*np.pi/v1FitPars[1]
print("T is", TFit1)
deltaTFit1 = 2*np.pi*v1DeltaPars[1]/v1FitPars[1]**2
print("its uncertainty is", deltaTFit1)

Oscillator 2

In [0]:

#initial estimates from picoscope measurements for oscillator 1
T_Delta = 
omega_Delta = 2*np.pi/T_Delta
T = 
omega= 2*np.pi/T
phi0 = 
phi1 = 1
phi2 = 1
C = v2.mean()
startPars=(phi0, omega, omega_Delta, phi1, phi2, C)

#do the fitting thing
v2Fit, v2FitPars, v2DeltaPars = fitCurve(beatingOscillators, t, v2, startPars)

In [0]:

plt.plot(ADD X AND Y VALUES HERE, label="expt")
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
plt.savefig("YOURFILENAMEHERE.pdf", dpi=300, orientation="landscape")

In [0]:

#calculate a fitted value of T_Delta and T and their uncertainties from the fitted parameters
T_DeltaFit2 = 2*np.pi/v2FitPars[2]
print("T_delta is", T_DeltaFit1)
deltaT_DeltaFit2 = 2*np.pi*v2DeltaPars[2]/v2FitPars[2]**2
print("its uncertainty is", deltaT_DeltaFit2)

TFit2= 2*np.pi/v2FitPars[1]
print("T is", TFit2)
deltaTFit2 = 2*np.pi*v2DeltaPars[1]/v2FitPars[1]**2
print("its uncertainty is", deltaTFit2)

Step 12: If you couldn't manage to fit a curve to both oscillators equally well, use the fitted values of $T_\Delta$ and $T$ , $\omega_\Delta$ and $\omega$ , that come from the more successful fit when proceeding with the next calculations mentioned in your notes.