Odeint nicht konsistent über mehrere Zeiträume, Modellierung angetriebenen Pendel

Question

Es scheint viele Fragen zu Pendelmodellierung und Odeint zu sein. Ich glaube, dass diese Frage spezifisch genug ist, um für sich selbst zu stehen. Es geht darum, ein Zeitfeld an Odeint zu übergeben.

Ich modelliere ein angetriebenes, gedämpftes Pendel. Ich erwarte, dass vorübergehendes Verhalten nach einigen Perioden absterben wird und habe meine Winkelgeschwindigkeit gegen die Zeit aufgetragen, um dies zu beobachten. Mein Problem ist, dass das Variieren der Anzahl der Perioden keine konsistenten Ergebnisse zu liefern scheint. Ich sehe nicht, wo der Code oder meine Annahmen versagen.

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) 

Ist die obige Zeile problematisch? Ist es schlecht, hier nicht-integrale Werte zu verwenden? Linspace sollte immer noch ein Array mit konstantem Abstand ergeben. Ich habe andere Beispiele gesehen, die das mit Erfolg machen ... Meine Idee war es, die Zeit für die Lenkzeit zugrunde zu legen und eine Anzahl von 256 Punkten pro Periode festzulegen. Ist das fehlerhaft?

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

Im Folgenden sind grafische Darstellungen für die Anzahl der Perioden (nT) gleich 140.180 und 200. Ich würde erwarten, dass die Fortsetzung des gleichen Verhaltens zu sehen, sondern das 180 Periodener verliert nicht seine vorübergehende und die 200 Ergebnis erreicht stationäres Verhalten am schnellsten! Wo ist der Fehler in meiner Logik?

Answer 1

Sie haben eine Lipschitz-Konstante von etwa L=1, die dT einen Fehler Vergrößerungsfaktor von exp(L*dT)=exp(dT) für Zeitdifferenzen gibt. Nur unter Berücksichtigung des normalen numerischen Rauschens von etwa 1e-16 muss nur dT=37 verwendet werden, um diesen anfänglichen Fehler auf einen Beitrag von etwa 1, wie exp(37)*1e-16 = 1.17, zu vergrößern.

Wie Sie sehen, wird in der Zeitspanne von 0 bis 1200 oder größer, selbst die geringste Variation in der Ausführung des Algorithmus zu scheinbar zufälligen Schwankungen der Flugbahn führen. Sie haben nur die Garantie von mindestens graphischer Ähnlichkeit unter diesen prozeduralen Variationen für eine Zeitspanne von 0 bis ungefähr 30.