We will start by reviewing simple harmonic motion (SHM) as it contains many of the important concepts that we will meet in waves. The most general form of the equation for a simple harmonic oscillator (SHO) including damping and driving forces can be written:

$$ m\frac{\partial^2 \psi}{\partial t^2} = -s\psi -b \frac{\partial \psi}{\partial t} + F_0 \cos \omega t, $$

where $\psi$ is the displacement, $m$ is the mass of the oscillator,
$s$ is the constant giving the restoring force (e.g. a spring
constant) sometimes known as the *stiffness*, $b$ is the damping
coefficient or *resistance* and $F_0$ gives the magnitude of the
driving force (which has *angular* frequency $\omega$).

If we set the damping and driving coefficients to zero, we recover the original, SHM equation:

$$m\frac{\partial^2 \psi}{\partial t^2} = -s \psi,$$

which can be shown to be solved using sinusoidal motion: $$ \psi = A \sin \omega_0 t + B \cos \omega_0 t = C \cos \left(\omega t + \phi\right), $$ where $\omega_0 = \sqrt{s/m}$, $\phi$ is a constant phase and $A = -C\sin \phi$ and $B = C \cos \phi$. As you have seen in both PHAS1245 (Mathematical Methods I) and PHAS1247 (Classical Mechanics), we can use De Moivre's theorem to write the sinusoidal terms as a complex exponential ($e^{i\theta} = \cos \theta + i \sin \theta$ where $i = \sqrt{-1}$), giving: $$ \psi = Re\left[A e^{i(\omega_0 t + \phi)}\right] = Re\left[D e^{i\omega_0 t}\right] $$ where $D = Ae^{i\phi}$.

We often call the argument of a trigonometric function the
*phase*, and $\phi$ is often called the phase difference when
there is more than one oscillation; it represents the offset in the
phase at $t=0$. Now that we have the representation of the
oscillation in terms of the complex exponential (or the sum of a
$\cos$ and $\sin$) we can see that there is an immediate link with
circular motion: the amplitude of the oscillation is just the
projection of circular motion onto the $x$-axis (or any other axis
that is chosen). (Note that from now on I will not write the need to
take the real part of $\psi$ explicitly, but will assume it.) All
these ideas are illustrated in the figure below.

In [4]:

```
from IPython.display import Image
Image(filename='SHMCircularMotion.png',width=300) #Size is in pixels
```

Out[4]:

We can also illustrate this more dynamically, with the python code below. The blue line shows the phasor at $t=0$ (with the angle relative to the x-axis being $\phi$. The red line shows the phasor at a later time $t$, with the real part shown as a dashed red line along the x-axis.

In [16]:

```
phi = 30.0*pi/180. # Phase in degrees
rad = 1.0
x = (0,rad*cos(phi))
y = (0,rad*sin(phi))
axis('scaled')
axis([-1.1,1.1,-1.1,1.1])
plot(x,y)
# Add arc to indicate phi
theta = linspace(0,phi,1000)
plot(0.5*rad*cos(theta),0.5*rad*sin(theta),'b-')
# Label
text(0.5, 0.18, r"$\phi$", fontsize=20, color="blue")
# Now draw full circle
theta = linspace(0,2.0*pi,1000)
plot(rad*cos(theta),rad*sin(theta),'g-')
# Plot the phasor
omega = 2.0
t = 2.5
psix = (0,rad*cos(omega*t+phi))
psiy = (0,rad*sin(omega*t+phi))
plot(psix,psiy,'r-')
# Projection on the x-axis
plot(psix,(0.,0.),'r--')
# Label
text(0.4,-0.2, r"$\psi$", fontsize=20, color="red")
# Arc showing omega t
theta = linspace(phi,omega*t+phi,1000)
plot(0.6*rad*cos(theta),0.6*rad*sin(theta),'r:')
# Label
text(-0.5, 0.6, r"$\omega t$", fontsize=20, color="red")
```

Out[16]:

The representation of a simple harmonic oscillator at a single point in time in the complex plane (e.g. figure above) is often called a *phasor* diagram, and the arrow representing the amplitude and phase of the oscillator relative to a fixed time or phase is called a *phasor*; note that for a phasor we just need
the arrow, not the associated circle. We will see later that we can combine two oscillations using phasors (or using complex arithmetic - the two are completely equivalent).

A simple way to think about phase differences is to consider the velocity and acceleration of the oscillator:

\begin{eqnarray} \psi &=& A e^{i(\omega_0 t + \phi)}\\ \dot{\psi} &=& \frac{\partial \psi}{\partial t} = i\omega_{0} A e^{i(\omega_0 t + \phi)} = i\omega_{0} \psi\\ \ddot{\psi} &=& \frac{\partial^2 \psi}{\partial t^2} = -\omega_{0}^2 A e^{i(\omega_0 t + \phi)} = -\omega_{0}^2 \psi \end{eqnarray}

Notice that both of these quantities vary harmonically and with the same frequency as the oscillation, though with different amplitudes; more importantly, however, note that the velocity is a factor of $i$ different to the displacement, while the acceleration is a factor of $-1$ different. These are easier to understand when we write them in complex exponential form: $i = e^{i\pi/2}$ and $-1 = e^{i\pi}$. So we can write:

\begin{eqnarray} \psi &=& A e^{i(\omega_0 t + \phi)}\\ \dot{\psi} &=& \omega A e^{i(\omega_0 t + \phi + \pi/2)} = \omega_{0} \psi e^{i\pi/2}\\ \ddot{\psi} &=& \omega^2 A e^{i(\omega_0 t + \phi + \pi)} = \omega_{0}^2 \psi e^{i\pi} \end{eqnarray}

We see that the velocity leads the displacement by a phase factor of $\pi/2$ and the acceleration leads the velocity by a phase factor of $\pi/2$. This phase relationship is shown in the next figure.

In [3]:

```
from IPython.display import Image
Image(filename='SHMDispVelAcc.png',width=300) #Size is in pixels
```

Out[3]:

The simple harmonic oscillator, as with all mechanical systems, has two forms of energy: kinetic and potential. If we have $W$ as the total energy and $T$ as kinetic and $V$ as potential, we can write:

\begin{eqnarray} T &=& \frac{1}{2} m\dot{\psi}^2\\ V &=& \frac{1}{2} s\psi^2\\ W &=& T + V \end{eqnarray}

The form of the potential energy follows directly from the form of the force (and is why the motion is called harmonic). As there is no form of dissipation in the motion, we can assert that the total energy does not change with time, just exchanging between potential (at a maximum when the displacement is at a maximum and the velocity is at a minimum) and kinetic (at a maximum when the displacement is at a minimum). They are out of phase - as we expect from the phase differences between displacement and velocity.

As the total energy is constant, we can write:

\begin{eqnarray} \frac{d W}{d t} = \frac{d T}{d t} + \frac{d V}{dt} &=& 0\\ \Rightarrow m\dot{\psi}\ddot{\psi} + s\psi\dot{\psi} &=&0\\ \Rightarrow m\ddot{\psi} &=& -s\psi \end{eqnarray}

which is of course just the original equation for simple harmonic motion.

Restoring the damping term changes the equation and its solution. We will use $\gamma = b/2m$ in the equations below. These solutions were derived in PHAS1247 and are discussed in many textbooks, so I will not derive them here. We find:

\begin{eqnarray} m\frac{\partial^2 \psi}{\partial t^2} &=& -s \psi -b \frac{\partial \psi}{\partial t}\\ \psi &=& \begin{cases} A e^{-\gamma t}e^{i\omega t} &\gamma < \omega_0\\ \hskip 0.5cm\omega = \sqrt{\omega^2_0 - \gamma^2} &\, \\ A e^{-\mu_+ t} + B e^{-\mu_-t} & \gamma > \omega_0\\ \hskip 0.5cm\mu_\pm = \gamma\mp \sqrt{\gamma^2 - \omega_0^2} &\, \\ A(1+\omega_0 t)e^{-\omega_0 t} & \gamma = \omega_0 \end{cases} \end{eqnarray}

The cases for damped harmonic motion correspond to light (or underdamped, $\gamma < \omega_0$), heavy (or overdamped, $\gamma > \omega_0$) and critical ($\gamma=\omega_0$) damping respectively. The total energy is \emph{not} conserved in this case (contrast the undamped oscillator) as the damping force opposes the motion at all times. We write the energy as before, and find the change with time:

\begin{eqnarray} W &=& \frac{1}{2}m\dot{\psi}^2 + \frac{1}{2}s\psi^2\\ \frac{dW}{dt} &=& \frac{d W}{d \dot{\psi}}\frac{d \dot{\psi}}{d t} + \frac{d W}{d \psi}\frac{d \psi}{d t} = (m\ddot{\psi} + s\psi)\dot{\psi}\\ &=& -b\dot{\psi}^2 \end{eqnarray}

where the last line comes from using the SHM equation for damped motion. Notice that the change in energy is always less than zero (unless $\dot{\psi}=0$) and so the total energy of the system decreases with time (as we would expect).

We illustrate this behaviour in the interactive figure below: test the effect of changing `b`

and `k`

(and re-run the cell to render the figure).

In [27]:

```
t = linspace(0,10.0,1000)
m = 1.0
k = 2.0
b = 5.0
omega0 = sqrt(k/m)
gamma = b/(2.0*m)
print "omega0 and gamma are ",omega0,gamma
A = 1.0
if gamma<omega0 :
print "Light damping"
omega = sqrt(omega0*omega0 - gamma*gamma)
plot(t,A*exp(-gamma*t)*cos(omega*t))
elif gamma>omega0 :
print "Heavy damping"
mum = gamma + sqrt(gamma*gamma - omega0*omega0)
mup = gamma - sqrt(gamma*gamma - omega0*omega0)
B = 1.0
plot(t,A*exp(-mum*t) + B*exp(-mup*t),label='sum')
plot(t,A*exp(-mum*t),'r--',label=r"$\mu_{-}$")
plot(t,B*exp(-mup*t),'g--',label=r"$\mu_{+}$")
legend()
```

Now we will introduce a driving term to the system. As we have a harmonic oscillator, it will make sense to use a harmonic driving term (though we could, for instance, use impulses at regular or irregular intervals). Note that the driving term will have an angular frequency $\omega$ which is \emph{different} to the natural frequency of the system, $\omega_0 = \sqrt{s/m}$. Also remember that $\omega = 2\pi \nu$ where $\nu$ is the frequency. We will retain the damping term to give a damped, driven oscillator (again derived in PHAS1247 and textbooks).

\begin{eqnarray} m\frac{\partial^2 \psi}{\partial t^2} &=& -s \psi -b \frac{\partial \psi}{\partial t} + F_0 \cos \omega t\\ \psi &=& A\cos (\omega t + \phi)\\ A &=& \frac{F_0}{m} \left( \frac{1}{(\omega_0^2 - \omega^2)^2 + 4\gamma^2 \omega^2}\right)^{\frac{1}{2}}\\ \tan \phi &=& \frac{-2\gamma\omega}{\omega_0^2 - \omega^2} \end{eqnarray}

The response of the system (including both the phase and the amplitude) depends strongly on the frequency; we can consider three regimes (though there isn't really time to do this in depth).

- When $\omega$ is small (i.e. $\omega \ll \omega_0$), the amplitude can be shown to be $A \simeq F_0/m\omega^2_0 = F_0/s$ and the motion is dominated by the spring constant (
*stiffness controlled*). - When $\omega$ is large (i.e. $\omega \gg \omega_0$), then the amplitude is $A \simeq F_0/m\omega^2$ and the motion is dominated by the mass (
*mass controlled*). - When $\omega \sim \omega_0$, we are at
*resonance*, and the response is controlled by the drag term $\gamma$ (also known*resistance limited*).

The behaviour of $A$ and $\phi$ are plotted below; you should try changing the values of `b`

(damping) and `s`

(stiffness) to see the effects. In particular, you may need to increase or decrease `F0`

to see the graph on a sensible scale. (Remember to run the cell by pressing the play button in the toolbar above.)

In [30]:

```
# Set up parameters of oscillator
m = 1.0
k = 2.0
b = 0.2
gamma = b/(2.0*m)
# Find resonant frequency
omega0 = sqrt(k/m)
print "omega0 and gamma are ",omega0,gamma
F0 = 1.0
# Create array of values of omega from omega0-1 to omega0+1
omega = linspace(omega0-1.0,omega0+1.0,200)
# Calculate phi
phi = arctan(-2.0*gamma*omega/(omega0*omega0 - omega*omega))
# arctan in NumPy is defined to lie between -pi/2 and pi/2 so we correct
# by -pi to give physical reasonableness
phi[100:200] -=pi
# Calculate amplitude
fac = (omega0*omega0 - omega*omega)**2 + 4*gamma*gamma*omega*omega
A = F0/(sqrt(fac)*m)
plot(omega,phi,label='phi')
plot(omega,A,label='A')
legend()
```

Out[30]:

Let's think a little about the power *absorbed* by the oscillator; to do that, we must consider the
work done *against* the drag. If the displacement changes from
$\psi$ to $\psi + \Delta \psi$ then the work done is $-F_d
\Delta\psi$, where $F_d = -b\dot{\psi}$ is the work done against the
drag. If that takes a time $\Delta t$ then the rate of work is $-F_d
(\Delta \psi/\Delta t)$ which tends to $-F_d \dot{\psi}$ as $\Delta t
\rightarrow 0$. So the instantaneous power adsorption becomes:
\begin{eqnarray}
P &=& -F_d \dot{\psi} = b\dot{\psi}^2
\end{eqnarray}

We are not actually interested in the instantaneous power, but the time averaged power (often written $\langle P\rangle$), which for a harmonic force can be shown to be: \begin{eqnarray} \langle P \rangle &=& b\langle \dot{\psi}^2 \rangle = \frac{1}{2}b\omega^2A^2\\ &=& \frac{F_0^2}{m^2}\frac{m\gamma\omega^2}{(\omega_0^2 - \omega^2)^2 + 4\gamma^2\omega^2} \end{eqnarray} where we have substituted $b = 2m\gamma$. This has a maximum value of $F^2_0/(4m\gamma)$ when $\omega = \omega_0$. Note that this is inversely proportional to the resistive force.

Finally we introduce the idea of *impedance* which is a measure
of the resistance to the motion of the oscillator. It is defined as
the amplitude of the driving force divided by the complex amplitude of
the velocity, $Z(\omega) = F_0/\dot{\psi} = b + i(m\omega -
s/\omega)$. At resonance, $Z(\omega) = b$; in the resonance region,
it only departs very slightly from this value. Away from resonance,
it includes a phase lag or lead, which reflects the relation of
velocity to the driving force. We will encounter impedance again with
waves.

We plot below the power absorbed, imaginary part of the impedance and the magnitude of the impedance, using the system given above. Note how the impedance has a minimum at the point that the maximum power is absorbed (and is purely real at that point). Away from that, the power is mainly dissipated (by the impedance).

In [35]:

```
# We use the values of b, m, k, omega and A defined above
power = 0.5*b*omega*omega*A*A
impi = m*omega - k/omega
Z = sqrt(b*b+impi*impi)
plot(omega,power,label='Power')
plot(omega,impi,label='Imaginary part of impedance')
plot(omega,Z,label=r'$\vert Z\vert$')
legend(loc=4)
```

Out[35]:

Note that so far we have discussed a *steady state*
solution: there is a transient behaviour at the beginning of the
oscillation which takes the form of a damped oscillator
added to a solution of the homogenous equation. So a real system will
respond at its natural frequency, but damping will cause that solution
to die off, while the driven oscillation persists into the steady
state:

\begin{equation}
\psi(t) = A_s \cos(\omega t +\phi) + A_t e^{-\gamma t} \cos(\omega_0 t + \phi_t)
\end{equation}
where the subscripts s and t stand for *steady state* and
*transient* respectively. For lightly damped system near
resonance we can approximate and find that $A_t \simeq A_s$ and
$\phi_t \simeq \phi_s + \pi$, giving:
\begin{equation}
\psi(t) = A_s \left(\cos(\omega t +\phi) - e^{-\gamma t} \cos(\omega_0 t + \phi_s)\right)
\end{equation}

We plot an example of this below, with the total response of the system over time in blue, and the transient in green: notice how the initial response of the system is at its natural frequency, $\omega_0$, and how this dies down and leads to the steady-state, driven response at $\omega$. You should try changing the value of the driving frequency `omegaD`

(though you will also have to change the amplitude of the driving force, `F0`

, and the initial response, `F1`

).

In [44]:

```
t = linspace(0,50.0,1000)
# Driving frequency
omegaD = 7.0
# Amplitude of driving (F0) and response (F1)
F0 = 10.0
F1 = 0.10
fac = (omega0*omega0 - omegaD*omegaD)**2 + 4*gamma*gamma*omegaD*omegaD
A = F0/(sqrt(fac)*m)
# Define transient
transient = 2*(F1/m)*exp(-gamma*t)*cos(omega0*t+pi/2)
steady = A*cos(omegaD*t+pi/2)
plot(t,transient+steady)
plot(t,transient)
#plot(t,A*cos(omegaD*t+pi/2)+2*(F1/m)*exp(-gamma*t)*cos(omega0*t+pi/2))
#plot(t,2*(F1/m)*exp(-gamma*t)*cos(omega0*t+pi/2))
```

Out[44]:

Instead of combining a single driving force with a damped response,
what happens if we have an oscillator with *two* driving forces ?
The full equation looks like this:

$$ \begin{equation} m\frac{\partial^2 \psi}{\partial t^2} = -s \psi -b \frac{\partial \psi}{\partial t} + F_1 \cos (\omega_1 t + \phi_1) + F_2 \cos (\omega_2 t + \phi_2), \end{equation} $$

We know that the SHO equation is *linear*, so we can try solving
this problem by summing the steady state solutions arising from each
driving force independently. We will consider three specific cases
(where it is easier to understand the behaviour) with the general case
more complex:

- Same amplitude and frequency ($F_1 = F_2 = F$, $\omega_1 = \omega_2 = \omega$), different phase ($\phi_1 \ne \phi_2$)
- Same frequency ($\omega_1 = \omega_2$), different phase and amplitude ($F_1 \ne F_2, \phi_1 \ne \phi_2$)
- Same amplitude and phase ($F_1 = F_2=F, \phi_1 = \phi_2$), different frequency ($\omega_1 \ne \omega_2$)

In all cases, we assume that we can write the displacement as a superposition of the responses to the individual oscillations: $$ \begin{equation} \psi = A_1 \cos(\omega_1 t + \phi_1) + A_2 \cos(\omega_2 t + \phi_2) \end{equation} $$

For simplicity, we will set
$\phi_1 = 0$ and set $\phi_2 = \phi$. Then we can perform the
following manipulations, using the complex exponential form for
simplicity:
$$
\begin{eqnarray}
\psi &=& Ae^{i\omega t} + Ae^{i(\omega t + \phi)}\\
&=& Ae^{i\omega t}\left(1 + e^{i\phi}\right)\\
&=& Ae^{i\omega t}\left(e^{i\phi/2}e^{-i\phi/2} + e^{i\phi/2}e^{i\phi/2}\right)\\
&=& Ae^{i\omega t}e^{i\phi/2}\left(e^{-i\phi/2} + e^{i\phi/2}\right)\\
&=& 2Ae^{i(\omega t+\phi/2)}\cos(\phi/2) = 2A\cos(\phi/2)e^{i(\omega t+\phi/2)}
\end{eqnarray}
$$
where the identification of $\cos (\phi/2)$ follows from De Moivre's
theorem. So the *resultant* oscillation has magnitude
$2A\cos(\phi/2)$ and phase $\phi/2$. The same result can be found
using trigonometry and a phasor diagram

In [28]:

```
from IPython.display import Image
Image(filename='PhasorSameFreqAmp.png',width=200)
```

Out[28]:

We can write a solution as before, by using superposition:

$$ \begin{equation} \psi = A_1 \cos(\omega t + \phi_1) + A_2 \cos(\omega t + \phi_2) = A \cos(\omega t + \theta) \end{equation} $$ This is illustrated in a phasor diagram below.

In [27]:

```
from IPython.display import Image
Image(filename='Phasors.png',width=200)
```

Out[27]:

This is just an oscillation with frequency $\omega$ but a total amplitude and phase which depend on the amplitudes and phases of the two driving forces. If, for instance, the phase difference is $\phi_1

- \phi_2 = \pi$ then the two driving forces are out-of-phase and we
will get destructive interference. The use of phasors allows a simple
visualisation of the resultant. Using trigonometry from the phasor
diagram (or the cosine rule), we can write the amplitude of the
resultant as:
$$
\begin{equation}
A^2 = A_1^2 + A^2_2 + 2A_1A_2\cos(\phi_2 - \phi_1)
\end{equation}
$$
Now we can see that the amplitude will vary between $A_1+A_2$ if
$\phi_1=\phi_2$ and $|A_1 - A_2|$ if $|\phi_1 - \phi_2| = \pi$. If
the driving forces are
*in phase*then we have the maximum amplitude, while if they are in*antiphase*we have the minimum amplitude.

The phase is also found using trigonometry; by projecting onto thereal and imaginary axes, we can write: $$ \begin{equation} \tan \theta = \frac{A_1\sin \phi_1 + A_2\sin \phi_2}{A_1 \cos \phi_1 + A_2 \cos \phi_2} \end{equation} $$ We can get the same result using complex notation. Start by finding that, for two general complex numbers $z_1$ and $z_2$: $$ \begin{eqnarray} \left| z_1 + z_2\right|^2 &=& (z_1+z_2)(z_1 + z_2)^\star = (z_1+z_2)(z^\star_1+z^\star_2)\\ &=& |z_1|^2 + |z_2|^2 + (z_1z^\star_2 + z^\star_1z_2)\\ &=& |z_1|^2 + |z_2|^2 +2\mathrm{Re}(z_1z^\star_2) \end{eqnarray} $$ Now we have $z_1 = A_1 e^{i(\omega t +\phi_1)}$ and $z_2 = A_2 e^{i(\omega t +\phi_2)}$ with $A_1$ and $A_2$ real. This gives: $$ \begin{eqnarray} \left| z_1 + z_2\right|^2 &=& |z_1|^2 + |z_2|^2 +2\Re(z_1z^\star_2)\\ &=& A_1^2 + A_2^2 + 2\mathrm{Re}(A_1 e^{i(\omega t +\phi_1)}A_2 e^{-i(\omega t +\phi_2)})\\ &=& A_1^2 + A_2^2 + 2\mathrm{Re}(A_1A_2 e^{i(\phi_1 - \phi_2)})\\ &=& A_1^2 + A_2^2 + 2A_1A_2 \cos(\phi_2 - \phi_1) \end{eqnarray} $$ while we can find the phase from: $$ \begin{eqnarray} \mathrm{arg}(z_1+z_2) &=& \tan^{-1} \left[ \frac{\mathrm{Im}(z_1+z_2)}{\mathrm{Re}(z_1+z_2)} \right]\\ \Rightarrow \theta &=& \tan^{-1} \left[\frac{A_1\sin \phi_1 + A_2\sin \phi_2}{A_1 \cos \phi_1 + A_2 \cos \phi_2}\right] \end{eqnarray} $$

For this case, we use the principle of superposition to write: $$ \begin{equation} \psi = A \cos \omega_1 t + A \cos \omega_2 t \end{equation} $$ Note that we have assumed that the common phase can be set to zero for simplicity. But we can rearrange this using a standard trigonometrical formula for the sum of two cosines: $$ \begin{eqnarray} \psi &=& A \left(\cos \omega_1 t + \cos \omega_2 t\right)\\ &=& 2A \cos\left(\frac{\omega_1 + \omega_2}{2} t\right) \cos\left(\frac{\omega_1 - \omega_2}{2} t\right) \end{eqnarray} $$

Now let's define two \emph{new} angular frequencies, and rewrite the solution: $$ \begin{eqnarray} \omega &=& \frac{1}{2}\left(\omega_1 + \omega_2\right)\\ \Delta \omega &=& \frac{1}{2}\left(\omega_1 - \omega_2\right)\\ \psi &=& 2A\cos (\omega t) \cos (\Delta\omega t) \end{eqnarray} $$ If the two angular frequencies are relatively close, then the form of the resulting oscillation is rather simple. There is an oscillation at the average frequency (the term $\cos \omega t$) whose amplitude is \emph{modulated} by a slow oscillation at the difference frequency (the term $\cos \Delta\omega t$ - also known as the envelope). This is a phenomenon known as \emph{beats}, and is illustrated in the figure below. It is important to understand that, while the angular frequency of the modulation is $\Delta \omega = \frac{1}{2}|\omega_1 - \omega_2|$, the angular frequency at which \emph{peaks} of activity occur is $|\omega_1 - \omega_2|$ (or equivalently zeroes of activity). The number of minima per second is $\Delta \omega/\pi$. The perceived effect (say for sound) will be of a sound at the average angular frequency with its amplitude varying according to the envelope.

In [48]:

```
fig = figure(figsize=[12,3])
# Add subplots: number in y, x, index number
ax = fig.add_subplot(121,autoscale_on=False,xlim=(0,65),ylim=(-2,2))
ax2 = fig.add_subplot(122,autoscale_on=False,xlim=(0,65),ylim=(-2,2))
t = linspace(0.0,65.0,1000)
o1 = 1.3
o2 = 1.0
omega = 0.5*(o1+o2)
domega = 0.5*(o1-o2)
# Plot the resulting pattern
ax.plot(t,2.0*cos(omega*t)*cos(domega*t))
# Show the envelope in a dashed red line - we need two lines
ax.plot(t,2.0*cos(domega*t),'r--')
ax.plot(t,2.0*cos(domega*t+pi),'r--')
# On a separate plot, show the two individual oscillations
ax2.plot(t,cos(o1*t))
ax2.plot(t,cos(o2*t))
```

Out[48]:

The result of two driving forces on one SHO with different frequencies but the same amplitude. Left: two forces shown together. Right: the resultant motion with the envelope superimposed in dashed lines. The original frequencies are $\omega = 1.0$ and $\omega=1.2 s^{-1}$; on ce how the peaks coincide at regular intervals, and this leads to the maxima in the envelope function.

The figure above shows an illustration of exactly this behaviour for the frequencies $\omega_1 = 1.2 s^{-1}$ and $\omega_2 = 1.0 s^{-1}$ which gives a resulting motion at frequency $\omega = 1.1 s^{-1}$ modulated by an envelope with frequency $\Delta \omega = 0.1 s^{-1}$.

The resulting motion will show a true periodicity if the ratio of the frequencies can be written as a ratio of integers (i.e. $\omega_1/\omega_2 = n_1/n_2$). The motion from two frequencies which are almost non-periodic is shown in the figure below (I used 25/3 and 7.1 here but we could have made it properly non-periodic if we'd tried harder).

Note that if there is a phase difference between the two oscillations then beats still result, but the phase difference will appear in both sum and difference phases. We can write $\cos(\omega_1 t + \phi_1) + \cos(\omega_2 t + \phi_2) = 2\cos(\omega t + \phi)\cos(\Delta \omega t + \Delta \phi)$ where $\omega$ and $\Delta \omega$ have the same meaning as before and $\phi = (\phi_1 + \phi_2)/2$ and $\Delta \phi = (\phi_1 - \phi_2)/2$.

In [24]:

```
fig = figure(figsize=[12,3])
# Add subplots: number in y, x, index number
xmax = 15.0
ax = fig.add_subplot(121,autoscale_on=False,xlim=(0,xmax),ylim=(-2,2))
ax2 = fig.add_subplot(122,autoscale_on=False,xlim=(0,xmax),ylim=(-2,2))
t = linspace(0.0,xmax,1000)
o1 = 25./3
o2 = 7.1
omega = 0.5*(o1+o2)
domega = 0.5*(o1-o2)
# Plot the resulting pattern
ax.plot(t,2.0*cos(omega*t)*cos(domega*t))
# Show the envelope in a dashed red line - we need two lines
ax.plot(t,2.0*cos(domega*t),'r--')
ax.plot(t,2.0*cos(domega*t+pi),'r--')
# On a separate plot, show the two individual oscillations
ax2.plot(t,cos(o1*t))
ax2.plot(t,cos(o2*t))
```

Out[24]:

If the amplitudes and frequencies all differ, then the phasor diagram must be dynamic as illustrated below.

In [21]:

```
from IPython.display import Image
Image(filename='PhasorsDiffFreq2.png',width=200) #Size is in pixels
```

Out[21]:

The
tip of the resultant will trace out a shape in time called a
*cycloid*. For instance, the path followed by the tip of the
resultant vector in the complex plane for different amplitudes and
frequencies differing by a factor of two is plotted below - you should try playing with the parameters, and seeing what effects they have.

In [15]:

```
# Define figure size to get square
fig = figure(figsize=[5,5])
# The range of t may need to change with angular frequencies
t = linspace(0,10*pi,200)
# Original frequencies: 1.2 and 0.6
o1 = 1.2
o2 = 0.6
# Original amplitudes: 1.5 and 1.3
A1 = 1.5
A2 = 1.3
# Calculate x and y projections of the two oscillations when added
x = A1*cos(o1*t) + A2*cos(o2*t)
y = A1*sin(o1*t) + A2*sin(o2*t)
# Plot path of the tip of the vector
plot(x,y)
```

Out[15]: