Harmonic Oscillator

Harmonic Oscillator#

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

Show code cell content

Hide code cell content

def draw_spring(x0, y0, x1, y1, coils=10, amplitude=0.1):
    """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

def draw_spring_system(spring_loc):
    xs,ys=draw_spring(-2, 0,spring_loc ,0,coils=10)
    plt.plot([-2.0,-2.0],[-0.25,0.25],color='red',linewidth=3,zorder=1)
    plt.plot(xs,ys,color='blue',linewidth=2,zorder=0)
    radius = 0.05
    circle = patches.Circle((spring_loc,0), radius, color='blue', fill=True)
    plt.gca().add_artist(circle)
    plt.axis('equal')
    plt.xlim(-3, 3)
    plt.show()

def animate_spring_system(x):
    fig, ax = plt.subplots()
    ax.plot([-2.0,-2.0],[-0.25,0.25],color='red',linewidth=3,zorder=1)
    line, = ax.plot([], [], lw=2, color='blue',zorder=0)  # Initialize an empty line object

    radius = 0.05
    circle = patches.Circle((0,0), radius, color='blue', fill=True)
    ax.add_artist(circle)
    top_val=np.max(np.abs(x))*1.1
    ax.set_xlim(-3, top_val)    
    ax.set_ylim(-1, 1)  

    def animate(i):
        x1 = x[i]  # x-coordinate changes over time
        xs, ys = draw_spring(-2, 0, x1, 0.0, coils=8, amplitude=0.10) 
        line.set_data(xs, ys)
        circle.center = (x1,0)
        return line,circle,

    ani = animation.FuncAnimation(fig, animate, frames=len(x), blit=True, interval=20,repeat=False) 
    ###save animation as html
    #ani.save('spring_system.html', writer='html')
    display(HTML(ani.to_jshtml()))  
    plt.close() 

Our goal in this assignment is to understand the Harmonic Oscillator.

Exercise 1. A Single Spring#

a. Simple Harmonic Oscillator#

In this part, we will be working with a simple harmonic oscillator. For the moment, we should think of this as a spring connected to a wall on one end and to a mass \(m\) on the other. We choose the natural equilibrium length of the spring to be 2 and let the x-coordinate at the equilbrium length to be \(x=0\) (so the deviation from the equilibrium and the x-coordinate match up). To draw this spring go ahead an plot

draw_spring_system(0)

where the 0 represents that it is at \(x=0\).

###ANSWER HERE

The relevant differential equation for the simple Harmonic oscillator is

\[ m \ddot{x} = -k x \]

where \(x\) represents the amount the spring is displaced away from its equilibrium position.

To solve this with python, it will be convenient to instead write it out as two first-order equations,

\[\begin{split}\begin{align} \dot{x} &= v \\ \dot{v} &= -\frac{k}{m} x \end{align} \end{split}\]

where \(v\) is the velocity of the mass.

These equations can be directly translated into python:

def equations(t,state,k,m):
  x,v = state
  x_dot = v
  v_dot = -k/m*x
  return [x_dot,v_dot]

We can now get python to solve this differential equation for us. First, we need to give it some initial conditions. For example, we can stretch the spring out to \(x=1\) and set the initial velocity to \(v=0\) by writing

initial_state = [1, 0]

Then we need to specify for which times \(t\) we want to evaluate our differential equation with respect to. Here let’s choose a grid of 400 time points between 0 and \(25\pi\) as

t = np.linspace(0, 25*np.pi, 400)

Finally, we need to pick some values for the parameters

k = 1.3
m = 2.25

We can now solve our differential equation with

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

Go ahead and solve this differential equation and plot it with

fig, ax1 = plt.subplots()
ax1.plot(t,x) #<-- This is the important line which plots.
ax1.set_xlabel("t")
ax1.set_ylabel("x(t)",color='b')
ax1.tick_params('y', colors='b')  
plt.show()

Notice that your wave-length \(\lambda = \frac{2\pi}{\omega_0}\) where the natural frequency \(\omega_0 = \sqrt{k/m}\).

Also go ahead and plot the velocity. You can plot it separately or if you want to plot it on the same figure so you can match the velocity and time before the plt.show() go ahead and do

ax2 = plt.gca().twinx()  
ax2.plot(t, v, 'r-',alpha=0.2)  #<--- you can turn up the alpha if it's too dim but I find it annoying to have them both be strong  
ax2.set_ylabel('v(t)', color='r') 
ax2.tick_params('y', colors='r') 

###ANSWER HERE

b. Computing the energy#

Compute the energy

\[ E = \frac{1}{2}mv^2 + \frac{1}{2}kx^2 \]

of your system and show that it is constant. Beware of the plot scale; it helps if you keep your y-limit from zero to one - i.e. ax1.set_ylim(0,1).

###ANSWER HERE

c. Animating your spring#

Now animate your spring. You can do this by calling

animate_spring_system(x)

on your x-coordinate.

In some cases throughout this notebook, animating might take an annoying while. So instead, you may animate with a lower “frame rate” by calling this function on x[0::2], for example, which samples every 2 elements of \(x\).

###ANSWER HERE

d. An underdamped Harmonic Oscillator#

Now let’s introduce some damping. The equation of motion for a damped harmonic oscillator is

\[ m \ddot{x} = -k x - c v \]

Rewrite this as a series of first-order differential equations and modify your equations(state,t,x,k,m,c) function to solve this differential equation using the same starting conditions and values for \(k\) amd \(m\) as earlier. Start by using \(c=0.2\). Recall that the dissipation \(\gamma=c/(2m\omega_0)\), and \(\omega_0 = \sqrt{k/m}\) is the undamped frequency. When \(\gamma<1\), as is this case, we have an underdamped oscillator.

Plot the position and velocity as a function of time. Also animate your spring.

###ANSWER HERE

Now, we would like to check two things explicitly about the damped oscillator:

We expect that the period of the damped oscillator should be \(2\pi/\omega_d\) where \(\omega_d \equiv \omega_0 \sqrt{1-\gamma^2}\) and \(\gamma = c/(2m\omega_0)\). Verify this by extracting the peak locations of the undamped oscillator.

To do this, we can use my_peaks=find_peaks(x)[0].tolist() and then use t[my_peaks] to get the time of the respective peaks. Check that those times are all equidistant and that they happen with the prescribed period. The answer you will actually get if you do this will be reasonable but not super-accurate. You can do better by (a) increasing the number of points in your linspace to 4000 (but don’t rerun the animation becuase it will be very slow) and/or instead of using the difference in time between two peaks, computing the slope of peaks x time elapsed using np.polyfit(range(0,len(my_peaks)),t[my_peaks],1).
Finally, verify that the damped oscillator’s strength decays exponentially. We can do this by looking at the peak values (x[my_peaks]) and plot them. To tell that something decreases exponentially, you should plot it on a semi-log scale (ax1.set_yscale('log')) and see that the plot is linear.

###ANSWER HERE

###ANSWER HERE

###ANSWER HERE

e. An overdamped oscillator#

Now we turn to the overdamped oscillator. Find \(c\) that gives \(\gamma=1.2\) and solve the corresponding ODE. Then consider some value where \(\gamma \gg 1\). For both scenarios plot the position and velocity vs time and animate them.

###ANSWER HERE

###ANSWER HERE

f. Gravity#

Now suppose you turn your system vertically in such a way that gravity acts as a constant force on your oscillator. Modify your equations of motion to include both the force from gravity and a dissipation term.

By running your simulation, compute the equilibrium value for masses masses=[1.0,2.0,3.0,4.0]. Then plot the equilibrium spring location versus the mass and show that the equilbrium position is linear with the mass. You may find the equilibrium simply by averaging the position, but be sure to disregard any transient effects.

How could you then use this information as a scale?

For concreteness, you may use \(g=10,k=1.3,c=0.2\), but that should make no qualitative difference. In fact, your choice of \(c\) should make no quantitative difference (why?).

###ANSWER HERE

###ANSWER HERE

Using this information, by calculating the equilibrium position we immediately determine the mass. So we have a scale.

Exercise 2. Two springs#

In this assignment, we are going to work on the simple Harmonic oscillator for two springs. We will be using slightly different animation and spring-drawing functions in this assignment.

a. Equation of motion for forced displacement#

Let us consider the scenario where we have two springs with 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

\[ m \frac{d^2 x}{dt^2} + (k_1+k_2)x = k_2 y(t) \tag{1} \]

Recall that the natural frequency of this system is \(\omega_0 = \sqrt{\frac{k_1+k_2}{m}}\). In this exercise, let’s generically use a displacement

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

There will be many quantities being used in this exercise, so we choose to condense the relevant parameters by defining the following classes:

class SpringParams(NamedTuple):
  k1: float
  m: float
  k2: float
  c: float = 0.0

class HarmonicDisplaceParams(NamedTuple):
  Y0: float
  omega: float 
  phi: float = 0.0

The parameters can be set up as

params = SpringParams(k1=1.0, m=1.0, k2=1.0)
D=DisplaceParams(Y0=4.0, omega=0.25*2*np.pi)
initState=[1.0,1.0]

which corrresponds to

\(k1=k2=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.k2 or D.Y0

Start by writing a function

def GetDisplacement_harmonic(D,t):
    ### do stuff
    return D_t

which takes the displacement parameters and an array of times t = np.linspace(0, 100, 400) and returns an array of displacements.

Go ahead and plot the displacements for the target Harmonic displacements.

###ANSWER HERE

Turn our second order differential equation into a system of coupled first order differential equations in \((x,v)\) and write a function

def equations(t, state, springParams, D, GetDisplacement=GetDisplacement_harmonic):

    return [xdot,vdot]

which takes a time t, the state, the spring and the displacement parameters and returns returns \(\dot{x}\) and \(\dot{v}\). By default it uses the function GetDisplacement_harmonic to generate the displacement.

###ANSWER HERE

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,D,))
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).

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_harmonic(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 exactly on resonance.

###ANSWER HERE

Now 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.D,t) from above.

###ANSWER HERE

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

Let \(\omega=\omega_0\), where \(\omega_0\) is the natural frequency defined as \(\omega_0=\sqrt{(k_1+k_2)/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\).

###ANSWER HERE

###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 damping into our system. The new equations of motion will be

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

Modify equations(t,state,params,D,GetDisplacement=GetDisplacement_harmonic) to solve the differential equations for this system.

Then apply it to

params = SpringParams(k1=1.0, m=1.0, k2=1.0, c=0.5)
D=DisplaceParams(Y0=4.0, omega=0.25*2*np.pi)
initial_state=[1.0, 1.0]

Plot \(x(t)\) and the forcing \(y(t)\), then animate this system out to a time of \(t=20\).

###ANSWER HERE

###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 (\(c \gg m\)) and the mass controlled regime (\(m \gg c\)).

c. Green’s Function \(G(\omega)\)#

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

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

after an initial transient is given by

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

where the Green’s function

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

with

\[ \gamma = \frac{c}{2m\omega_0} \]

Write a function G(omega,k,omega_0,gamma) 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 within an order of magnitude of \(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

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

Now investigate what happens as \(c\to 0\). Redo the previous plot with \(c=0.001\). It might be helpful to plot \(-\textrm{angle}[\tilde{G}(\omega)]\) separately, so you can see its features clearly.

###ANSWER HERE

From their behavior, you should be able to tell that the magnitude of \(G\) is (proportional to) the derivative of the angle of \(G\). As the derivative blows up, the undifferentiated function – here the angle – approaches a discontinuity! Although continuity is predominant in physics, and we’re often taught to disregard discontinuities as pathologies, this exact behavior happens in many places and is behind many interesting and deep phenomena such as phase transitions and shock waves.

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 power input goes as

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

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

Also, this is a good place to play with the \(c\) parameter and see and try to understand what is happening. Increase the range of \(t\) if you feel you need to.

###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 periods by doing

find_peaks(-getDisplacement(params.D,t))[0]

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

\[ \overline{P}=\frac{c\omega^2 Y_0^2 G^2}{2} \]

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

###ANSWER HERE

Acknowledgements:

Bryan Clark (original)