Two Springs and Forcing

Two Springs and Forcing#

import numpy as np
#from scipy.integrate import odeint
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

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


class DisplaceParams(NamedTuple):
  """A docstring"""
  Y0: float
  omega: float 
  phi: float = 0.0

class SHOParams(NamedTuple):
  """A docstring"""
  k: float
  m: float
  kp: float
  D: DisplaceParams
  c: float = 0.0
  initial_state: list = [0.0, 0.0]

Exercise 1. Two springs#

a. Equation of motion for forced displacement#

Let us now consider the scenario where we have two springs with springs constants $k$ (left) and $k'$ (right) with 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

\[ m \frac{d^2 x}{dt^2} + (k+k')x = k'y(t) \]

Recall that the natural frequency of this system is $\omega_0 = \sqrt{\frac{k+k'}{m}}$ In this exercise, let’s generically use a displacement

\[ y(t) = Y_0 \cos(\omega t) \]

Work out the coupled first order differential equations and write a function equations(t,state,params) which returns $\dot{x}$ and $\dot{v}$. The parameters can be set up as

params = SHOParams(k=1.0, m=1.0, kp=1.0, D=DisplaceParams(Y0=4.0, omega=0.25*2*np.pi), initial_state=[1.0, 1.0])

which corrresponds to

$k=k'=1.0$
$m=1.0$
$Y_0 = 4.0$
$L=0.1$
$ \omega = 0.25 (2\pi) \approx \omega_0+0.15$
initial state of $x=1$ and $v=1$

They can be accessed like params.kp or params.D.Y0

Once you’ve produced this function you can go ahead and numerically compute the solution of the differential equation using

sol = solve_ivp(equations, [t[0], t[-1]], params.initial_state, t_eval=t, args=(params,))
x = sol.y[0]
v= sol.y[1]

and plot both $x(t)$ and $Y(t)$ against $t$ for times t = np.linspace(0, 100, 400).

In order to do this, it will be useful to write a function

def GetDisplacement(displaceParams, t):
### do stuff
   return Y_t

which you call as GetDisplacement(params.D) and which will return a numpy array of the displacement.

To plot them both on the same axis you can then use

fig, ax1 = plt.subplots()
ax1.plot(t,x,zorder=2,label='x(t)',color='blue')
ax1.plot(t,GetDisplacement(params.D, t),zorder=1,label='y(t)',color='red',alpha=0.5)
ax1.set_xlabel("t")
ax1.set_ylabel("x(t)",color='b')
ax1.tick_params('y', colors='b')
plt.legend(loc='upper right')
plt.show()

You should see some clear beating going on due to being near but not on resonance.

###ANSWER HERE

Also go ahead and animate this system using

fig, ax, line, line2, circle = SetupAnimation()
ani = animation.FuncAnimation(fig, animate, frames=len(x), fargs=(x, Y_t,line,line2,circle,20), blit=True, interval=40, repeat=False)  # Pass parameters using fargs
display(HTML(ani.to_jshtml()))  
plt.close() 

To get the parameter Y_t you will need to call Y_t=GetDisplacement(params.D0) from above.

###ANSWER HERE

Finally, it’s interesting to consider tuning your system actually to resonance.

Let $\omega$ be the natural frequency $\sqrt{(k+k')/m}$.

Plot $x(t)$ over time $t$ for times t = np.linspace(0, 100, 400). Here you should see the oscillations simply get larger and larger over time.

In addition to generating this plot, let’s try to figure out how large $x(t)$ gets for various values of $\omega$. Loop over all values of $\omega$ in omegas=np.linspace(0.1,2,100) and plot the largest value of $x(t)$ as a function of $\omega$. To replace $\omega$ in your parameters you can do

new_params=params._replace(D=DisplaceParams(Y0=4.0, omega=omega))

###ANSWER HERE

b. Damping#

Our previous system was doing something unreasonable particularly as our driving frequency approached the natural frequency. This is partially because there was no damping in the system. Here, we will now introduce some dampling into our system. The new equations of motion will be

\[ m \frac{d^2 x}{dt^2} -cv (k+k')x = k' Y_0 \cos(\omega t) \]

Go ahead and modify equations(state,t,k,kp,m,Y0,omega,c) to solve the differential equations for this system. Then apply it to params = SHOParams(k=1.0, m=1.0, kp=1.0, c=0.5, D= D=DisplaceParams(Y0=4.0, omega=0.25*2*np.pi), initial_state=[1.0, 1.0])

Plot $x(t)$ and animate this system out to a time of $T=20$.

###ANSWER HERE

###ANSWER HERE

Now go ahead and plot the maximum $x(t)$ versus $\omega$. You should see that the damping now prevents the maximum height you reach from diverging.

###ANSWER HERE

While not required for credit, it is interesting to change the relevant parameters above and put it into the stiffness controlled regime and the mass controlled regime.

c. $G(\omega)$#

We have learned in class that the general solution of a damped driven harmonic oscillator of the form

\[ Y= F_0 \cos(\omega t+\phi) \]

after an initial transient is given by

\[ x(t) = F_0 |\tilde{G}(\omega)| \cos(\omega t + \phi + \textrm{angle}(\tilde{G}(\omega))) \]

where the greens function

\[ \tilde{G}(\omega) \equiv \frac{1}{k} \frac{1}{\left[1- \left(\frac{\omega}{\omega_n} \right)^2 + 2i \zeta \left(\frac{\omega}{\omega_n} \right) \right]} \]

with

\[ \zeta = \frac{c}{2m\omega_n} \]

Write a function G(omega,k,kp,omega_n,zeta) which computes $\tilde{G}(\omega)$.

Then use it to plot the steady state of $x(t)$ the driven damped oscillator above on top of the exact solution as a function of time. You expect to find that the transient roughly ends after time $2m/c$. Plot a vertical line there on your plot - i.e. plt.axvline(2m/c).

Notice that the non-transient values of $x(t)$ do not depend on the initial conditions in any way - i.e. it doesn’t show up in the formula at all.

###ANSWER HERE

It is also interesting to look at $\tilde{G}(\omega)$ as a function of $\omega/\omega_n$. Go ahead and make this plot both for $|\tilde{G}(\omega)|$ and $-\textrm{angle}[\tilde{G}(\omega)]$ for $\omega$ between 0 and 5.

###ANSWER HERE

d. Power Flow in the Steady State#

We are interested in solving for the power flow in the steady state. The power dissipated should go as

\[ P_\textrm{out} = cv^2 \]

over a single cycle.

The average power input goes as

\[ P_\textrm{in} = Fv \]

Using the same parameters as before **except with $c=1$, plot both the input and output power as a function of time using t=np.linspace(0, 20, 20000).

###ANSWER HERE

Now, we would like to evaluate the average power dissipated and put into the system. We want to average the power over one (or an integer) number of cycles.

For the paramaters above, we can get the index of an integer number of wavelengths by doing

t_idx_start=np.argmin(np.abs(t-10.0))
t_idx_stop=np.argmin(np.abs(t-20.))

Go ahead and evaluate the average of the power over these times. You can then compare it to the analytical formula

\[ \frac{c F_0^2 \omega^2 G^2}{2} \]

You should find that all three of these terms agree to 2-3 digits.

###ANSWER HERE

Exercise 2. 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.

\[ Y(t) = Y_1(t) + Y_2(t) \]

with an initial condition of $x=1.0; v=1.0$

We will let

$Y_1(t) = 4 \cos (2\pi (0.25) t)$
$Y_2(t) = 8 \cos (2\pi (0.1) t + 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’s let our initial conditions be $[0.3,0.6]$ and $[0.7,0.4]$ respectively.

For the other parameters, continue using

params = SHOParams(k=1.0, m=1.0, kp=1.0, c=0.5, D=..., initial_state=...)

with times t= np.linspace(0, 20, 10000)

Solve for and plot $x_1(t)$, $x_2(t)$ and $x(t)$. Your answer should look like ….

###ANSWER HERE

Suppose we just had the curve

### <-- you don't get to see this next line. Someone just hands you D_t
D_t=GetDisplacement(DisplaceParams(Y0=4.0, omega=0.25*2*np.pi, phi=0), t)+GetDisplacement(DisplaceParams(Y0=8.0, omega=0.1*2*np.pi, phi=1.0), t)  
###
plt.plot(D_t)
plt.show()

on a plot but we didn’t know what frequencies it came from.

In order to solve for the frequency $x(t)$ using our previous techniques, we would need to decompose our displacement $Y(t)$ into a linear combination of cosines with different frequences and phi. - i.e. $$ Y(t) = \sum_i Y_{0,i} \cos (\omega_i t + \phi_i ) $$

We can do this using a FFT. Let’s call

fft_freqs=fftfreq(len(D_t),t[1]-t[0])
the_fft=fft(D_t)

Then we can 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),'.')
plt.xlim(-1.2,1.2)
plt.show()
mask=np.where(np.abs(the_fft)>cutoff)
print("Frequencies: ",fft_freqs[mask])

These last two lines are just to indicate which frequencies are large. You should see (both from the graph and printing) that we get back $\pm 0.25$ and $ \pm 0.1$ as large frequencies.

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

###ANSWER HERE

From these graphs above we can extract that our curve was

\[ 4.0 \exp[2 \pi i t \times 0.25] + 4.0 \exp[-2 \pi i t \times 0.25] + 2.0 \exp[-2\pi i t \times 0.1 + 1] + 2.0 \exp[-2 \pi i t \times 0.1 - 1] \]

\[ = 8.0 \cos (2\pi i t \times 0.25) + 4.0 \sin (2\pi i t \times 0.1 + 1) \]

This is what was expected given what we put in.

Exercise 3. 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#

First, similar to DisplaceParams let us define a DeltaParams as

class DeltaParams(NamedTuple):
  """A docstring"""
  Y0: float
  t1: float 
  t2: float

which creates a impulse of height $Y_0$ for times $t_1 < t < t_2$.

Then we can create another parameter class

class SHOParams_delta(NamedTuple):
  """A docstring"""
  k: float
  m: float
  kp: float
  D: DeltaParams
  c: float = 0.0
  initial_state: list = [0.0, 0.0]

which is similar to our last class but takes a DeltaParams.

Now modify equations(t, state, shoParams) to instead implement the first order differential equations using this “delta function” $Y(t)$.

Defining $dt=0.1$, let us compute $x(t)$ for a potential with $Y_0 =1/dt$ and $0<t<dt$ on a set of t=np.arange(0,20,dt). For technical reasons, we need to add a max_step parameter to our solve_ivp function - i.e.

sol = solve_ivp(equations, [t[0], t[-1]], params.initial_state, t_eval=t, args=(params,),max_step=dt/10.)

Recall that we know the delta function potential should give a

\[ G(\omega)= \exp\left(\frac{-ct}{2m}\right) \frac{\sin(\omega_d t)}{m\omega_d} \]

where

\[ \omega_d = \omega_0 \sqrt{1-\zeta^2} \]

Plot this on top of your differential equation solution and show that they largely agree.

Note: After you use $dt=0.2$ you can rerun this exercise with $dt=0.1$ and it will look somewhat better but I think the $dt=0.2$ let’s you see what’s going on a little better when you decomose your forcing function.

###ANSWER HERE

Now let’s go ahead and instead put a delta function impulse at $t=10$.

###ANSWER HERE

b. Linear Combination of Delta Functions#

Now we want to write $D_t$ as a linear combination of these $\delta$ functions. To do this, we will generate a list D of DeltaParams so that for every time

Y_0 = D_t(t)
t_1 = t
t_2 = t+dt

Check that this produced 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(t,D_t)
PlotD(D)
plt.show()

###ANSWER HERE

Now we want to go ahead and generate a $x_i(t)$ for every forcing function in the list separately and then let $x(t) = \sum_i x_i(t)$. You need the initial conditions of all of your terms to sum to the initial conditions for $x(t)$ (i.e. [1,1]); I think the easiest way to do this is to make the first term have these initial conditions and the other to have $[0,0]$. Go ahead and compute all the $x_i(t)$ and plot them to give $x(t)$. Check that the $x(t)$ that you get out is the same $x(t)$ when you decomposed it into periodic forcing functions.

###ANSWER HERE

Exercise 4. Square Well#

In this exercise, we are going to work with a square well potential. We will see how both to decompose it into a series of periodic forcing functions 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):
    # Parameters
    #frequency =   0.2  # Frequency in Hz
    amplitude = 1
    #duration = 20 # Duration in seconds
    #sampling_rate = 500  # Samples per second
    duty_cycle = 0.5 # Duty cycle (0 to 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 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

\[ A \exp[2\pi i \times f] = |A|/N \cos(2\pi f + arg(A)) + |A|/N i \sin (...) \]

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 now 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].) Also remember to use SHOParams and not SHOParams_delta and to use the equations for the Harmonic displacement.

###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