python - Alternatives to Power Spectral Density for Response of an Oscillator - Stack Overflow

admin管理员组

文章数量:1201992

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?

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?

• python
• numpy
• matplotlib
• physics

Share Improve this question Follow edited Jan 21 at 13:52 Sourin Dey asked Jan 21 at 13:42 Sourin DeySourin Dey 5577 bronze badges 3

• Anti-resonance can happen with only one oscillator ? – rehaqds Commented Jan 21 at 16:15
• @rehaqds Yes. The parametric drive (as in the question) will interfere with the external drive, causing antiresonance. – Sourin Dey Commented Jan 21 at 19:00
• You are looking at the response for one forcing frequency. To see a damped frequency you would need to either sweep the frequencies, or force with a broadband input signal (Gaussian white noise is often a reasonable choice) – Jody Klymak Commented Jan 22 at 21:17

Add a comment  | 

1 Answer 1

Sorted by: Reset to default 2

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.

本文标签： pythonAlternatives to Power Spectral Density for Response of an OscillatorStack Overflow