# Q 3A # d2x/dt2 + gamma * dx/dt + kx = Fcos wt # Exact Sol of Amplitude = F / sqrt((k-w)^2 + gamma^2*w^2) # take gamma = 0.5 , k = 5 , F = 1 import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint k = 5 gam = 0.5 F = 1. def forced(x ,t): return [x[1], F*np.cos(w*t) - gam*x[1] - k*x[0]] x0 = [0., 1.] t = np.linspace(0, 100, 1001) freq = np.linspace(0.5, 5, 1001) A = [] for w in freq: sol = odeint(forced, x0, t) X = sol[:,0] amp = np.max(X[-200:]) A.append(amp) exact = F/(np.sqrt((k-freq**2)**2 + gam**2*freq**2)) plt.plot(freq[::2], A[::2], '--', label="soultion by odeint") plt.plot(freq, exact, 'o', ms=0.5, label ="exact solution") plt.xlabel("frequency") plt.ylabel("Amplitude") plt.legend() plt.show()