Problem 3.1: Damped, driven harmonic oscillator described by:

$m\ddot{x} + \gamma \dot{x} + kx = e^{i\omega t}$

A. Under what conditions will the governing equations for small displacements of a particle around an arbitrary 1D potential minimum be simple undamped harmonic motion?

when:

$e^{i\omega t} = \gamma \dot{x}$

we can reduce the equation to a simple, undamped harmonic oscillator:

$m\ddot{x} + kx = 0$

assuming x is of the form:

$x = Ae^{i\omega t}$

we have:

$\dot{x} = Aiwe^{i\omega t}$

if:

$\gamma = \dfrac{1}{Aiw}$

then this condition is satisfied.

I misunderstood the question, the answer is that the potential needs to be quatratic for small motions.

taylor series:

$V = V_0 + V_1x + 1/2V_2x^2 ...$

$dV/dt = V_1 + V_2x ...$

$m\ddot{x} = -V_1 - V_2x ...$

V_1 goes to zero at min

$m\ddot{x} = -V_2x ...$

assume no higher order terms

$m\ddot{x} = -V_2x$

$m\ddot{x} + V_2x = 0$

simple harmonic oscillator as long as $V_2 > 0$.

B. Find the solution to the homogeneous equation, and comment on the possible cases. How does the amplitude depend on the frequency?

homogenous equation:

$m\ddot{x} + \gamma \dot{x} + kx = 0$

assuming x of the form:

$x = Ae^{\omega t}$

we have:

$\dot{x} = A\omega e^{\omega t}$

$\ddot{x} = A\omega^2e^{\omega t}$

$m(A\omega^2e^{\omega t}) + \gamma(A\omega e^{\omega t}) + k(Ae^{\omega t}) = 0$

$m\omega^2e^{\omega t} + \gamma \omega e^{\omega t} + ke^{\omega t} = 0$

$m\omega^2 + \gamma \omega + k = 0$

using the quadratic equation:

$\omega = \dfrac{-\gamma \pm \sqrt{\gamma^2 - 4mk}}{2m}$

Back to $x(t)$, we have:

$x(t) = Ae^{(\dfrac{-\gamma \pm \sqrt{\gamma^2 - 4mk}}{2m} t)}$

$x(t) = Ae^{-\gamma t/2m}e^{(\pm \sqrt{\gamma^2 - 4mk}t/2m)}$

let:

$a = \gamma^2 - 4mk$

case 1:

if $a = 0$:

$x(t) = Ae^{-\gamma t/2m}$

this is exponential decay due to critical damping:

$0 = \gamma^2 - 4mk$

$\gamma = 2\sqrt{mk}$

case 2:

if $a < 0$:

$x(t) = Ae^{-\gamma t/2m}e^{(\pm i\sqrt{|\gamma^2 - 4mk}|t/2m)}$

Applying Euler's identity:

$x(t) = Ae^{-\gamma t/2m}(cos(\pm\sqrt{|\gamma^2 - 4mk}|t/2m) + isin(\pm\sqrt{|\gamma^2 - 4mk}|t/2m))$

taking only the real component we have an exponentially decaying oscillator, the underdamped case:

$x(t) = Ae^{-\gamma t/2m}cos(\pm\dfrac{\sqrt{|\gamma^2 - 4mk}|t}{2m})$

case 3:

if $a > 0$:

$x(t) = Ae^{-\gamma t/2m}e^{(\pm \sqrt{|\gamma^2 - 4mk}|t/2m)}$

this will always simplify to a negative exponential because $\gamma > \sqrt{|\gamma^2-4mk|}$, it represents the overdamped case.

C. Find a particular solution to the inhomogeneous problem by assuming a response at the driving frequency, and plot its magnitude and phase as a function of the driving frequency for $m = k = 1$, $\gamma = 0.1$.

$m\ddot{x} + \gamma \dot{x} + kx = e^{i\omega t}$

particular solution:

$x(t) = Ae^{i\omega t}$

$-A\omega^2 m e^{i\omega t} + Ai\omega \gamma e^{i\omega t} + Ake^{i\omega t} = e^{i\omega t}$

$Ae^{i\omega t} (-\omega^2 m + i\omega \gamma + k) = e^{i\omega t}$

$A = \dfrac{1}{-\omega^2 m + i\omega \gamma + k}$

import matplotlib.pyplot as plt
import numpy as np

m = 1
k = 1
gamma = 0.1

w = np.linspace(0, 2, 10000)
A = 1/(-w*w*m + w*1j*gamma + k)

plt.plot(w, np.absolute(A))

plt.ylabel('amplitude')
plt.xlabel('omega')
plt.show()

plt.plot(w, np.angle(A))
plt.ylabel('phase')
plt.xlabel('omega')
plt.show()

D For a driven oscillator the Q or Quality factor is defined as the ratio of the center frequency to the width of the curve of the average energy (kinetic + potential) in the oscillator versus the driving frequency (the width is defined by the places where the curve falls to half its maximum value). For an undriven oscillator the Q is defined to be the ratio of the energy in the oscillator to the energy lost per radian (one cycle is 2π radians). Show that these two definitions are equal, assuming that the damping is small. How long does it take the amplitude of a 100 Hz oscillator with a Q of 109 to decay by 1/e?

E Now find the solution to the driven oscillator by using Laplace transforms. Take the initial condition as $x(0) = \dot{x}(0) = 0$.

import numpy as np
import math
from scipy.integrate import odeint

#z = [x', x]
#z' = [x'', x']
m = 1
k = 1
gamma = 0.1
i = 1j
omega = 2
z0 = [0, 0] #[x', x]


t = np.linspace(0, 50, 10000)

def deriv(z, t):
    return np.array((
         1/m*(-gamma*z[0] - k*z[1] + math.cos(omega*t)), # math.exp(i*omega*t)),  # this is z[0]'
         z[0] # this is z[1]'
       ))

dx, x = odeint(deriv, z0, t).T

plt.plot(t, x)

plt.ylabel('x(t)')
plt.xlabel('t')
plt.ylim([-1,1])
plt.show()

F For an arbitrary potential minimum, work out the form of the lowest-order correction to simple undamped unforced harmonic motion.

$V(x) = V_0 + 1/2V_2x^2 + 1/6V_3x^3$

$m\ddot{x} + kx + \epsilon x^2 = 0$

let:

$x = x_0 + \epsilon x_1 + O(\epsilon^2)$

$x_0 = e^{i\omega t}$

$m(\ddot{x_0}+ \epsilon \ddot{x_1}) + k(x_0+ \epsilon x_1) + \epsilon (x_0+ \epsilon x_1)^2 = 0$

remove $\epsilon^2$ terms:

$m\ddot{x_0}+ m\epsilon \ddot{x_1} + kx_0+ k\epsilon x_1 + \epsilon (x_0^2+ 2x_0x_1\epsilon) = 0$

$m\ddot{x_0}+ m\epsilon \ddot{x_1} + kx_0+ k\epsilon x_1 + x_0^2\epsilon = 0$

$-m\omega^2e^{i\omega t}+ m\epsilon \ddot{x_1} + ke^{i\omega t}+ k\epsilon x_1 + e^{2i\omega t}\epsilon = 0$

Problem 3.2: Explicitly solve (and try to simplify) the system of differential equations for two coupled harmonic oscillators (see Figure 3.2; don’t worry about the initial transient), and then find the normal modes by matrix diagonalization.

$F_1 = -kx_1 + k(x_2-x_1)$

$F_2 = -kx_2 - k(x_2-x_1)$

$F_1 = k(-2x_1 + x_2)$

$F_2 = k(x_1-2x_2)$

$\begin{bmatrix} F_1 \\ F_2 \end{bmatrix} = \begin{bmatrix} -2k & k \\ k & -2k \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} $

$\begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix}\begin{bmatrix} \ddot{x_1} \\ \ddot{x_1} \end{bmatrix} = \begin{bmatrix} -2k & k \\ k & -2k \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} $

$\begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix}\ddot{x} - \begin{bmatrix} -2k & k \\ k & -2k \end{bmatrix} x = 0$

$\ddot{x} - \begin{bmatrix} -2k & k \\ k & -2k \end{bmatrix} \begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix}^{-1}x = 0$

$\ddot{x} - \begin{bmatrix} -2k & k \\ k & -2k \end{bmatrix} \begin{bmatrix} 1/m_1 & 0 \\ 0 & 1/m_2 \end{bmatrix}x = 0$

$\ddot{x} - \begin{bmatrix} -2k/m_1 & k/m_2 \\ k/m_1 & -2k/m_2 \end{bmatrix} x = 0$

$A = \begin{bmatrix} 2k/m_1 & -k/m_2 \\ -k/m_1 & 2k/m_2 \end{bmatrix}$

$m_1 = m_2$

$\omega_0^2 = k/m$

$A = \begin{bmatrix} 2\omega_0^2 & -\omega_0^2 \\ -\omega_0^2 & 2\omega_0^2 \end{bmatrix}$

$0 = |A-\lambda I|$

$0 = det(\begin{bmatrix} 2\omega_0^2 -\lambda_0 & -\omega_0^2 \\ -\omega_0^2 & 2\omega_0^2 -\lambda_1 \end{bmatrix})$

let: $a = 2k/m_1$, $b = -k/m_2$, $c = -k/m_1$, $d = 2k/m_2$, $D=ad-bc$, $T=a+d$

$D = \begin{bmatrix} T/2 + (T^2/4-D)^{1/2} & 0 \\ 0 & T/2 - (T^2/4-D)^{1/2} \end{bmatrix}$

$M = \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}$

Problem 3.3: A common simple digital filter used for smoothing a signal is:

$y(k) = \alpha y(k-1) + (1-\alpha)x(k)$

where $\alpha$ is a parameter that determines the response of the filter. Use z-transforms to solve for $y(k)$ as a function of $x(k)$ (assume $y(k < 0) = 0$).

$Z\{y(k)\} = \sum_{k=0}^{\infty}(\alpha y(k-1) + (1-\alpha)x(k) )z^{-k}$

$Z\{y(k)\} = \sum_{k=0}^{\infty}\alpha y(k-1)z^{-k} + \sum_{k=0}^{\infty}(1-\alpha)x(k)z^{-k}$

$Z\{y(k)\} = \alpha \sum_{k=0}^{\infty}y(k-1)z^{-k} + (1-\alpha)\sum_{k=0}^{\infty}x(k)z^{-k}$

let $k' = k-1$:

$Z\{y(k)\} = \alpha \sum_{k'=-1}^{\infty}y(k')z^{-k'}z^{-1} + (1-\alpha)\sum_{k=0}^{\infty}x(k)z^{-k}$

$Z\{y(k)\} = \dfrac{\alpha}{z} \sum_{k'=-1}^{\infty}y(k')z^{-k'} + (1-\alpha)\sum_{k=0}^{\infty}x(k)z^{-k}$

$Z\{y(k)\} = \dfrac{\alpha}{z} (\sum_{k'=0}^{\infty}y(k')z^{-k'} + y(-1)z) + (1-\alpha)\sum_{k=0}^{\infty}x(k)z^{-k}$

$y(-1) = 0$

$Z\{y(k)\} = \dfrac{\alpha}{z} \sum_{k'=0}^{\infty}y(k')z^{-k'} + (1-\alpha)\sum_{k=0}^{\infty}x(k)z^{-k}$

$Y(z) \equiv Z\{y(k)\}$ , $X(z) \equiv Z\{x(k)\}$

$Y(z) = \dfrac{\alpha}{z}Y(z) + (1-\alpha)X(z)$

$Y(z) (1 - \dfrac{\alpha}{z}) = (1-\alpha)X(z)$

$Y(z) = \dfrac{(1-\alpha)}{(1 - \dfrac{\alpha}{z})}X(z)$

What is the amplitude of the frequency response?

$H(z) = \dfrac{Y(z)}{X(z)}$

$H(z) = \dfrac{(1-\alpha)}{(1 - \dfrac{\alpha}{z})}$