# 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

(1)

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

(2)

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}}

(3)

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

(4)

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)

(5)

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.