python - Odeint 在多个时期内不一致，建模驱动摆

Question

关于摆模型和 odeint 似乎有很多问题。我相信这个问题本身就足够具体。它涉及将时间数组传递给 odeint。

我正在模拟一个受驱动的阻尼摆。我预计 transient 行为会在一定数量的周期后消失，并绘制了角速度与时间的关系图来观察这一点。 我的问题是，改变周期数似乎并不能提供一致的结果。我看不出代码或我的假设在哪里失败。

from numpy import *
from scipy.integrate import odeint
import matplotlib.pyplot as plt

#pendulum diff eq
def pendulum(y,t,b,gamma,drivefreq):
    phi,omega=y
    dydt = [omega,-b*omega - sin(phi) + g*cos(drivefreq*t)]
    return dydt

#pendulum parameters: dampening, force amplitude, drivefreq
b=0.05;     g=0.4;      drivefreq=0.7 
args=(b,g,drivefreq) 

#num pts per period, num periods, time array
N=256;  nT=200;
t=linspace(0,nT*2*pi/drivefreq,nT*N)

上面这行有问题吗？在这里使用非整数值是一种不好的形式吗？ Linspace 仍然应该给出一个恒定间隔的数组。我见过其他例子成功地做到了这一点...我的想法是以驾驶时间为基础的时间，并为每个时间段设置一些数字，256 个点。这是有问题吗？

#initial conditions
y0= [0,0] #[phi0,omega0]

#run odeint
out=odeint(pendulum, y0,t,args)
omega  = out[:,1]

#plot ang velocity vs time
fig=plt.figure('ang velocity vs time')
plt.plot(t,omega)

下面是周期数 (nT) 等于 140,180 和 200 的图。我希望​​看到相同行为的延续，但 180 周期结果不会失去其 transient ，而 200 周期结果达到稳定状态行为最快!我的逻辑错在哪里？

Answer 1

您的 Lipschitz 常数约为 L=1，它给出的时间差误差放大系数为 exp(L*dT)=exp(dT) dT。仅考虑大约 1e-16 的正常数值噪声，只需 dT=37 即可将该初始误差放大到大约 1 的贡献，为exp(37)*1e-16 = 1.17。

如您所见，在 0 到 1200 或更大的时间跨度内，即使算法执行中最微小的变化也会导致轨迹看似随机的变化。您只能保证在这些程序变化下，在 0 到大约 30 的时间范围内至少具有图形相似性。