numerical methods - Power spectral density in python

Question

I have analytical and numerical solutions for differential equations:

#  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, 100, 1000)
plt.plot(ww, spectrum_density_analytical(ww, x0, p0, omega0, gamma, F, Fm, OMEGA), label='Exact PSD')
plt.title('Spectral density for $gamma$=' + str(gamma) + ', $omega_0$=' + str(omega) + ', $F=$' + str(F) +
     ', $F_m=$' + str(Fm) + ', $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()

And it gives: solutions for x(t) and power spectral densities. As you can see, when omega = 1 analytical solution is correct. But numerical PSD is shifted. Where is my problem?

question from:https://stackoverflow.com/questions/65882894/power-spectral-density-in-python

Answer 1

Waitting for answers