在数字上获得阻尼驱动振荡器的响应在错误的频率下达到峰值。 我正在试图绘制一个定期驱动的振荡器的响应，该振荡器的动力受到，该振荡器的动态受到， x'' + 2gx' + f0^2 x = f cos（ft） 常数表示以下内容。 G：阻尼系数...

Question

G：阻尼系数

F0：固有频率

F：经常驾驶

F：驱动力的力量

为此，我求解了x（t）的上述微分方程。接下来，我从x（t）中提取了稳态部分，进行了傅立叶变换，并绘制了其幅度以可视化振荡器的响应。

试图实现它的代码。

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)

由于振荡器被迫并达到稳定状态，我期望在驾驶频率周围达到峰值的响应。但是，上述代码使一个峰远处靠近f。我在做什么错？

您的频率f0和f1在rad/s的有限差异模型中应用。这可能是您的意图，也可能不是您的意图。 nove，您的FFT频率为循环/s.

由于您使用的是符号f，而不是欧米茄（Omega），我想您希望它们以周期/s为单位。在您的有限差异模型中，您必须在以前放置F的两个位置使用2.pi.f。特别是在线

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

然后您以5 Hz的频率获得峰值能量。 （修剪X轴刻度。）

您的湿度非常严重，顺便说一句。另外，您并不是严格绘制PSD。

Answer 1

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()

有几件事需要改进：

首先，模型。如已经提到的，微分方程是使用频率编写的，但它们应该是角频率。更重要的是，谐波振荡器的二阶微分方程通常用固有频率W0参数化，并且具有前比例2*W0的阻尼比。在您的方程式中，W0缺少了前因子，因此它已包含在阻尼G中（从而成为谐振频率依赖性）。

第二，解决方案方法。该代码实现了集成的“矩形规则”，但性能较差，因此需要非常大的N。一种更好的方法是采用第四阶runge kutta方法。由于您已经使用了Scipy，请尝试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()

第三，为了表明，传递函数在固有频率下达到峰值，您需要以不同的频率运行上述模拟并绘制幅度。尽管这是可能的，但从理论上做到这一点要容易得多。从我的头顶上，用于力量驱动的振荡器：

Answer 2

• 如果以上图是在日志尺度上制作的，则可以看出，在较高的频率下，振幅以1/w ** 2（每十年的频率为2个数量级）衰减，除非使阻尼变高并成为主要项，在这种情况下，初始斜率将为1/w（每十年一个数量级）。