CODEWITHSUNDEEP
pythonphysicsodeint
zahbaz
zahbaz

Reputation: 183

Odeint non-consistent over multiple periods, modeling driven pendulum

There seem to be many questions addressing pendulum modelling and odeint. I believe this question is specific enough to stand by itself. It's concerned with passing a time array to odeint.

I'm modelling a driven, damped pendulum. I expect transient behavior to die off after some number of periods and have been plotting my angular velocity vs time to observe this. My problem is that varying the number of periods does not seem to provide consistent results. I don't see where the code or my assumptions fail.

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)

Is the above line problematic? Is it bad form to use non-integral values here? Linspace should still give a constantly spaced array. I've seen other examples do this with success... My idea was to base the time off of the driving period and set some number, 256, points per period. Is this faulty?

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

Below are plots for number of periods (nT) equalling 140,180,and 200. I'd expect to see the continuation of the same behavior, but instead the 180 period result doesn't lose its transient and the 200 result reaches steady state behavior the most quickly! Where is the fault in my logic?

pendulum

Upvotes: 1

Views: 158

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

You have a Lipschitz constant of about L=1, which gives an error magnification factor of exp(L*dT)=exp(dT) for time differences dT. Only considering the normal numerical noise of about 1e-16 it only needs dT=37 to magnify this initial error to a contribution of about 1, as exp(37)*1e-16 = 1.17.

As you see, over the time span from 0 to 1200 or larger, even the slightest variation in the execution of the algorithm will lead to seemingly random variations in the trajectory. You only have the guarantee of at least graphical similarity under these procedural variations for a time span from 0 to about 30.

Upvotes: 1

Related Questions