我需要你的帮助，因为我想编码一个正弦激励塔的运动。问题是，当我绘制结果时，有一个看起来不正常的正弦噪声，我不知道它是从哪里来的……我确实期待一个更平滑的曲线，因为它通常是驱动阻尼谐振子的情况。下面是运动方程:
ddx1 + (f1/m1)*dx1 + (k1/m1)*x1 = omega^2*Em*sin(omega*t)
初始条件:x0=0 m, v0=dx0=0 m/s
下面是我的代码:from math import *
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#params
m1=264000000. # kg
f1 = 5000000. # kg/s
k1=225000000. # N/m
#initial displacement of the tower:
x0 = 0. # m
dx0 = 0. # m/s
N=1000000
duration=200
time = np.linspace(0, duration, N)
# Creating the excitation
#sinusoidal excitation
def entry(Em,f,t):
omega = 2*np.pi*f
return -omega**2*Em*np.sin(omega*t)
# Equation: ddx1 + (f1/m1)*dx1 + (k1/m1)*x1 = omega^2*Em*sin(omega*t)
# Solving
def dX(X,t):
#X = [x1, dx1]
A=np.array([[ 0 , 1 ],
[-k1/m1, -f1/m1]])
B=np.array([0, -entry(1,50,t)])
dX=np.dot(A,X)+B
return dX
result = odeint(dX,[x0,dx0],time)
plt.plot(time, result[:, 0])
plt.show()
下面是一些图片:第一张图片当我放大
你能告诉我我的代码有什么问题吗?
提前感谢您的帮助!
[EDIT]我已经尝试了较小频率的代码，它更符合我的期望。我没有想到的是，正如JustLearning指出的那样，固有频率和驱动频率之间的差异非常重要，因此实际上有这些微振荡是正常的。至于这些参数的值，它们确实很重要，因为它们是台北塔的参数。但由于每次都有所有这些量的一个比例，我认为(但我可能是错的)python不会费心去做计算。我真的是新手，所以谢谢你这么快的回答和帮助我。
纯粹根据情节来评估ODE的问题可能是不明智的。为了检查你的代码在运行时是否有意义，你可能应该进行收敛测试:选择一个已知的解析解，并检查|numerical - analytic|的L2范数是否随着你使时间步长变小而减少。
也就是说，仅仅从你的图中观察振荡之上的振荡，你所看到的只不过是非强迫阻尼振荡和强迫项之间的叠加。这个叠加频率之所以如此微小，是因为它大约是50 !比振荡器的自然+阻尼频率大几倍。如果你把50中的entry改变一个更小的值，比如1，你会发现自然+阻尼振荡和强迫都将以大致相似的频率叠加。试试0.1，所有振荡的叠加效果更有可比性。