[Python] Numerically obtaining response of a damped driven oscillator gives peak at wrong...

I am trying to plot the response of a periodically-driven damped oscillator whose dynamics is governed by,

x''+ 2Gx' + f0^2 x = F cos(ft)

where the constants denote the following.

G: Damping coefficient

f0: Natural frequency

f: Driving frequently

F: Strength of the drive

To do so, I solved the above differential equation for x(t). Next, I extracted the steady-state part from x(t), took its Fourier transform, and plotted its magnitude to visualize the response of the oscillator.

Here is the code that attempts to achieve it.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq

G=1.0
f0=2
f1=5
F=1

N=500000
T=50
dt=T/N

t=np.linspace(0,T,N)
u=np.zeros(N,dtype=float) # Position
v=np.zeros(N,dtype=float) # Velocity

u[0]=0
v[0]=0.5

for i in range(N-1):
u[i+1] = u + v*dt
v[i+1] = v - 2*G*v*dt - (f0*f0)*u*dt + F*np.cos(f1*t)*dt

slice_index=int(20/dt)
U=u[slice_index:]

X_f = fft(U)
frequencies = fftfreq(len(U), dt)
psd = np.abs(X_f)
positive_freqs = frequencies[frequencies > 0]

plt.plot(positive_freqs, psd[frequencies > 0], label="Simulated PSD")
plt.plot(frequencies, psd)



Since the oscillator is forced and reaches a steady state, I expect the response to peak around the driving frequency. However, the above code gives a peak located nowhere near f. What am I doing wrong?

Continue reading...