Why Does My FFT of a Damped Driven Oscillator Peak at an Unexpected Frequency?

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.

PYTHONCopied
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()

There are several things that need to be improved:

• First, the model. As already mentioned the differential equation is written using frequencies but these should be angular frequencies. More importantly, the second order differential equation of a harmonic oscillator is typically parametrized with a natural frequency w0 and a damping ratio with a pre-factor 2*w0. In your equation the w0 is missing from the pre-factor, therefore it is now included in damping G (which thereby becomes resonant-frequency dependent).
• Second, the solution method. The code implements a 'rectangle rule' for integration but that has poor performance and thus needs very large N. A much better approach would be to employ a 4th order Runge Kutta method. Since you already use scipy, try odeint:

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

pi2 = 2 * np.pi

# input:
G=.2  # the damping
f0=2  # natural frequency
fd=5  # driving frequency
F=1  # driving force amplitude

N=5000
T=5

# compute w0 and wd
w0 = pi2 * f0
wd = pi2 * fd

# x''+ 2Gx' + w0^2 x = F cos(wd*t)
# Write this a a system of first order ODE's by defining v = dx/dt:
# x' = v
# v' = F * cos(wd*t) - 2 G*w0*v - w0**2 * x 

def osc(y, t, F, G, wd):
    x, v = y
    dydt = [v, F*np.cos(wd*t) - 2*G*w0*v - w0**2*x]
    return dydt

# the time steps:
t=np.linspace(0,T,N)

# the initial conditions for x and v
y0 = [0, 0.5]

# solve the ODE
sol = odeint(osc, y0, t, args=(F, G, wd))

# and plot the result:
plt.plot(t, sol[:, 0], 'g', label='x(t)')
plt.plot(t, sol[:, 1], 'r', label='v(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()

• Third, to show that the transfer function peaks at the natural frequency you need to run the above simulation for different frequencies and plot the amplitude. Although that is possible, it is much easier to do this from theory. From the top of my head, for a force driven oscillator the differential equation is:

PYTHONCopied
import numpy as np
import matplotlib.pyplot as plt

G=.2  # the damping
w0=2 # natural frequency

w = np.linspace(0, 10, 200)
transfer_function = 1 / (w0**2 + 2*1j*G*w0*w -w**2)
plt.plot(w, abs(transfer_function), 'b', label='|H(w)|')
plt.legend(loc='best')
plt.xlabel('w')
plt.grid()
plt.show()

If the plot above is made on a log scale, it can be seen that at higher frequencies the amplitude decays with 1/w**2 (2 orders of magnitude per frequency decade) unless the damping is made high and becomes the dominant term, in that case the initial slope will be 1/w (one order of magnitude per decade).