Alternatives to Power Spectral Density for Response of an Oscillator

Question

Background:

Antiresonance is the phenomenon when the amplitude of a system (say an oscillator) shows a significant dip. Here is an article in Wikipedia discussing the same in a coupled system. The article contains an analytical treatment of the system suggesting the existence of antiresonance.

Motivation:

The following equation describes another system where I'm trying to look for the same phenomena.

x'' + 2 Γ x' + f0^2( 1 + λ sin(2ft) ) x = F sin(ft)

where the constants denote the following.

Γ: Damping coefficient

f0: Natural frequency

f: Driving frequently

F: Strength of the drive

λ: Strength of the parametric drive

Implementing a similar treatment, one can show that the above system too exhibits antiresonance. (My derivation can be found here.)

I want to test this numerically by solving the differential equation and obtaining plots similar to the ones depicted in the article (the first image of the article). Therefore, my question is how do I generate the response numerically that will capture the dip?

Attempt:

I've successfully obtained responses of damped driven (periodically and stochastically) oscillators in the past by computing their Power Spectral Density (PSD). Numerically, to achieve that, one has to take a fast Fourier transform (FFT) of the solution x(t). I followed that same route here, solved the equation for x(t), and plotted the PSD by FFTing the solution.

Here is the code that does the above.

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])

However, I only obtain a peak at f=2, the driving frequency, and no dip anywhere. In fact, no choices of parameters showed a dip.

Question:

Am I computing the response wrong? If so, how to do it? Why does PSD fail to capture it?

Answer 1

A log-log plot reveals a sharp drop at 5Hz, and again at ~18Hz. You can switch to logarithmic axes using plt.xscale('log').

For solving ODEs, scipy.integrate has various solvers that might capture temporal dynamics more accurately.

I think it's worth trying your method out on a system for which you know the derivation and values to be correct (e.g. the article you referenced). Once it's confirmed to be working on a simpler reference system, you can try it on your experimental system. I find that breaking it down like that helps with troubleshooting and validation.

#... code from OP

plt.plot(positive_freqs, psd[frequencies > 0])
plt.xscale('log')
plt.yscale('log')
plt.grid(visible=True, which='both', lw=0.5, color='gray', linestyle=':')
plt.gcf().set_size_inches(10, 3)
plt.xlim(0.1, 200)
plt.xlabel('frequency/Hz')
plt.ylabel('PSD')

---

FYI, you might like to see this answer that also uses a finite-difference approach, but in 2D.