Die numerisch erhaltene Reaktion eines gedämpften angetriebenen Oszillator ⇐ Python

[Anonymous]

Ich versuche, die Reaktion eines regelmäßig getriebenen gedämpften Oszillators zu zeichnen, dessen Dynamik durch, < /p>
x '' + 2gx ' + f0^2 x = f cos (ft ) < /p>
wobei die Konstanten Folgendes bezeichnen. /p>
f: Häufig fahren < /p>
f: Stärke des Antriebs < /p>
, um dies zu tun, habe ich das obige gelöst Differentialgleichung für x (t). Als Nächst Erreichen Sie es. < /p>
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 + v*dt
v[i+1] = v - 2*G*v*dt - (f0*f0)*u*dt + F*np.cos(f1*t)*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)

< /code>
Da der Oszillator gezwungen ist und einen stationären Zustand erreicht, erwarte ich, dass die Reaktion auf die Spitze um die Antriebsfrequenz ist. Der obige Code bietet jedoch einen Spitzenwert beispielsweise bei F. Was mache ich falsch?