作者热门文章
- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我正在使用 scipy ode 来解决钟摆问题。
from scipy import *
import matplotlib.pyplot as plt
from scipy.integrate import ode
def pendulumdot(t, y, gamma, omega0, Fd):
theta, v = y[0], y[1]
return array([v, -2*gamma*v-omega0**2*sin(theta)+Fd*cos(2*pi*t)])
def pendulum_sample(theta0, thetadot0, gamma, omega0, Fd, t):
Theta0 = array([theta0, thetadot0])
r = ode(pendulumdot)
r.set_integrator('dopri5')
r.set_initial_value(Theta0)
r.set_f_params( gamma, omega0, Fd)
#r.set_jac_params()
theta = zeros(len(t))
thetadot = zeros(len(t))
theta[0] = theta0
thetadot[0] = thetadot0
for n in range(1, len(t)):
r.integrate(t[n])
assert r.successful()
theta[n] = (r.y)[0]
thetadot[n] = (r.y)[1]
return theta, thetadot
def pendulum_demo():
gamma = 0.1
omega0 = 10.0
theta0 = 0.0
thetadot0 = 0.0
Fd = 50.0
t1 = linspace(0, 200, 10000)
theta1, thetadot1 = pendulum_sample(theta0, thetadot0, gamma, omega0, Fd, t1)
plt.plot(t1, theta1)
t2 = linspace(0, 150, 10000)
theta2, thetadot2 = pendulum_sample(theta0, thetadot0, gamma, omega0, Fd, t2)
plt.plot(t2, theta2)
plt.show()
pendulum_demo()
我绘制了不同时间范围内的两个 θ 与时间的关系图,一个是 (0, 150),一个是 (0, 200)。我期望两个数字在时间范围 (0, 150) 内应该相同,但是,这不是我观察到的。我的脚本有什么问题吗?谢谢。
最佳答案
您可能应该使用完整的初始化,如下所示
r.set_integrator('dopri5', atol=1, rtol = 1e-8)
r.set_initial_value(Theta0, t=t[0])
它使您可以控制绝对和相对误差阈值并显式设置初始时间。
<小时/>您的系统有一个 Lipschitz 常数 L
大约或更大omega0**2=100
。误差传播主要由因子 exp(L*(t_final-t))
决定。 。
例如,从时间 50 到时间 150,此因子为 exp(100*100)
大约是 10^4343
。出于所有实际目的,这意味着最终值对初始值没有明显的依赖性,即系统是混沌的。
实际情况看起来比这种悲观的估计要好一些,因为两条曲线大约重合于 t=10
。这意味着,假设误差容限为 1e-8
,即exp(L*10)<=1e+8
或L=1.84...
。对于大小为 100 的间隔,误差放大系数为 exp(L*100)=1e+80
它仍然足够大,可以说结果是困惑的。
误差阈值实验表明初始分歧点约为 t=8.4
对相对误差不敏感。此外,分歧在 dt=1
中产生(在我的实验中,不是你的图片)从 1 到 24 的差异增长。这给出了大约L=3
,仍然与最后的估计一致,并且远小于理论L=100
.
为了更好地理解发生的情况,您还应该研究导数图,
plt.plot(t1,thetadot1,t2,thetadot2)
plt.show()
这引人注目地展示了困惑行为的起源。
关于python - 涉及 scipy 颂歌的钟摆,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/29204533/
我是一名优秀的程序员,十分优秀!