Damped Harmonic Oscillator

dho_weak.py

# Output graph of a weakly damped harmonic oscillator

import numpy as np
import matplotlib.pyplot as plot

def x(A_0, gamma, omega_d, phi_0, t):
  return A_0 * np.exp(-gamma*t) * np.cos(omega_d*t + phi_0)

def exp(A_0, gamma, t):
  return A_0 * np.exp(-gamma*t)

x_0, v_0 = 5.0, 15.0
m, c = 0.5, 0.8
k = 5.0
omega_0 = np.sqrt(k/m)
gamma = 0.5*c/m

if gamma**2 - omega_0**2 < 0: # weak damping
  omega_d = omega_0**2 - gamma**2
  phi_0 = np.arctan(-1/omega_d*(gamma + v_0/x_0))
  A_0 = x_0 / np.cos(phi_0)
  
  t = np.arange(0.0, 5.0, 0.01)
  plot.xlabel('t')
  plot.ylabel('x')
  plot.plot(t, x(A_0, gamma, omega_d, phi_0, t), lw=2)
  plot.plot(t, exp(A_0, gamma, t))
  plot.plot(t, -exp(A_0, gamma, t))
  plot.savefig('dho_weak.png')
else:
  print('Not a condition for weak damping')