python中的功率谱密度

如何解决python中的功率谱密度

我有微分方程的解析解和数值解：

#  define system of motion equations 
def diff_eq(z,t):
    q,p = z
    return [p,F*np.cos(OMEGA * t) + Fm*np.cos(3*OMEGA * t) - 2*gamma*p + k*omega0**2 * q - beta * q**3]

#  linear for beta = 0
#  analytical
x0 = 1
p0 = 0
omega0 = 1
gamma = 0.0
F = 0.0
Fm = 0.0
OMEGA = 1.4

t = np.linspace(0,100,1000)
plt.plot(t,motion_analytical(t,x0,p0,omega0,gamma,F,Fm,OMEGA)[0],label='Exact',lw=3)

#  numerical
T = 1 / OMEGA
omega = -omega0
beta = 0.0

z0 = [x0,p0]
sol1 = odeintw(diff_eq,z0,t,atol=1e-13,rtol=1e-13,mxstep=1000)
num_q,num_p = sol1[:,0],sol1[:,1]
a = plt.plot(t,num_q,label='ODEint')
#plt.legend(bBox_to_anchor=(1.05,1),loc='upper left')
plt.legend()
#plt.xlim([10,20])
plt.title('Time series for $\gamma$=' + str(gamma) + ',$omega_0$=' + str(omega) + ',$F=$' + str(F) +
     ',$F_m=$' + str(Fm) + ',$\Omega=$' + str(OMEGA))
plt.grid()
plt.show()

#  spectral density
#  analytical
ww = np.linspace(-2,2,1000)
t = np.linspace(0,1000)
plt.plot(ww,spectrum_density_analytical(ww,OMEGA),label='Exact PSD')
plt.title('Spectral density for $\gamma$=' + str(gamma) + ',$\Omega=$' + str(OMEGA))
plt.grid()

#  numerical
signal = num_q
fourier = sp.fft.fft(signal)
n = signal.size
timestep = 0.1
freq = sp.fft.fftfreq(n,d=timestep)
plt.plot(freq,abs(fourier)**2,label='SciPy PSD')
plt.legend()
plt.show()

它给出：solutions for x(t) 和 power spectral densities。如您所见，当 omega = 1 解析解正确时。但是数字PSD发生了变化。我的问题在哪里？