scipy odeint Ergebnis hängt von Eingangszeit Array

Question

Ich versuche, die folgende eindimensionale ODE (2. Ordnung) unter Verwendung der scipy odeint Funktion zu lösen:

z'' = 2*F*cos(vt-kz)*(z*cos(vt-k*z) - k*(1+z^2)*sin(vt-kz))/(1+z^2)^2 

Und mein Python-Skript

import numpy as np 
from scipy.integrate import odeint 
import matplotlib.pyplot as plt 

def func(var,t,F,k,v): 
    z,Vz = var 
    tmp = v*t-k*z 
    cos_tmp = np.cos(tmp) 
    denom = 1+z**2 

    return [Vz, 
      2*F*cos_tmp*(z*cos_tmp - k*denom*np.sin(tmp))/(denom*denom)] 

ist

Ich verstehe, dass Odein einen adaptiven Zeitschritt verwendet und es die Lösung zu den Zeitpunkten ausgezahlt, die durch die Eingabe t Array angegeben werden. Ich habe jedoch beobachtet, dass die Lösung divergiert, wenn das Zeitfeld nicht gut genug ist, was für mich keinen Sinn ergibt, wenn Odein den internen Zeitschritt automatisch wählt.

k= 12184. 
F = -0.000125597154176 
v = 0.3 
y0 = [-1.0,0.] 

t1 = np.linspace(0,10,10000) 
t2 = np.linspace(0,10,100) 
t3 = np.linspace(0,10,10) 

result1 = odeint(func,y0,t1,args=(F,k,v)) 
result2 = odeint(func,y0,t2,args=(F,k,v)) 
result3 = odeint(func,y0,t3,args=(F,k,v)) 

fig,ax = plt.subplots(1,2,figsize=(14,5)) 

ax[0].plot(t1,result1[:,0],label='t1') 
ax[0].plot(t2,result2[:,0],'r.',label='t2') 
ax[1].plot(t3,result3[:,0],'r.',label='t3') 

ax[0].legend(frameon=False,loc='upper left',numpoints=1) 
ax[1].legend(frameon=False,loc='upper left',numpoints=1) 

plt.show() 

Ein Bild der Lösung ist hier:

Answer 1

Wenn Sie den Code lief, haben Sie die Warnung bemerken, die gedruckt wird? Es sagt

lsoda-- at current t (=r1), mxstep (=i1) steps 
     taken on this call before reaching tout  
     in above message, i1 =  500 
     in above message, r1 = 0.3622415764602D+00 
[...]/scipy/integrate/odepack.py:218: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. 
    warnings.warn(warning_msg, ODEintWarning) 

es ein wenig kryptisch, aber was es bedeutet, dass odeint den Grenzwert von 500 internen Zeitschritten treffen, bevor eine angeforderte Zeit erreicht.

Sie können dieses Maximum erhöhen, indem Sie das Argument mxstep setzen. Ich habe gerade versucht, Ihren Code mit

result3 = odeint(func,y0,t3,args=(F,k,v), mxstep=2000) 

und es funktionierte.