Utilizzo di dimensioni di passo adattive con scipy.integrate.ode

Question

La (breve) documentazione per scipy.integrate.ode dice che due metodi (dopri5 e dop853) hanno controllo di passi e output denso. Guardando gli esempi e il codice stesso, vedo solo un modo molto semplice per ottenere l'output da un integratore. Vale a dire, sembra che appena un passo in avanti l'integratore da alcuni dt fisso, ottiene il valore della funzione (s) in quel momento, e ripetere.

Il mio problema ha tempi piuttosto variabili, quindi vorrei solo ottenere i valori in qualunque momento passi che deve valutare per raggiungere le tolleranze richieste. Cioè, presto, le cose stanno cambiando lentamente, quindi i passi del tempo di uscita possono essere grandi. Ma dato che le cose diventano interessanti, i passi del tempo di output devono essere minori. In realtà non voglio output denso a intervalli uguali, voglio solo i passi temporali che la funzione adattiva utilizza.

EDIT: uscita Dense

Una nozione relativa (quasi l'opposto) è "uscita densa", per cui le misure adottate sono grandi come stepper preoccupa di prendere, ma i valori della funzione sono interpolate (solitamente con precisione paragonabile all'accuratezza dello stepper) a qualsiasi cosa tu voglia. FORTRAN sottostante scipy.integrate.ode è apparentemente capace di questo, ma ode non ha l'interfaccia. odeint, d'altra parte, si basa su un codice diverso e fa evidentemente un output denso. (È possibile eseguire l'output ogni volta che viene richiamato il proprio lato destro per vedere quando ciò accade e vedere che non ha nulla a che fare con i tempi di output.)

Così ho potuto sfruttare l'adattabilità, purché Potrei decidere i passi temporali di uscita che voglio in anticipo. Purtroppo, per il mio sistema preferito, io non so nemmeno che cosa i tempi approssimativi sono in funzione del tempo, fino a quando ho eseguito l'integrazione. Quindi dovrò combinare l'idea di prendere un passo dell'integratore con questa nozione di output denso.

Answer 1

Da SciPy 0.13.0,

I risultati intermedi della dopri famiglia di risolutori ODE possono ora accessibili da una funzione solout callback.

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


def logistic(t, y, r): 
    return r * y * (1.0 - y) 

r = .01 
t0 = 0 
y0 = 1e-5 
t1 = 5000.0 

backend = 'dopri5' 
# backend = 'dop853' 
solver = ode(logistic).set_integrator(backend) 

sol = [] 
def solout(t, y): 
    sol.append([t, *y]) 
solver.set_solout(solout) 
solver.set_initial_value(y0, t0).set_f_params(r) 
solver.integrate(t1) 

sol = np.array(sol) 

plt.plot(sol[:,0], sol[:,1], 'b.-') 
plt.show() 

Risultato:

Il risultato sembra essere leggermente diverso da Tim D's, anche se entrambi utilizzano lo stesso backend. Sospetto che questo abbia a che fare con la proprietà FSAL di dopri5. Nell'approccio di Tim, penso che il risultato k7 del settimo stadio venga scartato, quindi k1 viene calcolato di nuovo.

Nota: c'è un noto bug with set_solout not working if you set it after setting initial values. È stato risolto a partire da SciPy 0.17.0.

Answer 2

Il metodo integrate accetta un argomento booleano step che indica al metodo di restituire un singolo passaggio interno. Tuttavia, sembra che i risolutori 'dopri5' e 'dop853' non lo supportino.

Il codice seguente mostra come è possibile ottenere i passi interne adottate dal risolutore quando viene utilizzato il risolutore 'vode':

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


def logistic(t, y, r): 
    return r * y * (1.0 - y) 

r = .01 

t0 = 0 
y0 = 1e-5 
t1 = 5000.0 

backend = 'vode' 
#backend = 'dopri5' 
#backend = 'dop853' 
solver = ode(logistic).set_integrator(backend) 
solver.set_initial_value(y0, t0).set_f_params(r) 

sol = [] 
while solver.successful() and solver.t < t1: 
    solver.integrate(t1, step=True) 
    sol.append([solver.t, solver.y]) 

sol = np.array(sol) 

plt.plot(sol[:,0], sol[:,1], 'b.-') 
plt.show() 

Risultato:

Answer 3

Ho cercato in questo per cercare di ottenere lo stesso risultato. Si scopre che è possibile utilizzare un hack per ottenere i risultati passo-passo impostando nsteps = 1 nell'istanza dell'ode. Genererà un UserWarning ad ogni passo (questo può essere catturato e soppresso).

import numpy as np 
from scipy.integrate import ode 
import matplotlib.pyplot as plt 
import warnings 


def logistic(t, y, r): 
    return r * y * (1.0 - y) 

r = .01 
t0 = 0 
y0 = 1e-5 
t1 = 5000.0 

#backend = 'vode' 
backend = 'dopri5' 
#backend = 'dop853' 

solver = ode(logistic).set_integrator(backend, nsteps=1) 
solver.set_initial_value(y0, t0).set_f_params(r) 
# suppress Fortran-printed warning 
solver._integrator.iwork[2] = -1 

sol = [] 
warnings.filterwarnings("ignore", category=UserWarning) 
while solver.t < t1: 
    solver.integrate(t1, step=True) 
    sol.append([solver.t, solver.y]) 
warnings.resetwarnings() 
sol = np.array(sol) 

plt.plot(sol[:,0], sol[:,1], 'b.-') 
plt.show() 

risultato:

Answer 4

Ecco un'altra opzione che dovrebbe funzionare anche con dopri5 e dop853. In sostanza, il solutore chiamerà la funzione logistic() la frequenza necessaria per calcolare i valori intermedi così che è dove abbiamo memorizzare i risultati:

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

sol = [] 
def logistic(t, y, r): 
    sol.append([t, y]) 
    return r * y * (1.0 - y) 

r = .01 

t0 = 0 
y0 = 1e-5 
t1 = 5000.0 
# Maximum number of steps that the integrator is allowed 
# to do along the whole interval [t0, t1]. 
N = 10000 

#backend = 'vode' 
backend = 'dopri5' 
#backend = 'dop853' 
solver = ode(logistic).set_integrator(backend, nsteps=N) 
solver.set_initial_value(y0, t0).set_f_params(r) 

# Single call to solver.integrate() 
solver.integrate(t1) 
sol = np.array(sol) 

plt.plot(sol[:,0], sol[:,1], 'b.-') 
plt.show() 

Answer 5

Cordiali saluti, anche se una risposta è stata già accettata, vorrei sottolineare per il record storico tale output denso e campionamento arbitrario da ovunque lungo la traiettoria calcolata è supportato nativamente in PyDSTool.Ciò include anche una registrazione di tutte le fasi temporali determinate in modo adattivo utilizzate internamente dal risolutore. Questo si interfaccia con entrambi dopri853 and radau5 e genera automaticamente il codice C necessario per interfacciarli con essi piuttosto che fare affidamento su (molto più lenti) callback di funzioni python per la definizione del lato destro. Nessuna di queste funzionalità è fornita in modo nativo o efficiente in nessun altro solutore focalizzato su Python, per quanto ne so.