Note

This notebook can be downloaded here: 04_ODE_Harmonic_Oscillator.ipynb

Tutorial 2: Driven Harmonic Oscillator

In this example, you will simulate an harmonic oscillator and compare the numerical solution to the closed form one.

Theory

Read about the theory of harmonic oscillators on Wikipedia

Mechanical oscillator

The case of the one dimensional mechanical oscillator leads to the following equation:

\[m \ddot x + \mu \dot x + k x = m \ddot x_d\]

Where:

• \(x\) is the position,
• \(\dot x\) and \(\ddot x\) are respectively the speed and acceleration,
• \(m\) is the mass,
• \(\mu\) the
• \(k\) the stiffness,
• and \(\ddot x_d\) the driving acceleration which is null if the oscillator is free.

Canonical equation

Most 1D oscilators follow the same canonical equation:

\[\ddot x + 2 \zeta \omega_0 \dot x + \omega_0^2 x = \ddot x_d\]

Where:

• \(\omega_0\) is the undamped pulsation,
• \(\zeta\) is damping ratio,
• \(\ddot x_d = a_d\sin(\omega_d t)\) is the imposed acceleration.

In the case of the mechanical oscillator:

\[\omega_0 = \sqrt{\dfrac{k}{m}}\]

\[\zeta = \dfrac{\mu}{2\sqrt{mk}}\]

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
%matplotlib nbagg

omega0 = 2. * np.pi * 1.
zeta   = .1
ad = 1.
omegad = 2. * np.pi * 1.2

Part 1: Coding

First, you need to code: the harmonic oscillator using the standarde ODE formulation:

\[\dot X = f(X, t)\]

def f(X, t):
    """
    Driven Harmonic oscillator.
    """
    return

Part 2: Solving

Solve the ODE using odeint and plot it.

Part 3: Steady state amplitude

Determine the steady state amplitude \(a(\omega_d)\) of the oscillator.

Part 4: Amplitude vs. drive pulsation

Plot the evolution of the amplitude of the driven oscillator and compare it with the theory.