Numerically obtaining response of a damped driven oscillator gives peak at wrong frequency

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[i] + v[i]*dt
    v[i+1] = v[i] - 2*G*v[i]*dt - (f0*f0)*u[i]*dt + F*np.cos(f1*t[i])*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?

Your frequencies f0 and f1 are applied in the finite-difference model in rad/s. This may or may not have been your intention.

However, your frequencies from the FFT are in cycles/s.

Since you are using the symbol f, rather than omega, I would guess that you want them in cycles/s. In your finite-difference model then you would have to use 2.PI.f in both locations where you put an f before. Specifically in the line

v[i+1] = v[i] - 2*G*v[i]*dt - (2 * np.pi * f0 ) ** 2 * u[i]*dt + F*np.cos( 2 * np.pi * f1*t[i] ) * dt

Then you get peak energy at a frequency of 5 Hz. (Trim the x-axis scale.)

You are very heavily damped, BTW. Also, you aren't strictly plotting PSD.

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[i] + v[i]*dt
    v[i+1] = v[i] - 2*G*v[i]*dt - (2 * np.pi * f0 ) ** 2 * u[i]*dt + F*np.cos(2 * np.pi * f1*t[i] ) * 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.xlim(0,10)
plt.show()