Notebook

In [ ]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

Simple Harmonic Motion¶

Lecture 6¶

Examples of oscillatory motion.

First year: the pendulum.

The Simple Pendulum¶

Review of theory...

$\Omega = \sqrt{g/L}$

This is a very simple model.

Simple numerical approach¶

As we did with the both exponential growth/decay, and the projectile motion problems, we'll start with the basic equations and solve them numerically. This way we can compare with the analytical solution from theory. Later, we'll deal with slightly more complicated oscillatory problems for which exact results are not available.

Equations of motion:

\begin{align} \frac{d\omega}{dt} &= - \frac{g}{L} \theta\;, \\ \frac{d\theta}{dt} &= \omega\;, \end{align}

where $\omega$ is the angular velocity of the pendulum. You should recognize the approach of going for a second-order differential equation into a pair of first-order differential equations. We did the same thing when dealing with projectile motion.

We can proceed to solve this problem with Euler's method as we have now down several times already.

\begin{align} \omega_{i+1} &= \omega_i - \frac{g}{L} \theta_i \Delta t \\ \theta_{i+1} &= \theta_i + \omega_i \Delta t \end{align}

In [ ]:

# parameters and constants

L = 1 # m
g = 9.81 # m/s^2

Ω = np.sqrt(g/L)

ω0 = 0.0
θ0 = 0.2

In [ ]:

def calculate_euler(dt=0.04, tmax=10):
    """
    Solve the simple pendulum equations using Euler's method.
    Takes time step dt and maximum time tmax as arguments.
    
    Returns arrays for t, ω, θ
    """
    
    # allocate space
    num_steps = round(tmax/dt)
    t = np.zeros(num_steps)
    ω = np.zeros(num_steps)
    θ = np.zeros(num_steps)

    # initial values
    t[0] = 0
    ω[0] = ω0
    θ[0] = θ0

    # apply Euler's method
    for i in range(num_steps-1):
        t[i+1] = t[i] + dt
        ω[i+1] = ω[i] - g/L*θ[i]*dt
        θ[i+1] = θ[i] + ω[i]*dt
    
    return t, ω, θ

In [ ]:

def plot(t, θ, label=''):
    """
    Plot position vs time
    """
    plt.xlabel('time (s)')
    plt.ylabel(r'$\theta$ (rad)')
    plt.plot(t, θ, linewidth=2, label=label)

In [ ]:

# main driver
fig, axes = plt.subplots()

dt = 0.1
tmax = 10
t, ω, θ = calculate_euler(dt=dt, tmax=tmax)

θ_theory = θ0 * np.cos( Ω*t )
plot(t, θ, label='Numerical, dt={:.2f}'.format(dt))
plot(t, θ_theory, label='Theoretical')
plt.legend(loc='upper left')
plt.show()

Exercise:¶

The exact solution and the numerical answer do not agree. Let's explore this by rerunning the

Change the time step to dt=0.01. Does that improve the problem?

Keeping dt = 0.01, change tmax to be 100.

Turns out that there is no small enough value of dt that will help us in this case. While the motion is necessarily periodic, Euler's method always leads to the solution growing in amplitude with time.

The issue is that the Euler method is an approximation to the differential equation. The key word here is approximate. While for some problems reduce dt can reduce the error in this approximation, for this problem the energy will always grow in time for any nonzero value of $\Delta t$.

Euler's Method of Solving ODE's¶

Notes...

Energy¶

The model doesn't contain a source of energy nor any form of friction so the total energy should remain constant.

For a simple harmonic oscillator, Euler's method an example of a numerical method that is inherently unstable. To understand this instability, we should look at the total energy of the pendulum system.

The energy, $E$, is given by \begin{align} E = \frac{1}{2} m L^2 \omega^2 + m g L ( 1 - cos(\theta) ) \end{align}

The first term is the kinetic energy \begin{align} \frac{1}{2} m v^2 = \frac{1}{2} m (\omega L)^2 \end{align}

The second term is the gravitational potential energy

\begin{align} m g h = m g (L (1-cos\theta) ) \end{align}

where $h$ is the height above the lowest point.

If we assume $\theta$ is small then then the energy can be approximated as \begin{align} E &= \frac{1}{2} m L^2 \omega^2 + m g L \frac{1}{2} \theta^2 \\ &= \frac{1}{2} m L^2 \left( \omega^2 + \frac{g}{L} \theta^2 \right) \end{align}

In [ ]:

def plot_energy(t, ω, θ):
    # calculate the total energy per unit mass
    E = 0.5 * L**2 *(ω**2 + g/L*θ**2)
    
    fig, axes = plt.subplots(2,1, figsize=(8,6))
    
    # plot solution
    plt.sca(axes[0]) # SetCurrentAxes
    plot(t, θ)
    
    # plot energy
    plt.sca(axes[1]) # SetCurrentAxes
    plt.plot(t, E)
    plt.xlabel('t (s)')
    plt.ylabel('E/m (J/kg)')        

In [ ]:

plot_energy(t, ω, θ)
plt.show()

Substitute in Euler's method... and show energy grows with time.

\begin{align} E_{i+1} &= \frac{1}{2} m L^2 \left( \omega^2 + \frac{g}{L} \theta^2 \right) \\ &= \frac{1}{2} m L^2 \left( \left(\omega_i - \frac{g}{L} \theta_i \Delta t \right)^2 + \frac{g}{L} \left(\theta_i + \omega_i \Delta t \right)^2 \right) \\ & = \frac{1}{2} m L^2 \left( \left(\omega_i^2 - 2 \frac{g}{L} \theta_i \omega_i \Delta t + (\frac{g}{L} \theta_i \Delta t)^2\right) + \left( \frac{g}{L}\theta_i^2 + 2 \frac{g}{L}\theta_i \omega_i \Delta t + \frac{g}{L}(\omega_i \Delta t)^2 \right) \right) \\ & = \frac{1}{2} m L^2 \left( \omega_i^2 + (\frac{g}{L} \theta_i \Delta t)^2 + \frac{g}{L}\theta_i^2 + \frac{g}{L}(\omega_i \Delta t)^2 \right) \\ & = \frac{1}{2} m L^2 \left( \omega_i^2 + \frac{g}{L}\theta_i^2 \right) + \frac{1}{2} m L^2 \left( (\frac{g}{L} \theta_i \Delta t)^2 + \frac{g}{L}(\omega_i \Delta t)^2 \right) \\ & = \frac{1}{2} m L^2 \left( \omega_i^2 + \frac{g}{L}\theta_i^2 \right) + \frac{1}{2} m g L \left( \frac{g}{L}\theta_i^2 + \omega_i^2 \right)(\Delta t)^2 \\ E_{i+1} &= E_i + \frac{1}{2} m g L \left( \frac{g}{L}\theta_i^2 + \omega_i^2 \right)(\Delta t)^2 \end{align}

Euler's method itself add energy every time step.

Turns out we can avoid this issue of energy increasing with time by making one small change and using a modified Euler's method: \begin{align} \omega_{i+1} &= \omega_i - \frac{g}{L} \theta_i \Delta t \\ \theta_{i+1} &= \theta_i + \omega_{i+1} \Delta t \end{align}

Exercise¶

Can you spot the difference between this and the original Euler's method?

This modified method is call the Euler-Cromer method. While the value from the previous time step is used to predict the velocity in the next time step, the position is updated using the velocity from the current time step. It is also known as the semi-implicit Euler method.

This one small change gives a significant advantage.

In [ ]:

def calculate_euler_cromer(dt=0.04, tmax=10):
    """
    Solve the simple pendulum equations using Euler-Cromer method.
    Takes time step dt and maximum time tmax as arguments.
    
    Returns arrays for t, ω, θ
    """
    
    # allocate space
    num_steps = round(tmax/dt)
    t = np.zeros(num_steps)
    ω = np.zeros(num_steps)
    θ = np.zeros(num_steps)

    # initial values
    t[0] = 0
    ω[0] = ω0
    θ[0] = θ0

    # apply Euler-Crommer method
    for i in range(num_steps-1):
        t[i+1] = t[i] + dt
        ω[i+1] = ω[i] - g/L*θ[i]*dt
        θ[i+1] = θ[i] + ω[i+1]*dt
    
    return t, ω, θ

In [ ]:

fig, axes = plt.subplots(figsize=(8, 6))

dt = 0.04
tmax = 10
t, ω, θ = calculate_euler_cromer(dt=dt, tmax=tmax)

θ_theory = θ0 * np.cos(Ω*t)
plot(t, θ, label='Numerical, dt={:.2f}'.format(dt))
plot(t, θ_theory, label='Theoretical')
plt.legend(loc='upper left')
plt.show()

In [ ]:

plot_energy(t, ω, θ)
plt.ylim(0, 0.5)
plt.show()

In the Euler-Cromer method, energy (when averaged over one period of an oscillation) is conserved. This energy conservation is a key property of the numerical stability of this method.