Численное получение отклика затухающего ведомого генератора ⇐ Python

[Anonymous]

Я пытаюсь построить график реакции затухающего генератора с периодическим приводом, динамика которого определяется соотношением:
x''+ 2Gx' + f0^2 x = F cos(ft )
где константы обозначают следующее.
G: Коэффициент затухания
f0: Собственная частота< /p>
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. Что я делаю не так?

Подробнее здесь: https://stackoverflow.com/questions/793 ... oscillator