# 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()
Comments
Post a Comment