Driven Harmonic Oscillator

Hide code cell source 
# IF USING JUPYTER NOTEBOOK:
#%matplotlib notebook
# IF USING JUPYTER LAB:
%matplotlib widget
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from scipy.integrate import odeint

Driven Harmonic Oscillator#

Author: Ludovic Charleux (ludovic.charleux@univ-smb.fr)

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}} \]
omega0 = 2.0 * np.pi * 1.0
zeta = 0.1
ad = 1.0
omegad = 2.0 * np.pi * 1.2

Part 1: Coding#

First, you need to code: the harmonic oscillator 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.

Contents