Let's start with the problem statement:

Consider the motion of a damped, driven harmonic oscillator (such as a mass on a spring, a ball in a well, or a pendulum making small motions):

$$m\ddot{x} + \gamma \dot{x} + k x = e^{i\omega t} \tag{4.58} \label{eq:problem}$$

(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?

In physics, forces and potentials follow this relationship:

$$F(\mathbf{x}) = - \nabla U(\mathbf{x}) \tag{1}$$

in which $\mathbf{x}$ can be expressed in 3D, and $\nabla$ is the gradient operator. In 1D, this simplifies to:

$$F(x) = - \frac{\rm d}{{\rm d}x} U(x) \tag{2}$$

Combining this with Newton's second law, one obtains:

$$- \frac{\rm d}{{\rm d}x} U(x) = m \ddot{x} \tag{3} \label{eq:pot_newton}$$

Without loss of generality, we can assume the coordinate system is centered around the local minimum of $U(x)$. If this is not the case, a simple change of coordinate $x'=x-x_0$ can be used at the end of this reasoning, if $x_0$ is the local minimum.

Combining equations (4.58) and (3) assuming $\gamma = 0$ and no driving force, we get:

$$- \frac{\rm d}{{\rm d}x} U(x) = m \ddot{x} = -kx$$$$\Rightarrow U(x) = \frac{k x^2}{2} + U(0) \tag{4}$$

where $U_0$ is an arbitrary additive constant, not affecting the physics of the system.

In a more general case, we can show that this result holds as long as $U(0)$ is a local minimum and $x$ is small. We start with a Taylor expansion of $U(x)$:

$$U(x) \approx U(0) + x U'(0) + \frac{x^2}{2} U''(0) + \frac{x^3}{6} U'''(0) + \dots \tag{5}$$

which is consistent with equation (4) only if $U'(0)=0$, $U''(0)=k$, and all terms of power higher than 2 can be neglected, which is reasonable under a small $x$. Intuitively, the constraint $U'(0)=0$ ensures no constant force exists across the entire domain, which would lead to strange results anyway; it would make $x$ drift to infinity in the absence of a spring.

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

The homgeneous equation is:

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

This is a linear system of derivatives of $x$, so we can assume a solution of the form $x=e^{rt}$ with $r \in \mathbb{C}$, leading to:

$$m r^2 e^{rt} + \gamma r e^{rt} + ke^{rt} = 0 \tag{7}$$

This equation is true $\forall t$ iff:

$$m r^2 + \gamma r + k = 0 \tag{8}$$

Which is a second degree polynomial with solutions:

$$r = \frac{-\gamma \pm \sqrt{\gamma^2 - 4mk}}{2m} \tag{9}$$

Reorganizing this and defining $\lambda = \frac{\gamma}{2m}$, we obtain:

$$r = \lambda \left(-1 \pm \sqrt{1-\frac{k}{\lambda^2}} \tag{10}\right)$$

Let's break down the different scenarios:

case 1: $\gamma = 0$, undamped oscillation

In this case, equation (10) has two purely imaginary solutions:

$$r = \pm i \sqrt{\frac{k}{m}} = \pm i \omega_0 \tag{11}$$

where $\omega_0$ is the natural frequency of the system. General solutions to the ODE are then:

$$x(t) = A e^{i\omega_0 t} + B e^{-i\omega_0 t} \tag{12}$$

with $A, B \in \mathbb{C}$ are arbitrary constants to be determined from initial conditions. We are only looking for the real part of this family of solutions; by redifining new constants $A'$ and $B'$, we can rewrite this as:

$$x(t) = A' \frac{e^{i\omega_0 t}+e^{-i\omega_0 t}}{2} + B' \frac{e^{i\omega_0 t}-e^{-i\omega_0 t}}{2i} = A' \cos(\omega_0 t) + B' \sin(\omega_0 t)\tag{13}$$

Which is a familiar oscillator. Let's plot a solution for $x(0)=1,~\dot x(0)=0$ using numpy and matplotlib:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation
import matplotlib
import math
from IPython.display import HTML

matplotlib.rcParams["animation.bitrate"] = 1600

t_s = 10
dt_s = 0.016

k = 5
m = 1

omega_0 = math.sqrt(k/m)

t = np.arange(0, t_s, dt_s)
x = np.cos(omega_0 * t)

fig, axs = plt.subplots(1, 2)
axs[0].axis([0,t_s,-1.5,1.5])
axs[0].set_xlabel("t [s]")
axs[0].set_ylabel("x(t) [m]")
axs[1].axis([-1, 1, -1.5, 1.5])
l1, = axs[0].plot([], [], 'b-')
l2, = axs[1].plot([], [], 'b-')
l3, = axs[1].plot([], [], 'b.', ms=40)

def animate(i):
    l1.set_data(t[:i], x[:i])
    l2.set_data([0, 0], [0, x[i]])
    l3.set_data([0], [x[i]])

ani = matplotlib.animation.FuncAnimation(fig, animate, frames=len(t), interval=int(dt_s*1e3))

plt.close()
HTML(ani.to_html5_video())

case 2: $1 - \frac{k}{\lambda^2} > 0$, overdamped

This happens when $\gamma > 2m\sqrt{k}$. In those cases, equation (10) leads to two real solutions:

$$r = -\lambda \pm \sqrt{1-\frac{k}{\lambda^2}} = \lambda \pm \rho \tag{14}$$

Leading to general solutions:

$$x(t) = A e^{(-\lambda + \rho) t} + B e^{(-\lambda - \rho) t} = e^{-\lambda t} \left( A e^{\rho t} + B e^{-\rho t}\right) \tag{15}$$

This family of solutions can only converge if $A = 0$. Let's plot one example with the same settings as earlier and $\gamma = 5$:

k = 5
m = 1
gamma = 5

lambd = gamma / (2*m)
rho = math.sqrt(1-k/(lambd*lambd))

t = np.arange(0, t_s, dt_s)
x = np.exp(- (lambd + rho) * t)

fig, axs = plt.subplots(1, 2)
axs[0].axis([0,t_s,-1.5,1.5])
axs[0].set_xlabel("t [s]")
axs[0].set_ylabel("x(t) [m]")
axs[1].axis([-1, 1, -1.5, 1.5])
l1, = axs[0].plot([], [], 'b-')
l2, = axs[1].plot([], [], 'b-')
l3, = axs[1].plot([], [], 'b.', ms=40)

def animate(i):
    l1.set_data(t[:i], x[:i])
    l2.set_data([0, 0], [0, x[i]])
    l3.set_data([0], [x[i]])

ani = matplotlib.animation.FuncAnimation(fig, animate, frames=len(t), interval=int(dt_s*1e3))

plt.close()
HTML(ani.to_html5_video())

case 3: $1 - \frac{k}{\lambda^2} < 0$, damped oscillations

This implies $\gamma < 2m\sqrt{k}$. Two complex solutions for $r$ can be found:

$$r = -\lambda \pm \sqrt{1 - \frac{k}{\lambda^2}} = -\lambda \pm \sqrt{-1}\sqrt{\frac{k}{\lambda^2}-1} := -\lambda \pm i \omega_d \tag{16}$$

where $\omega_d$ is the natural frequency of the damped oscillator. The general solutions are now:

$$x(t) = A e^{(-\lambda + i\omega_d) t} + B e^{(-\lambda - i\omega_d) t} = e^{-\lambda t} \left( A e^{i \omega_d t} + B e^{-i \omega_d t}\right) \tag{17}$$

Using new factors $A'$ and $B'$, this can be rewritten as:

$$x(t) = e^{-\lambda t} \left( A' \cos(\omega t) + B' \sin(\omega t)\right) \tag{18}$$

Let's plot these with the same settings as before and $\gamma = 1$, $A'=1$ and $B'=0$:

k = 5
m = 1
gamma = 1

lambd = gamma / (2*m)
omega_d = math.sqrt(k/(lambd*lambd) - 1)

t = np.arange(0, t_s, dt_s)
x = np.exp(- lambd * t) * np.cos(omega_d * t)

fig, axs = plt.subplots(1, 2)
axs[0].axis([0,t_s,-1.5,1.5])
axs[0].set_xlabel("t [s]")
axs[0].set_ylabel("x(t) [m]")
axs[1].axis([-1, 1, -1.5, 1.5])
l1, = axs[0].plot([], [], 'b-')
l2, = axs[1].plot([], [], 'b-')
l3, = axs[1].plot([], [], 'b.', ms=40)

def animate(i):
    l1.set_data(t[:i], x[:i])
    l2.set_data([0, 0], [0, x[i]])
    l3.set_data([0], [x[i]])

ani = matplotlib.animation.FuncAnimation(fig, animate, frames=len(t), interval=int(dt_s*1e3))

plt.close()
HTML(ani.to_html5_video())

One can see that the amplitude does not depend on the frequency, but only on the initial conditions.

(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, γ = 0.1.

Let's assume the solution is of the form:

$$x(t) = A e^{i\omega t} \tag{19}$$

With $A \in \mathbb{C}$ to account for the initial amplitude and phase. Combining this with equation (4.58) yields:

$$-m A \omega^2 e^{i\omega t} + i \gamma A \omega e^{i\omega t} + k A e^{i\omega t} = e^{i\omega t} \tag{19}$$

Which must be true $\forall t$, implying:

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

Isolating $A$ yields:

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

Which is a general complex number. Let's use sympy to find the real and imaginary parts:

import sympy
from IPython.display import display

k, m, omega, gamma = sympy.symbols("k m omega gamma", real=True)
A = 1/(k-m*omega*omega+sympy.I * gamma)
display(sympy.Eq(sympy.S("A"), sympy.expand_complex(A)))

$\displaystyle A = - \frac{i \gamma}{\gamma^{2} + \left(k - m \omega^{2}\right)^{2}} + \frac{k - m \omega^{2}}{\gamma^{2} + \left(k - m \omega^{2}\right)^{2}}$

Finally, let's plot the magnitude and phase as a function of $\omega$ for $m=k=1$, and $\gamma=0.1$:

m = 1
k = 1
gamma = 1

n = 512
omega = np.linspace(0, 3, n)

A = 1/(k-m*omega*omega+gamma*1j)

A_abs = np.abs(A)
A_phase = np.angle(A)

plt.figure()
plt.plot(omega, A_abs)
plt.grid("on")
plt.xlabel("$\omega$")
plt.ylabel("$|A(\omega)|$")

plt.figure()
plt.plot(omega, np.degrees(A_phase))
plt.grid("on")
plt.xlabel("$\omega$")
plt.ylabel("arg($A(\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 defini- tions 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?

TBD

(e) Now find the solution to equation (4.58) by using Laplace transforms. Take the initial condition as $x(0) = ̇\dot{x}(0) = 0$

TBD

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

TBD