Springs with Forced Displacement#
Basic Imports:#
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.animation as animation
from IPython.display import HTML
from scipy.signal import find_peaks
from typing import NamedTuple
from scipy.fft import fft, fftfreq, ifft, fftshift
from scipy import signal
Animation Functions:#
def draw_spring(x0, y0, x1, y1, coils=10, amplitude=0.1, linewidth=2, color='blue'):
"""Draws a spring between (x0, y0) and (x1, y1)."""
num_points=500
dx = x1 - x0
dy = y1 - y0
length = np.sqrt(dx**2 + dy**2)
angle = np.arctan2(dy, dx)
x_line = np.linspace(x0, x1, num_points)
y_line = np.linspace(y0, y1, num_points)
coil_x = amplitude * np.cos(np.linspace(0, coils * 2 * np.pi, num_points))
coil_y = amplitude * np.sin(np.linspace(0, coils * 2 * np.pi, num_points))
x_spring = x_line + coil_x * np.cos(angle + np.pi/2)
y_spring = y_line + coil_y * np.sin(angle + np.pi/2)
return x_spring, y_spring # Return the spring coordinates
# Animation setup
def SetupAnimation():
fig, ax = plt.subplots()
ax.plot([0.0,0.0],[-0.25,0.25],color='red',linewidth=4)
line, = ax.plot([], [], lw=2, color='blue') # Initialize an empty line object
line2, = ax.plot([], [], lw=2, color='red') # Initialize an empty line object
center_x=0
center_y=0
radius = 0.1
circle = patches.Circle((center_x, center_y), radius, color='blue', fill=True)
ax.add_artist(circle)
ax.set_xlim(-6, 45) # Set appropriate limits for x
ax.set_ylim(-1, 1) # Set appropriate limits for y. Adjust as needed.
#ax.set_aspect('equal') # Important for proper spring visualization
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Spring Animation')
return fig, ax, line, line2, circle
def animate(i, x, Y_t, line,line2,circle, L=20):
x1 = x[i] # x-coordinate changes over time
xs, ys = draw_spring(0, 0, x1 + L, 0.0, coils=8, amplitude=0.10, color='blue')
line.set_data(xs, ys)
xs, ys = draw_spring(x1 + L, 0, Y_t[i] + 2 * L, 0.0, coils=8, amplitude=0.10, color='red')
line2.set_data(xs, ys)
new_x = x1
circle.center = (new_x + L, 0)
return line, line2, circle
Exercise 0. Background#
In this assignment, we continue working with two coupled springs with damping and forced displacements. Recall that we store our spring parameters as
class SpringParams(NamedTuple):
k1: float
m: float
k2: float
c: float = 0.0
where we will consistently be using the parameters
params= SpringParams(k1=1.0, m=1.0, k2=1.0, c=0.5)
Our right spring is pulled by a Harmonic displacement of the form
and again we define a Harmonic displacement function and store our Harmonic Displacement parameters as
class HarmonicDisplaceParams(NamedTuple):
Y0: float
omega: float
phi: float = 0.0
def GetDisplacement_harmonic(D,t):
y = D.Y0 * np.cos(D.omega * t+D.phi)
return y
We are considering the scenario where the two springs have springs constants \(k_1\) (left) and \(k_2\) (right) connected to a mass \(m\) between them. The left spring is fixed to a wall and the right spring is displaced by a function \(y(t)\). The equation of motion for such a setup is
As in the “Classical Harmonic Oscillator” assignment, turn our second order differential equation into two first order differential equations writing
def equations(t, state, springParams, D,Displacement=GetDisplacement_harmonic):
x, v = state
### Do stuff
return [xdot,vdot]
You may just copy this from your previous assignment if you wish.
###ANSWER HERE
Exercise 1. Linear Combination of Periodic Forcing#
a. Two Harmonic Forces#
In this exercise, we’d like to understand what happens if we apply a driving force which is a sum of multiple Harmonic forces - i.e.
where
\(Y_1(t) = 4 \cos (2\pi (0.25) t)\)
\(Y_2(t) = 8 \cos (2\pi (0.1) t + 1.0)\)
with an initial condition of \(x=1.0; v=1.0\).
To solve for \(x(t)\), we can separately solve for \(x_1(t)\) (with Harmonic force \(Y_1(t)\)) and \(x_2(t)\) (with Harmonic force \(Y_2(t)\)) and add them together to get \(x(t)=x_1(t) +x_2(t)\). When we do this, we need to be careful about our initial conditions: the sum of the initial conditions for \(x_1(t)\) and \(x_2(t)\) must be equal to the initial conditions for \(x(t)\). Let our initial conditions be \([0.3,0.6]\) and \([0.7,0.4]\) respectively.
Using t= np.linspace(0, 20, 10000), solve for and plot \(x_1(t)\), \(x_2(t)\) and \(x(t)\)
###ANSWER HERE
b. Fourier Analysis#
In the previous section, we saw how we could solve our problem with a forced displacement of
when \(Y_1(t)\) and \(Y_2(t)\) are both Harmonic displacements that we knew. Suppose instead that we only knew about \(Y(t)\) (plot it below).
###ANSWER HERE
Then to solve for our differential equation we would need to figure out how to write \(Y(t)\) as a linear combination of Harmonic displacements (which turns out to be possible for any curve) decomposing is as
where we solve for \(\omega_i\) and \(\phi_i\) for every \(i\).
We can do this using a FFT. Defining Y_t as the command curve we can do
fft_freqs=fftfreq(len(y_t),t[1]-t[0])
the_fft=fft(D_t)
Then plot the absolute value of \(F(t)\) and this should tell us what frequencies (as a fraction of \(\pi\)) make up our curves
plt.plot(fft_freqs,np.abs(the_fft)/len(t),'.')
plt.xlim(-1.2,1.2)
plt.xlabel('Frequency $\omega$')
plt.ylabel('$|FFT(D(t))|$')
plt.show()
cutoff = 0.1
mask=np.where(np.abs(the_fft)>cutoff)
print("Frequencies: ",fft_freqs[mask])
Run this and you should see (both from the graph and printing) that we get back \(\pm 0.25\) and \( \pm 0.1\) as large frequencies.
As you can read in the documentation, scipy.fft does not take into account our dimensions of time. Because of that, to get the right amplitudes \(Y_i\), we need to normalize the_fft by len(t).
Also, due to floating-point errors, nothing in practice is exactly 0. For this reason we had to introduce a cutoff.
###ANSWER HERE
To get the phases of our curves, we should then plot
plt.plot(fft_freqs[mask],np.angle(the_fft[mask]),'.')
plt.xlim(-1,1)
plt.ylim(-np.pi,np.pi)
plt.xlabel('Frequency $\omega$')
plt.ylabel('Phase $\phi$')
plt.show()
###ANSWER HERE
From the graphs above we can extract that our curve was
This is what was expected given what we put in.
Exercise 2. Linear Combination of Delta Function Forcing#
In the previous problem we decomposed our signal into a sum of Harmonic signals. In this problem we are going to decompose our solution into a sum of (almost) delta functions.
a. Delta Function#
Our delta function will be parameterized as an input of height \(Y_0\) for times \(t_1 <t < t_2\) (where we will choose \(t_1\) and \(t_2\) to be close together such that \(t_2-t_1=dt\).
To use this type of forced displacement, we are going to define a new class to store the parameters for this delta function
class DeltaDisplaceParams(NamedTuple):
Y0: float
t1: float
t2: float
In addition, write a function
def GetDisplacement_delta(D,t):
### do stuff
return y
which returns \(Y_t\) array.
Plot the delta function with
\(Y_0 = 1/dt\) for \(0<dt\)
\(0\) otherwise.
using \(dt=0.1\) and t=np.linspace(0,20,20000)
###ANSWER HERE
Now we want to solve the differential equation for our Harmonic oscillator with a forced displacement of
\(Y_0 = 1/dt\) for \(0<dt\)
\(0\) otherwise.
using \(dt=0.1\) and t=np.linspace(0,20,20000). We will let the initial conditions be \(x(0)=v(0)=0\)
We need to make two minor changes to our solve_ivp call. First, we need to add a max_step=dt/10 keyword to our call (this is to make sure it sees the delta function which is only non-zero for a very little bit). Secondly, because we are using the GetDisplacement_delta function we will need to replace that using
sol = solve_ivp(equations, [t[0], t[-1]], init, t_eval=t, args=(params,D,GetDisplacement_delta),max_step=dt/10.)
Solve for and plot \(x(t)\) with the delta function displacement.
###ANSWER HERE
You may recall that in class, we worked out the Green’s function for the delta function potential as
where
and \(\gamma = \frac{c}{2m\omega_0}\) and \(\omega_0 = \sqrt{(k_1+k_2)/m}\).
Write a function G(t,c,m,omega), plot it on top of your differential equation solution and show that they largely agree.
Note: The Green’s function in this particular form breaks down in the overdamped case: \(\gamma>1\) gives you an imaginary frequency \(\omega\), which breaks down the cosine. There are ways around this, but you shouldn’t worry about it here.
###ANSWER HERE
Repeat the previous plots (for both the direct solve and the Green’s function) with a delta function impulse instead at \(t=10\). You will have to modify you’re greens function to be 0 when \(t<0\) and then adjust the call to it to be G(t-10,...)
###ANSWER HERE
b. Linear Combination of Delta Functions#
Now we want to write a general displacement function \(D_t\) as a linear combination of these \(\delta\) functions. To do this, we will generate from \(D_t\) a list D of DeltaParams so that for any time \(t\):
Y_0 = D_t(t)t_1 = tt_2 = t+dt
using dt=0.2 and t=np.arange(0,20,dt)
Apply this to the displacement from the previous exercise and check that it produces the correct thing by plotting the discretized version below.
def PlotD(D):
for d in D:
plt.bar(d.t1, d.Y0, width=d.t2 - d.t1, align='edge', alpha=0.5, edgecolor='black')
plt.ylabel("Y(t)")
plt.xlabel("t")
plt.plot(D_t)
PlotD(D)
plt.show()
###ANSWER HERE
Now we’ll indivudally solve for \(x_i(t)\) for each delta function directly and then some of the individual \(x_i(t)\) to make our final \(x(t)\).
Generate a \(x_i(t)\) for every delta forcing function in the list separately and then let \(x(t) = \sum_i x_i(t)\). Check that the \(x(t)\) that you get out is the same \(x(t)\) you got when you decomposed it into periodic forcing functions.
###ANSWER HERE
We can also solve for the linear combination of delta functions using our Green’s function. If the initial conditions were zero then we would have
where \(t_i\) and \(Y_{0i}\) are the time and height of the \(i\)’th delta function.
Since there are non-trivial initial conditions (\(x(0)=v(0)=1\)), we need to also solve the homogenous equation (e.g. the oscillator where there is no forced displacement) and add it to the solution with greens functions.
We can do this by evaluating with these intial conditions and using
def GetDisplacement_none(D,t):
return np.zeros_like(t)
Compute \(x(t)\) this way and compare against the result were it was computed directly.
###ANSWER HERE
Exercise 4. Square Well#
In this exercise, we are going to work with a square well potential. We will see how to decompose it into a series of periodic as well as delta function forcing functions.
a. The square well forcing function#
We will generate our square well using
def GetSquareWave(sampling_rate=500,frequency=0.2,duration=20,duty_cycle = 0.5,amplitude = 1):
# Time vector
t = np.linspace(0, duration, int(sampling_rate * duration), endpoint=False)
# Generate the square wave
square_wave = amplitude * signal.square(2 * np.pi * frequency * t, duty=duty_cycle)
return square_wave,t
Generate a square wave potential \(D_t\) and plot it over time.
###ANSWER HERE
b. Fourier transform of a square wave#
Now we’d like to take the fourier transform of the square wave. This can be done in the same way that you took the fourier transform of your sum of Harmonic forces above. This time, though, instead of just getting two frequencies, you should expect to get many frequencies. Go ahead and do this and plot both the absolute value and angle. When plotting the angle it might help to mask it by doing
mask=np.where(np.abs(the_fft)>1e-5*np.max(np.abs(the_fft)))
plt.plot(fft_freqs[mask],np.angle(the_fft[mask]),'o')
This just ensures that it doesn’t plot angles (which are poorly defined) for values of the fft that are essentially zero.
###ANSWER HERE
At this point, we need to figure out how to take this fft information and turn it into a sum over cosines.
Each frequency \(f\) corresponds to
where \(N\) is the sizes of the fft.
Use this information to produce a list of displace parameters myDisplacements that generate a set of cosines which sum to the square wave.
If you’ve done it correctly, you should get back the square wave by doing
myFunc=np.zeros(len(t),dtype=float)
for d in myDisplacements:
myFunc+=GetDisplacement(d,t)
plt.plot(t,myFunc)
plt.show()
###ANSWER HERE
c. Generating \(x(t)\) from many Harmonic Forces#
At this point, we have a list of Harmonic displacements. We now want to generate \(x(t)=\sum_i x_i(t)\) by solving the differential equation for each Harmonic forcing cosine function independently. Remember you need the initial conditions of all of your terms to sum to the initial conditions for \(x(t)\) (i.e. [1,1].)
###ANSWER HERE
d. Generate \(x(t)\) from many delta functions#
Now instead of decomposing our square wave into a sum of Harmonic waves, let’s instead decompose it into a sum of delta functions. Do this
Decomposing it into a sum of delta functions
Plotting your decomposition as a bunch of bars like you did above
Then computing \(x_i(t)\) for each of these delta functions
Sum them up to make \(x(t)=\sum_i x_i(t)\)
You should get the same results as in (c)
Make sure you use max_step/10 in your sol_ivp and you are using the appropriate parameters and equations for delta functions.
###ANSWER HERE
Acknowledgements:
Bryan Clark (original)