Numerically obtaining response of a damped driven oscillator gives peak at wrong frequency

Question

I am trying to plot the response of a periodically-driven damped oscillator whose dynamics is governed by,

x''+ 2Gx' + f0^2 x = F cos(ft)

where the constants denote the following.

G: Damping coefficient

f0: Natural frequency

f: Driving frequently

F: Strength of the drive

To do so, I solved the above differential equation for x(t). Next, I extracted the steady-state part from x(t), took its Fourier transform, and plotted its magnitude to visualize the response of the oscillator.

Here is the code that attempts to achieve it.

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)

Since the oscillator is forced and reaches a steady state, I expect the response to peak around the driving frequency. However, the above code gives a peak located nowhere near f. What am I doing wrong?

Answer 1

There are several things that need to be improved:

• First, the model. As already mentioned the differential equation is written using frequencies but these should be angular frequencies. More importantly, the second order differential equation of a harmonic oscillator is typically parametrized with a natural frequency w0 and a damping ratio with a pre-factor 2*w0. In your equation the w0 is missing from the pre-factor, therefore it is now included in damping G (which thereby becomes resonant-frequency dependent).

• Second, the solution method. The code implements a 'rectangle rule' for integration but that has poor performance and thus needs very large N. A much better approach would be to employ a 4th order Runge Kutta method. Since you already use scipy, try odeint:

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

pi2 = 2 * np.pi


# input:
G=.2  # the damping
f0=2  # natural frequency
fd=5  # driving frequency
F=1  # driving force amplitude

N=5000
T=5

# compute w0 and wd
w0 = pi2 * f0
wd = pi2 * fd

# x''+ 2Gx' + w0^2 x = F cos(wd*t)
# Write this a a system of first order ODE's by defining v = dx/dt:
# x' = v
# v' = F * cos(wd*t) - 2 G*w0*v - w0**2 * x 

def osc(y, t, F, G, wd):
    x, v = y
    dydt = [v, F*np.cos(wd*t) - 2*G*w0*v - w0**2*x]
    return dydt

# the time steps:
t=np.linspace(0,T,N)

# the initial conditions for x and v
y0 = [0, 0.5]

# solve the ODE
sol = odeint(osc, y0, t, args=(F, G, wd))

# and plot the result:
plt.plot(t, sol[:, 0], 'g', label='x(t)')
plt.plot(t, sol[:, 1], 'r', label='v(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()

• Third, to show that the transfer function peaks at the natural frequency you need to run the above simulation for different frequencies and plot the amplitude. Although that is possible, it is much easier to do this from theory. From the top of my head, for a force driven oscillator the differential equation is:

import numpy as np
import matplotlib.pyplot as plt

G=.2  # the damping
w0=2 # natural frequency

w = np.linspace(0, 10, 200)
transfer_function = 1 / (w0**2 + 2*1j*G*w0*w -w**2)
plt.plot(w, abs(transfer_function), 'b', label='|H(w)|')
plt.legend(loc='best')
plt.xlabel('w')
plt.grid()
plt.show()

If the plot above is made on a log scale, it can be seen that at higher frequencies the amplitude decays with 1/w**2 (2 orders of magnitude per frequency decade) unless the damping is made high and becomes the dominant term, in that case the initial slope will be 1/w (one order of magnitude per decade).

Answer 2

Your frequencies f0 and f1 are applied in the finite-difference model in rad/s. This may or may not have been your intention.

However, your frequencies from the FFT are in cycles/s.

Since you are using the symbol f, rather than omega, I would guess that you want them in cycles/s. In your finite-difference model then you would have to use 2.PI.f in both locations where you put an f before. Specifically in the line

v[i+1] = v[i] - 2*G*v[i]*dt - (2 * np.pi * f0 ) ** 2 * u[i]*dt + F*np.cos( 2 * np.pi * f1*t[i] ) * dt

Then you get peak energy at a frequency of 5 Hz. (Trim the x-axis scale.)

You are very heavily damped, BTW. Also, you aren't strictly plotting PSD.

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 - (2 * np.pi * f0 ) ** 2 * u[i]*dt + F*np.cos(2 * np.pi * 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.xlim(0,10)
plt.show()