Answer 3.1 (f)

_code/pset_02.py

@@ -5,9 +5,10 @@ m = 1
 k = 1
 gamma = 0.1
 
+n_points = 1000
+
 w_min = 0
 w_max = 2
-n_points = 1000
 w = np.linspace(w_min, w_max, n_points)
 A = 1 / (-w**2 * m + w * 1j * gamma + k)
 
@@ -15,14 +16,35 @@ fig1 = plt.figure()
 left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
 ax1 = fig1.add_axes([left, bottom, width, height])
 ax1.plot(w, np.absolute(A))
-ax1.set_ylabel("magnitude")
 ax1.set_xlabel('ω')
+ax1.set_ylabel("magnitude")
 fig1.savefig("../assets/img/02_amplitude.png", transparent=True)
 
 fig2 = plt.figure()
 left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
 ax2 = fig2.add_axes([left, bottom, width, height])
 ax2.plot(w, np.angle(A))
-ax2.set_ylabel("phase")
 ax2.set_xlabel('ω')
+ax2.set_ylabel("phase")
 fig2.savefig("../assets/img/02_phase.png", transparent=True)
+
+t_min = 0
+t_max = 12 * np.pi
+t = np.linspace(t_min, t_max, n_points)
+w = 1.5
+eps = 0.5
+x_0 = np.exp(1j * w * t)
+x_1 = (
+    (1 - eps / (3 * w**3) - eps / (6 * w**2)) * np.exp(1j * w * t) +
+    (eps / (3 * w**3) - eps / (6 * w**2)) * np.exp(-1j * w * t) +
+    (eps / (3 * w**2)) * np.exp(1j * 2 * t)
+)
+fig3 = plt.figure(figsize=[7.4, 4.8])
+left, bottom, width, height = 0.15, 0.1, 0.8, 0.9
+ax3 = fig3.add_axes([left, bottom, width, height])
+ax3.plot(t, np.real(x_0), label="x_0")
+ax3.plot(t, np.real(x_1), label="x_0 + ε x_1")
+ax3.set_xlabel('t')
+ax3.set_ylabel("real part")
+ax3.legend()
+fig3.savefig("../assets/img/02_perturbed_oscillator.png", transparent=True)

_psets/02.md

@@ -160,3 +160,101 @@ $$x(0) = \dot{x}(0) = 0$$.
 {:.question}
 For an arbitrary potential minimum, work out the form of the lowest-order correction to simple
 undamped unforced harmonic motion.
+
+First we add a small quadratic term to the harmonic oscillator ODE.
+
+$$
+m \ddot{x} + kx + \epsilon x^2 = 0
+$$
+
+Though to push around one less symbol, let's define $$\omega = \sqrt{k / m}$$ and absorb a factor of
+$$1/m$$ into $$\epsilon$$.
+
+$$
+\ddot{x} + \omega^2 x + \epsilon x^2 = 0
+$$
+
+We'll express our solution as
+
+$$
+x(t) = x_0(t) + \epsilon x_1(t) + \mathcal{O}(\epsilon^2)
+$$
+
+Plugging this into the ODE above, and dropping $$\mathcal{O}(\epsilon^2)$$ terms, we find that
+
+$$
+\begin{aligned}
+\ddot{x_0} + \epsilon \ddot{x_1} + \omega^2 x_0 + \omega^2 \epsilon x_1 + \epsilon (x_0 + \epsilon x_1)^2 &= 0 \\
+\ddot{x_0} + \epsilon \ddot{x_1} + \omega^2 x_0 + \omega^2 \epsilon x_1 + \epsilon x_0^2 &= 0
+\end{aligned}
+$$
+
+So we recover the original harmonic motion ODE, and gain a second ODE at first order in
+$$\epsilon$$.
+
+$$
+\ddot{x_0} + \omega^2 x_0 = 0 \\
+\ddot{x_1} + \omega^2 x_1 + x_0^2 = 0
+$$
+
+Again using the ansatz $$x_0(t) = e^{rt}$$ for the first equation, we find that $$r^2 = -\omega^2$$.
+So the two roots are $$r = i \omega$$ and $$r = -i \omega$$, and the general solution is
+
+$$
+x_0(t) = A e^{i \omega t} + B e^{-i \omega t}
+$$
+
+To simply algebra I'm going to assume that $$A = 1$$ and $$B = 0$$, so $$x_0(t) = e^{i \omega t}$$.
+This just amounts to choosing some initial conditions: $$x_0(0) = 1$$ and $$\dot{x_0}(0) = i
+\omega$$. The second equation becomes
+
+$$
+\ddot{x_1} + \omega^2 x_1 + e^{i 2 \omega t} = 0
+$$
+
+The general solution is the same as before. For the particular solution, use the ansatz $$x_1(t) = C
+e^{i 2 \omega t}$$.
+
+$$
+-4 \omega^2 C + \omega^2 C + 1 = 0 \\
+C = \frac{1}{3 \omega^2}
+$$
+
+So
+
+$$
+x_1(t) = A e^{i \omega t} + B e^{-i \omega t} + \frac{1}{3 \omega^2} e^{i 2 t}
+$$
+
+Thus our full solution is
+
+$$
+\begin{aligned}
+x(t)
+&= x_0(t) + \epsilon x_1(t) \\
+&= (\epsilon A + 1) e^{i \omega t} + \epsilon B e^{-i \omega t} + \frac{\epsilon}{3 \omega^2} e^{i 2 t}
+\end{aligned}
+$$
+
+To facilitate comparison of $$x$$ and $$x_0$$, let's solve for the same initial conditions we
+imposed earlier.
+
+$$
+\begin{aligned}
+1 &= x(0) \\
+&= \epsilon A + 1 + \epsilon B + \frac{\epsilon}{3 \omega^2} \\
+B &= -A - \frac{1}{3 \omega^2} \\
+i \omega &= \dot{x}(0) \\
+&= i \omega (\epsilon A + 1) - i \omega \epsilon ( 1 - A - \frac{1}{3 \omega^2}) + \frac{2 i \epsilon}{3 \omega^2} \\
+A &= -\frac{1}{6 \omega^2} - \frac{1}{3 \omega^3} \\
+B &= -\frac{1}{6 \omega^2} + \frac{1}{3 \omega^3} \\
+x(t) &= \left( 1 - \frac{\epsilon}{3 \omega^3} - \frac{\epsilon}{6 \omega^2} \right) e^{i \omega t}
+    + \left( \frac{\epsilon}{3 \omega^3} - \frac{\epsilon}{6 \omega^2} \right) e^{-i \omega t}
+    + \frac{\epsilon}{3 \omega^2} e^{i 2 t}
+\end{aligned}
+$$
+
+Now we can compare the two. In this plot $$\omega = 1.5$$ and $$\epsilon = 0.5$$ (which is rather
+large for the model to be quantitatively valid, but qualitatively helps accentuate the difference).
+
+![x_0 vs x_1](../assets/img/02_perturbed_oscillator.png)