How can I numerically capture antiresonance dips in a parametrically driven oscillator?

Below is an explanation of what’s likely happening with your numerical approach and how you might modify it to capture the antiresonance dip.

---

Why the PSD Doesn’t Show the Dip

1. Single-Run vs. Frequency Sweep:
   The FFT of one long simulation with fixed driving frequency f only reveals peaks at the frequencies present in that run (here, at f and its sidebands due to the parametric modulation). Antiresonance is observed as a dip in the amplitude when the driving frequency is swept across a range of frequencies. In a fixed‑frequency simulation you never “see” the dip because you aren’t changing f.

2. Destructive Interference Signature:
   Antiresonance comes from cancellation between different response components. When you compute a PSD, you measure the power at a given frequency, but you do not capture the subtle cancellation that appears when comparing responses at nearby driving frequencies.

3. Resolution and Modal Mixing:
   The parametric term “λ sin(2 f t)” not only shifts energy into additional frequency components (sidebands) but also makes the cancellation (dip) very sensitive. The PSD you computed is dominated by the drive; the cancellation does not show up as an “absent frequency” in the spectrum, but rather as a reduction in the amplitude of the primary response when examined as a function of f.

---

How to Capture the Dip Numerically

To observe the antiresonance, you need to map out the steady‑state amplitude versus the driving frequency. Here are two complementary strategies:

1. Frequency Sweep

Implement a loop over different driving frequencies (f) and—for each frequency—solve the ODE long enough for the transient to die out. For each run, extract the steady‑state amplitude (for example, by taking the peak-to-peak value or by filtering the time series) and then plot amplitude versus f.

Example Sketch:

PYTHONCopied
import numpy as np
from matplotlib import pyplot as plt

# Simulation parameters
G = 0.6       # damping coefficient
Om = 3        # natural frequency
L = 0.5       # lambda (parametric strength)
F = 1         # forcing strength

tf = 50
ti = 0
N = 1000000
dt = (tf - ti)/N
t = np.linspace(ti, tf, N)

# Range of driving frequencies to test
drive_freqs = np.linspace(1.5, 2.5, 50)  # adjust range as needed
steady_amp = []  # store measured amplitudes for each drive freq

for f in drive_freqs:
    u = np.empty(N, dtype=float)
    v = np.empty(N, dtype=float)
    u[0] = 1
    v[0] = 0
    
    # Integrate ODE using Euler method (or better: use a higher order integrator)
    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*f*t[i]) * u[i]*dt \
                 + F * np.sin(2*np.pi*f*t[i])*dt

    # Skip transient; use final part of simulation
    slice_index = int(20/dt)
    U = u[slice_index:]
    
    # Simple amplitude estimate (e.g., max-min over steady state)
    amp = (U.max() - U.min()) / 2
    steady_amp.append(amp)

plt.plot(drive_freqs, steady_amp, 'o-')
plt.xlabel('Driving Frequency (Hz)')
plt.ylabel('Steady-State Amplitude')
plt.title('Response Amplitude vs. Driving Frequency')
plt.show()

Notes:

• You might use a more accurate ODE solver (e.g., from scipy.integrate) to improve accuracy.
• The range of drive frequencies and the method to extract amplitude can be refined based on your specific system parameters.

2. Lock-In Amplifier Technique

Another approach is to “lock in” on the drive frequency. Multiply the solution by a reference sine and cosine of frequency f and then low-pass filter the product. This isolates the in-phase and quadrature components so you can compute the amplitude. Again, sweeping f over a range will reveal the dip.

Example (Concept Outline):

PYTHONCopied
# After obtaining u(t) for a given f
# Multiply by reference signals
ref_sin = np.sin(2*np.pi*f*t[slice_index:])
ref_cos = np.cos(2*np.pi*f*t[slice_index:])

# Low-pass filter via averaging (or using a proper filter)
in_phase = np.mean( U * ref_cos )
quad = np.mean( U * ref_sin )
amp = np.sqrt(in_phase**2 + quad**2)

Repeat this for a range of f as done in the frequency-sweep strategy.

---

Summary

• PSD Limitation: A single fixed‑frequency simulation and its FFT show the energy at that drive frequency and its sidebands but do not reveal the cancellation effect that manifests when you vary the driving frequency.
• Recommended Approach: Perform a frequency sweep (or use a lock-in amplification method) to record the steady‑state amplitude versus f. The antiresonance dip should appear as a marked reduction in amplitude at a specific f.

Using one of these approaches, you should be able to numerically capture the antiresonance dip in your system’s response.