# Description Objects which have a periodic behaviour described as _simple harmonic motion_ are called "Harmonic oscillators". This category includes a wide variety of contraptions including springs, pendulums etc. ## Simple harmonic motion Simple harmonic motion is motion following the form $m\ddot{x} + kx = k x_{eq},$ meaning that there is no "forcing" (terms which add additional energy to the oscillation over time) or "damping" (things such as friction, shock absorbers etc which would decrease the oscillations over time). For reasons which will hopefully become apparent this is often rewritten with $\omega^2 = k/m$ as $\ddot{x} + \omega^2 x = \omega^2 x_{eq} \tag{1}$ This has a general solution (which we will derive below) of the form $x(t) = B \cos \omega t + C \sin \omega t + x_{eq}, \tag{2}$ If given a set of initial conditions, it will usually be easiest to find a particular solution using form (2), however a slightly more useful equivalent version is $x(t) = A \cos(\omega t + \phi) + x_{eq} \tag{3}$ where: - $x$ is position, - $t$ is time, - $\omega$ is the _angular frequency_ (determining period and frequency), - $-\pi < \phi \leq \pi$ is the _phase angle_ of the oscillations - $x_{eq}$ is the _equilibrium position_, and - $A$, $B$ and $C$ are arbitrary (but related) constants, with $A>0$. Note that the following relationships hold - $A = \sqrt{B^2 + C^2}$ is the amplitude of the oscillation, - $\cos \phi = B/A$, and $\sin \phi = -C/A$, - the _period_(time for a full cycle) of the oscillation is $\tau = 2\pi/\omega$, - the _frequency_(number of cycles per second) is $f=1/\tau=\omega/2\pi$. The reason (3) is more useful is that the amplitude is apparent and also that since the phase angle is explicit and there is only a single term, the expression is much easier to sketch and reason about. If you have access to mathematica, here is a Manipulator that will allow you to play around with this functional form and see what all the parameters do. ```mathematica Manipulate[ Plot[ a Cos[\[Omega] t + \[Phi]] + x0, {t, 0, 2 \[Pi]}, PlotLabel -> "x(t)= a Cos(\[Omega]t + \[Phi]) + \!\(\*SubscriptBox[\(x\), \ \(0\)]\)", AxesLabel -> {t, "x(t)"} ], {{\[Omega], 1, "angular frequency (\[Omega])"}, -10, 10, Appearance -> "Labeled"}, {{a, 1, "amplitude (a)"}, -10, 10, Appearance -> "Labeled"}, {{\[Phi], 0, "phase angle (\[Phi])"}, -\[Pi], \[Pi], Appearance -> "Labeled"}, {{x0, 0, "equilibrium position (\!\(\*SubscriptBox[\(x\), \(0\)]\))"}, -3, 3, Appearance -> "Labeled"}] ``` ## Derivation Consider the motion of a mass $m$ on a sled on a smooth surface connected to a solid wall by a spring of stiffness $k$ and natural length $l_0$. The spring is extended (or compressed - it doesn't matter which) to length $l$ and then released with no initial velocity. *Note:* This isn't the only setup which could be used to derive these equations, but this is convenient in that the restorative force of the spring and the force due to weight are perpendicular to each other. Let - $x(t)$ be the position at time $t$ with the wall as the origin and $x$ increasing away from the wall, - $\mathbf{i}$ be parallel with the surface away from the wall and $\mathbf{j}$ perpendicular to $\mathbf{i}$ upwards, and - $\widehat{\mathbf{s}}$ be a unit vector in the direction from the mass to the spring. Then the setup is as follows: ![[harmonic oscillator 1.png]] Now let - $\mathbf{W}$ be the force due to weight, with $|\mathbf{W}|=m g$ where $g$ is the constant of gravitational acceleration, - $\mathbf{N}$ be the normal reaction of the surface, with $|\mathbf{N}|=N$, and - $\mathbf{H}$ be the restorative force of the spring. Model the mass as a particle and the spring as a model spring (in the sense that it doesn't deform and the force follows Hooke's law precisely) and given the surface is smooth, assume no friction. Assume also that all motion occurs in two dimensions and there is no wobble or lateral motion introduced by the spring. Then a diagram representing the forces at work in this situation is ![[harmonic oscillator 2.png]] Note that this is showing the case where the spring has been extended so the restorative force is a contraction towards the equilibrium position. $\mathbf{H}$ could equally be in the positive $\mathbf{i}$ direction if the spring had been compressed and was therefore exerting an expansive force. This would reverse the direction of the arrow for $\mathbf{H}$ in the diagram above, but the arrow for $\widehat{\mathbf{s}}$ would not change. It is conventional (I believe) to draw $\mathbf{H}$ in the same direction as $\widehat{\mathbf{s}}$ in the way shown above irrespective of the fact that the actual force may be antiparallel to this. The resultant force is thus $\mathbf{R}=\mathbf{H}+\mathbf{N}+\mathbf{W}$ with $ \begin{align*} \widehat{\mathbf{s}} &= -\mathbf{i}\\ \mathbf{W} &= -mg\mathbf{j}\\ \mathbf{N} &= N\mathbf{j}\\ \mathbf{H} &= k(l-l_0)\widehat{\mathbf{s}} &~\text{by Hooke's Law}\\ &= -k(l-l_0)\mathbf{i}\\[8px] \text{So}~\mathbf{R} &= (N-mg)\mathbf{j} -k(l-l_0)\mathbf{i}. \end{align*} $ See some context about [[#Sidebar on Hooke's Law|Hooke's Law]] below. Resolving first in the $\mathbf{j}$ direction, since the sled does not leave the surface, vertical acceleration must be 0 and hence $ \begin{align*} R_j &= 0\\ N-mg &= 0\\ N &= mg. \end{align*} $ Then, in the $\mathbf{i}$ direction, $ \begin{align*} R_i &= -k(l-l_0). \end{align*} $ So since there is no vertical movement, and we have assumed no lateral movement, we can simplify things and consider only horizontal motion from now on. From Newton's second law $ \begin{align*} R &= ma\\ \text{and thus}~a &= \frac{R}{m},\\ \text{so}~a &= -\frac{k}{m}(l-l_0). \end{align*} $ Now notice that $a = \ddot{x}$ (in [[Newton's Laws of Motion#Newton’s notation for derivatives|Newton's notation]]) and $l = x$, so $ \begin{align*} \ddot{x} &= -\frac{k}{m}(x-l_0)\\ \ddot{x} + \frac{k}{m} x &= \frac{k}{m}l_0\\ \text{Let}~\omega^2&=\frac{k}{m}\\ \text{so that}~ \ddot{x} + \omega^2 x &= \omega^2l_0 \tag{4} \end{align*} $ ### Derivation of (2) Now since this is an inhomogeneous constant-coefficient linear second-order differential equation, to find the solution the [[Solving ODEs#Inhomogeneous linear constant-coefficient second-order ODE|method]] is to first find the complementary function $x_c(t)$, then find a particular integral $x_p(t)$, and then use superposition to claim that $x(t) = x_c(t) + x_p(t)$ and that this is the general solution to this differential equation. To find $x_c$, which is a general solution to the associated homogeneous differential equation $\ddot{x} + \omega^2 x=0$, consider the auxiliary equation $ \begin{align*} \lambda^2+\omega^2 &= 0\\ \text{thus}~\lambda^2 &= -\omega^2\\ \lambda &= \pm \omega \;\mathrm{i} \end{align*} $ Thus $ \begin{align*} x_c(t) &= \mathrm{e}^{\alpha t}(B \cos\beta t + C \sin \beta t) \tag{5}\\ &= B \cos\omega t + C \sin \omega t\\ \end{align*} $ …since we have $\alpha=0$ and $\beta=\omega$. (See the [[#Derivation of (5)]] below) Then, to find $x_p$ since it is constant with respect to time, guess $x_p(t) = c$ so that $\dot{x}_p=\ddot{x}_p=0$, and (4) becomes $ \begin{align*} 0 + \omega^2x_p &= \omega^2l_0\\ \text{and thus}~x_p &= l_0 \end{align*} $ This is known as the _equilibrium position_, and often denoted $x_{eq}$. In our case it is simply the natural length of the spring, but in more complicated setups it will be the position in which the oscillator has not been displaced from rest so there is no restorative force and hence no acceleration. Thus, by superposition $ \begin{align*} x(t) &= x_c(t)+x_p(t)\\ &= B \cos\omega t + C \sin \omega t + x_{eq} \end{align*} $ ...where $B$ and $C$ are arbitrary constants, and we have derived (2). Note also this gives rise to an alternative formulation for the acceleration equation for simple harmonic motion, which is $m\ddot{x} +kx = kx_{eq}.$ It can be seen that this is entirely equivalent to (1) when $\omega^2 = k/m$. #### Check $ \begin{align*} \text{If}~x(t) &= B \cos\omega t + C \sin \omega t + x_{eq}\\ \text{then}~\dot{x}(t) &= -B \omega \sin \omega t + C \omega \cos \omega t\\ \text{and}~\ddot{x}(t) &= -B \omega^2 \cos \omega t - C \omega^2 \sin \omega t\\ LHS &= \ddot{x} + \omega^2 x \\ &=-B \omega^2 \cos \omega t - C \omega^2 \sin \omega t + \omega^2\bigl(B \cos\omega t + C \sin \omega t + x_{eq}\bigr)\\ &= \omega^2 x_{eq} = RHS\\ \end{align*} $ So $x(t) = B \cos\omega t + C \sin \omega t + x_{eq}$ is indeed a general solution to $\ddot{x} + \omega^2 x= \omega^2 x_{eq}$. ### Derivation of (3) Then, to derive (3), consider $ \begin{align*} A \cos(\omega t + \phi) &= A \cos \omega t \cos \phi - A \sin \omega t \sin \phi &\text{$\cos$ angle addition identity}\\ \text{so if}~A\cos(\omega t + \phi) + x_{eq}&= B\cos \omega t + C \sin \omega t + x_{eq},\\ \text{then}~A\cos \omega t \cos \phi - A \sin \omega t \sin \phi &= B \cos \omega t + C \sin \omega t\\ \end{align*} $ Considering first terms in $\cos \omega t$, $ \begin{align*} A \cos \omega t \cos \phi &= B \cos \omega t\\ \cos \phi &= \frac{B}{A} &~\text{dividing by $A\cos \omega t$} \end{align*} $ Considering then terms in $\sin \omega t$, $ \begin{align*} -A \sin \omega t\sin \phi &= C \sin \omega t\\ \sin \phi &= -\frac{C}{A} &~\text{dividing by $-A \sin \omega t$} \end{align*} $ Therefore, $ \begin{align*} B^2 + C^2 &= (A \cos \phi)^2 + (-A \sin \phi)^2\\ &= A^2 \cos^2 \phi + (-A)^2 \sin^2 \phi\\ &= A^2 (\cos^2 \phi + \sin^2 \phi)\\ &= A^2 &~\text{by the Pythagorean trigonometric identity}\\ \text{So}~A &= \sqrt{B^2 + C^2}. \end{align*} $ And hence, so long as $\cos \phi = B/A$ and $\sin \phi = -C/A$ an alternative form (3) is $x(t) = A \cos(\omega t + \phi) + x_{eq}$ where - $A=\sqrt{B^2 + C^2}$, - $B$ and $C$ are arbitrary constants, - $\omega = \sqrt{\frac{k}{m}}$ is the angular frequency, - $\phi$ is the phase angle, and - $x_{eq}$ is the equilibrium position *Note*: when solving problems involving harmonic oscillators, to find the correct value for $\phi$ it is generally necessary to first consider their signs of _both_ $\cos \phi$ and $\sin \phi$ and the resulting quadrant $\phi$ needs to be in before applying the requisite inverse trigonometric function(s). If you think you can just do $\phi = \arccos \frac{B}{A}$ you'll often find you've got the sign wrong and/or you should have added $\pi$. ### Derivation of (5) If we have complex conjugate eigenvalues of the form $\lambda = \alpha \pm \beta \;\mathrm{i}$, the general solution is worth deriving for that. We know that if we have roots of the auxiliary equation $\lambda_1$ and $\lambda_2$, then the general form of the solution is $x(t) = D\mathrm{e}^{\lambda_1 t} + E\mathrm{e}^{\lambda_2 t} $ ...for some complex constants $D$ and $E$. So $\begin{align*} x(t) &= D\mathrm{e}^{(\alpha+\beta\mathrm{i}) t} + E\mathrm{e}^{(\alpha-\beta\mathrm{i}) t} \\ &= D \mathrm{e}^{\alpha t} \mathrm{e}^{\beta\mathrm{i} t} + E\mathrm{e}^{\alpha t} \mathrm{e}^{-\beta\mathrm{i} t} \\ \end{align*} $ Now recall [[Complex Numbers#...using Taylor series|Euler's formula]] $\mathrm{e}^{\mathrm{i}\theta} = \cos \theta + \mathrm{i} \sin \theta$. This means $\begin{align*} x(t) &= D \mathrm{e}^{\alpha t} (\cos \beta t + \mathrm{i}\sin \beta t) + E\mathrm{e}^{\alpha t} \bigl(\cos (-\beta t) + \mathrm{i}\sin (-\beta t)\bigr) \end{align*} $ But $\cos -\theta = \cos \theta$ and $\sin -\theta = -\sin \theta$, so $\begin{align*} x(t) &= D \mathrm{e}^{\alpha t} (\cos \beta t + \mathrm{i}\sin \beta t) + E\mathrm{e}^{\alpha t} \bigl(\cos \beta t - \mathrm{i}\sin \beta t\bigr)\\ &= D \mathrm{e}^{\alpha t} \cos \beta t + \mathrm{i} D \mathrm{e}^{\alpha t}\sin \beta t + E\mathrm{e}^{\alpha t} \cos \beta t - \mathrm{i} E \mathrm{e}^{\alpha t} \sin \beta t\\ &= \mathrm{e}^{\alpha t} \bigl((D + E)\cos \beta t+\mathrm{i}(D-E)\sin \beta t\bigr)\\ \end{align*} $ So let $B=D+E$, and $C=\mathrm{i}(D-E)$ and we have a general solution. $ x(t)= \mathrm{e}^{\alpha t} \bigl(B\cos \beta t+C\sin \beta t\bigr)\\ $ Provided any initial conditions are real-valued, $D$ and $E$ are complex variables such that $B$ and $C$ are both real. ### External resources Here is Trefor Bazett's explanation of the harmonic oscillator using a very similar setup to the one above  ...and here is Steve Brunton solving the Harmonic Ocillator 4 ways, which is pretty instructive in general.  ### Sidebar on Hooke's Law Notice that [Hooke's Law](https://phys.libretexts.org/Courses/University_of_California_Davis/UCD%3A_Physics_7B_-_General_Physics/6%3A_Newton's_Laws_of_Motion/6.2%3A_The_Force_model) , which we have used here to get the force in the spring, is less of an actual law than an empirical observation relating spring force to displacement in a way that is compatible with Newton's Laws. It is only ever true over certain portions of the displacement of a real spring, depending on materials and other factors. For example if you take a weak spring and extend it too far, or repeatedly extend and compress a real spring then in actual practise it will deform and not spring back to the natural length it had previously. We work around these limitations by explicitly assuming a model spring and recognising this will in practise limit the applicability of our model to only certain displacements, and make them inaccurate in the face of very weak materials or repetitive strain from frequent use of a sprung gadget in real life. Additionally, you will often (as in both videos linked above) see the 1-D version of Hooke's Law written as $H = -kx.$ When giving the law in this form, the coordinates system has been chosen such that $x=0$ is the equilibrium position, meaning that $x$ gives the displacement from the equilibrium position given as $(l-l_0)$ in the formulation used above. In this system $x$ also gets larger as the spring extends, so the restorative force for positive displacements is in the negative $x$ direction. The formulation I have chosen $\mathbf{H} = k(l-l_0)\widehat{\mathbf{s}}$ makes explicit the direction of the force (by the use of $\widehat{\mathbf{s}}$) and allows for an arbitrary choice of coordinates at the expense of being a little more complex-looking for the most basic case. # Solving oscillators using a CAS Computer algebra systems such as maxima or Mathematica are often a little bit persnickety when it comes to solving oscillators and they need to be done in a specific way if you want consistent results. In wxMaxima solving a pendulum oscillator for example with a given length looks like this: ```maxima kill(solns,a,phi, B, C)$ params:[h=2.5, l=2.1, g=9.81]$ derived:at([ω=sqrt(g/l),τ=(2*%pi*sqrt(l))/sqrt(g)],params)$ params:append(params,derived)$ /* it's easier to find a particular solution using this form */ altθ:B*cos(ω*t)+C*sin(ω*t)$ daltθdt:diff(altθ,t)$ /*initial conditions */ ics:at([altθ =%pi/4,daltθdt=0],append([t=0],params)); solns:linsolve(ics,[B,C]); /* having found B and C, find the appropriate A and φ */ [a,phi]:at([sqrt(B^2+C^2),asin(-C/A)],solns)$ params:append(params,[A=a,φ=phi]); ``` Notice the pattern: 1) Create a list of parameters. This helps to prevent your local context to be too stuffed with assignments 2) Calculate derived parameters such as $\tau$ and $\omega$ and add them to the list 3) Solve the initial conditions using the form given at (2) above and then calculate $A$ and $\phi$ to put it into form (3), appending those to the parameters so they are available for use by other parts of whatever model. - It's helpful (as above) to start the cell by `kill`-ing any variables you're going to derive in the cell so `solve`/`linsolve` doesn't fail if you run the block multiple times while working on it. - Although I suppress a lot of output, I generally try to always print out the results of `solve`/`linsolve` to give immediate visual feedback if something is amiss because I will just see an empty list instead of my solutions. Here's that same example in mathematica ```mathematica Clear[\[Alpha], \[Beta], \[Gamma], a, b, c, \[Phi]]; (* pendulum motion *) \[Theta]0[ t_] := \[Beta] Cos[\[Omega] t ] + \[Gamma] Sin[\[Omega] t ]; ics = {\[Theta]0[0] == \[Pi]/4, \[Theta]0'[0] == 0}; (* solve for \[Beta] and \[Gamma] *) sols = First[g Assuming[\[Beta] > 0 && g > 0 && l > 0, Solve[ics /. \[Omega] -> Sqrt[g/l], {\[Beta], \[Gamma]}, Reals]]]; (* use that to find \[Alpha] and \[Phi], so we have \[Theta][t] in \ the standard form *) \[Alpha] = Sqrt[\[Beta]^2 + \[Gamma]^2] /. sols; \[Phi] = ArcSin[-(\[Gamma]/\[Alpha])] /. sols; \[Theta][t_] := \[Alpha] Cos[\[Omega] t + \[Phi]]; \[Theta][t] ``` - Notice how the `Assuming` ... `Solve[... , Reals]` works. That's a very useful general pattern for hinting `Solve`, `Reduce` etc by putting bounds on variables. - I don't try to name things using subscripts. I don't understand how that is supposed to work in mathematica but every time I try that it just breaks things.