Sinusoidally forced damped oscillator

We seek a particular solution of the differential equation

$ x'' + 2p x' + w_0^2 x = G_0 \cos(\omega t)$

and we guess

$ x = a \cos(\omega t) + b \sin(\omega t).$

We plug the guess in, and see what $a$ and $b$ need to be.

# sinusoidally_forced_damped_oscillator
from resources306 import *
from sympy import *
a,b,w,w0,p,t  = symbols('a,b,w,w0,p,t',real=True)
G0  = symbols('G0',positive=True)
x = a*cos(w*t) + b*sin(w*t)
#print('The guessed solution:')
#display(x)
de = Eq(diff(x,t,t) + 2*p*diff(x,t) + w0**2*x, G0*cos(w*t) )
dec = de.subs(t,0)
des = de.subs(t,pi/2/w)
print('The two equations to solve for a and b:')
display(dec)
display(des)

sol = solve([dec,des],[a,b])
print('The solution:')
display(sol)
amp = simplify( sqrt(a**2 + b**2).subs(sol) )
print('The amplitude of the solution:')
display( sp.simplify(amp) )

amp = amp.subs({G0:1,w0:1})
from numpy import *
import matplotlib.pyplot as plt
wvals = linspace(0,3,300)
for pval in [2,1,0.5,0.25,0.125,0.0625][::-1]:
             ampp = lambdify( w, amp.subs(p,pval), 'numpy' )
             plt.plot( wvals, ampp(wvals),label='damping p = '+str(pval) )
plt.legend()
plt.title('amplitude of sinusoidally forced damped oscillator with $G_0=1, \omega_0=1$')
plt.xlabel('forcing frequency, $\omega$')
plt.ylabel('amplitude of response');

The two equations to solve for a and b:

$\displaystyle - a w^{2} + a w_{0}^{2} + 2 b p w = G_{0}$

$\displaystyle - 2 a p w - b w^{2} + b w_{0}^{2} = 0$

The solution:

$\displaystyle \left\{ a : \frac{- G_{0} w^{2} + G_{0} w_{0}^{2}}{4 p^{2} w^{2} + w^{4} - 2 w^{2} w_{0}^{2} + w_{0}^{4}}, \ b : \frac{2 G_{0} p w}{4 p^{2} w^{2} + w^{4} - 2 w^{2} w_{0}^{2} + w_{0}^{4}}\right\}$

The amplitude of the solution:

$\displaystyle \frac{G_{0} \sqrt{4 p^{2} w^{2} + \left(w^{2} - w_{0}^{2}\right)^{2}}}{\left|{4 p^{2} w^{2} + w^{4} - 2 w^{2} w_{0}^{2} + w_{0}^{4}}\right|}$