数值求解受迫阻尼振子响应时峰值频率与预期不符的问题求助

我正在尝试绘制一个受周期驱动的阻尼振子的响应，它的动力学由以下微分方程支配：

x'' + 2Gx' + f₀²x = F cos(ft)

其中各常数的含义如下：

• G：阻尼系数
• f₀：固有频率
• f：驱动频率
• F：驱动强度

为了实现这个目标，我先求解了上述微分方程得到x(t)，接着提取了x(t)的稳态部分，对其进行傅里叶变换后绘制幅值，以此可视化振子的响应。

下面是我用来实现这个过程的代码：

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)

由于这个振子是受迫驱动的，最终会达到稳态，我预期响应会在驱动频率附近出现峰值。但运行上述代码后，得到的峰值位置完全不在f附近。请问我哪里出错了？

备注：内容来源于stack exchange，提问作者Sourin Dey