Альтернативы спектральной плотности мощности для отклика генератора ⇐ Python

[Anonymous]

Общая информация:
Антирезонанс — это явление, когда амплитуда системы (скажем, осциллятора) значительно падает. Вот статья в Википедии, в которой обсуждается то же самое в связанной системе. В статье содержится аналитическая трактовка системы, предполагающая существование антирезонанса.
Мотивация:
Следующее уравнение описывает другая система, в которой я пытаюсь искать те же явления.
x'' + 2 Γ x' + f0^2( 1 + λ sin(2ft)) x = F sin (футов)
где константы обозначают следующее.
Γ: Коэффициент демпфирования
f0: Собственная частота
f: Частое вождение
F: Сила привода
λ: Сила параметрического привода
Реализация аналогичной обработки , можно показать, что указанная система также обладает антирезонансом. (Мой вывод можно найти здесь.)
Я хочу проверить это численно, решив дифференциальное уравнение и получив графики, подобные тем, которые изображены в статье (первое изображение статьи). . Поэтому мой вопрос: как мне сгенерировать численный ответ, который будет отражать падение?
Попытка:
Я В прошлом я успешно получал отклики затухающих (периодических и стохастических) осцилляторов, вычисляя их спектральную плотность мощности (PSD). Численно, чтобы добиться этого, необходимо выполнить быстрое преобразование Фурье (БПФ) решения x(t). Я пошел по тому же пути, решил уравнение для x(t) и построил PSD, используя БПФ.
Вот код, который делает вышеописанное.

Код: Выделить всё

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

G=0.6       #Gamma
Om=3        #Omega_zero
L=0.5       #Lambda
F=1         #Forcing
Og=2        #Drive_Freq

N=1000000
tf=50
ti=0
dt=(tf-ti)/N

t=np.linspace(ti, tf, N)
u=np.empty(N, dtype=float)
v=np.empty(N, dtype=float)
u[0]=1
v[0]=0

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*Om)**2 *u[i]*dt - (2*np.pi*Om)**2 *L*np.sin(2*np.pi*2*Og*t[i])*u[i]*dt + F*np.sin(2*np.pi*Og*t[i])*dt

slice_index=int(20/dt) #Take the steady state part
U=u[slice_index:]
X_f = fft(U)
frequencies = fftfreq(len(U), dt)
psd = np.abs(X_f)**2

positive_freqs = frequencies[frequencies > 0]
plt.plot(positive_freqs, psd[frequencies > 0])

Однако я получаю пик только при f=2, частоте возбуждения, и никакого провала нигде. Фактически, ни один из вариантов выбора параметров не показал провала.
Вопрос:
Я неправильно рассчитал ответ? Если да, то как это сделать? Почему PSD не может его захватить?

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