涉及scipy od的摆锤

Question

我用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）应该是相同的，然而，这不是我观察到的。我的剧本有什么问题吗？谢谢。

Answer 1

您可能应该使用完整的初始化，如

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或{}。对于大小为100的区间，这给出了exp(L*100)=1e+80的误差放大因子，它仍然足够大，可以称结果为混沌。在

误差阈值实验表明，初始发散点在t=8.4附近对相对误差不敏感。此外，发散产生了（在我的实验中，而不是你的图片中）从1到24的差异增长，产生了大约{}，这仍然与最后的估计相一致，远远小于理论上的{}。在

为了更好地了解发生了什么，你还应该调查一下衍生工具图

^{pr2}$

这清楚地证明了混沌行为的起源。在