CODEWITHSUNDEEP
pythonscipyode
epx
epx

Reputation: 591

pendulum involving scipy ode

I am using scipy ode to solve the pendulum problem.

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

I plot two figure of theta versus time for different time range, one is (0, 150) and one is (0, 200). What I expect is that the two figure should be same in the time range (0, 150), however, that is not what I observe. Anything wrong with my script? Thanks.enter image description here

enter image description here

Upvotes: 1

Views: 539

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

You probably should use a complete initialization as in

r.set_integrator('dopri5', atol=1, rtol = 1e-8)
r.set_initial_value(Theta0, t=t[0])

which gives you control over the absolute and relative error thresholds and explicitly sets the initial time.

Your system has a Lipschitz constant L of about or larger omega0**2=100. The error propagation is dominated by the factor exp(L*(t_final-t)).

For example, from time 50 to time 150 this factor is exp(100*100) which is about 10^4343. For all practical purposes this means that there is no discernible dependence of the final value on the initial value, i.e., the system is chaotic.

The practical side looks somewhat better than this pessimistic estimate since both curves coincide to about t=10. Which means that, assuming the error tolerance is 1e-8, that exp(L*10)<=1e+8 or L=1.84.... This gives for the interval of size 100 an error magnification factor of exp(L*100)=1e+80 which is still large enough to call the result chaotic.

Experiments with the error thresholds show that the initial divergence point at about t=8.4 is insensitive to the relative error. Further, the divergence produces (in my experiments, not your picture) a difference growth from 1 to 24 in dt=1 which gives about L=3, which is still consistent with the last estimate and much less than the theoretical L=100.

For a better understanding what happens you should also investigate the derivatives plot,

plt.plot(t1,thetadot1,t2,thetadot2)
plt.show()

which strikingly demonstrates the origin of the chaotic behavior.

Upvotes: 2

Related Questions