Integrazione di ODE rigide con Python

Question

Sto cercando una buona libreria che integri ODE rigide in Python. Il problema è che l'odeint di Scipy mi dà buone soluzioni a volte, ma il minimo cambiamento nelle condizioni iniziali fa sì che cada e si arrenda. Lo stesso problema è risolto abbastanza felicemente dai risolutori rigidi di MATLAB (ode15s e ode23s), ma non posso usarlo (nemmeno da Python, perché nessuno dei binding Python per l'API MATLAB C implementa i callback, e ho bisogno di passare una funzione al risolutore ODE). Sto provando PyGSL, ma è terribilmente complesso. Ogni suggerimento sarà molto apprezzato.

MODIFICA: il problema specifico che sto riscontrando con PyGSL è la scelta della funzione del passo corretto. Ce ne sono diversi, ma non analoghi diretti a ode15s o ode23s (formula in bdf e Rosenbrock modificato se ciò ha senso). Quindi, qual è una buona funzione di passaggio da scegliere per un sistema rigido? Devo risolvere questo sistema per un tempo molto lungo per assicurarmi che raggiunga lo stato stazionario, ei risolutori GSL scelgono un passo temporale minuscolo o troppo grande.

Answer 1

Python può chiamare C. Lo standard di settore è LSODE in ODEPACK. È di dominio pubblico. È possibile scaricare il C version. Questi risolutori sono estremamente difficili, quindi è meglio usare un codice ben collaudato.

Aggiunto: Assicurati di avere un sistema rigido, ad esempio se i tassi (autovalori) differiscono di più di 2 o 3 ordini di grandezza. Inoltre, se il sistema è rigido, ma stai solo cercando una soluzione a stato stazionario, questi risolutori ti danno la possibilità di risolvere alcune equazioni algebricamente. In caso contrario, un buon solutore di Runge-Kutta come DVERK sarà una soluzione buona e molto più semplice.

Aggiunto qui perché non si adatterebbe in un commento: Questo è dal doc intestazione DLSODE:

C  T  :INOUT Value of the independent variable. On return it 
C     will be the current value of t (normally TOUT). 
C 
C  TOUT :IN  Next point where output is desired (.NE. T). 

Inoltre, sì cinetica di Michaelis-Menten è non lineare. L'accelerazione di Aitken funziona comunque con esso. (Se vuoi una breve spiegazione, considera innanzitutto il semplice caso in cui Y è uno scalare: esegui il sistema per ottenere 3 punti Y. Tocca una curva esponenziale attraverso di essi (algebra semplice), quindi imposta Y all'asintoto e Ora basta generalizzare a Y come un vettore: supponi che i 3 punti siano in un piano - va bene se non lo sono.) Inoltre, a meno che tu non abbia una funzione di forzatura (come una flebo IV costante), l'eliminazione di MM decadrà via e il sistema si avvicinerà alla linearità. Spero possa aiutare.

Answer 2

Attualmente sto studiando un po 'di ODE e dei suoi risolutori, quindi la tua domanda è molto interessante per me ...

Da quello che ho sentito e letto, per problemi stiff la strada giusta da percorrere è di scegliere un metodo implicito come funzione di passaggio (correggimi se sbaglio, sto ancora imparando i misteri dei risolutori ODE). Non posso citarti dove leggo questo, perché non ricordo, ma qui c'è un thread da gsl-help in cui è stata posta una domanda simile.

Quindi, in breve, sembra che il metodo bsimp valga la pena di scattare, anche se richiede una funzione jacobiana. Se non riesci a calcolare lo Jacobian, proverò con rk2imp, rk4imp o con uno qualsiasi dei metodi di cambio.

Answer 3

Se è possibile risolvere il problema con Matlab ode15s, si dovrebbe essere in grado di risolverlo con il risolutore vode di scipy. Per simulare ode15s, io uso le seguenti impostazioni:

ode15s = scipy.integrate.ode(f) 
ode15s.set_integrator('vode', method='bdf', order=15, nsteps=3000) 
ode15s.set_initial_value(u0, t0) 

e quindi si può tranquillamente risolvere il problema con ode15s.integrate(t_final). Dovrebbe funzionare abbastanza bene su un problema rigido.

(vedi anche http://www.scipy.org/NumPy_for_Matlab_Users)

Answer 4

PyDSTool avvolge il risolutore Radau, che è un ottimo integratore rigida implicito. Questo ha più setup di odeint, ma molto meno di PyGSL. Il più grande vantaggio è che la tua funzione RHS è specificata come una stringa (tipicamente, sebbene tu possa costruire un sistema usando manipolazioni simboliche) e sia convertita in C, quindi non ci sono chiamate callback lente di Python e l'intera cosa sarà molto veloce.