# How to plot forcing ocillation with damping correctly? [closed]

The forced oscillation with damping is described as

$$\frac{d^2x}{dt^2}+\gamma\frac{dx}{dt}+\omega_0^2x=\frac{F_0}{m}\cos{\omega t}$$

Its solution is $x=A\cos{(\omega t-\delta)}$, with

$$A(\omega)=\frac{F_0/m}{\sqrt{(\omega_0^2-\omega^2)^2+(\gamma\omega)^2}}$$ $$\tan{\delta(\omega)}=\frac{\gamma\omega}{\omega_0^2-\omega^2}$$

However, I'm unable to reproduce its plots. Set $\gamma=0.5, \omega_0=2$ and time and other constant to $1$, here are how I get: • $A(\omega)$: plot 1/sqrt((4-x^2)^2+(0.5x)^2  Where are my wrongs?

• Comments are not for extended discussion; this conversation has been moved to chat. – ACuriousMind Sep 9 '17 at 14:51

Your plot of the steady state amplitude as a function of frequency looks correct - although usually you would use a strictly positive range of frequencies, rather than a range from - 5 to 5 (which is why it's a bit unusual).

The plot of the phase shift looks strange because your axis wraps at $\pi/2$. At resonance, you have a $\pi/2$ phase difference between driving force and response amplitude - if you didn't wrap the plot, it would show this. As you figured out yourself after a hint from me, you could do this with plot atan2(0.5x,4-x^2) - the atan2 function takes account of the sign of both the x and y coordinate of the vector, and so returns a result over the domain $<-\pi, \pi]$.

The third plot, of the time evolution of the amplitude, is completely wrong. For a given driving frequency, the phase shift will be fixed; you would be left with a differential equation for which there are some complex-looking solutions, but I like to do numerical integration to get a feeling for things. Here's an example of the output of such an integration, and the Python code that generated it:

And the code (using very simple numerical integration... nothing fancy):

# example of numerical integration of 1D driven oscilla
from math import sin, sqrt
import matplotlib.pyplot as plt

# constants to give resonant frequency of about 1 /s
k = 2.0
m = 0.5
g = 0.05 # light damping
F = 1.0
w0 = sqrt(k/m)

w = 0.99 * w0
tmax = 150
dt = 0.01
t = 0
x = 0
v = 0
a = 0

# equation of motion
# m d2x/dt2 - 2 g dx/dt + kx = F sin wt

print('w0 = %f\n'%w0)
def force(x, v, t):
return -g*v - k*x + F*sin(w*t)

X = [0.]
T = [0.]
FF = [1.]

while t<tmax:
f = force(x,v, t)
a = f / m
x = x + v * dt + 0.5*a *dt*dt
v = v + a * dt
FF.append(F*sin(w*t))
t = t + dt
X.append(x)
T.append(t)

plt.figure()
plt.plot(T,X)
plt.plot(T, FF)
plt.title('amplitude response of lightly damped driven oscillator')
plt.legend(('displacement','force'))
plt.xlabel('time')
plt.show()


Play around with the driving frequency and you will get some unexpected results... try w=0.9*w0 for example and you will see some "beating" (which only happens with light damping, and sufficiently far from resonance).