Integre ODE rígidos con Python

Question

Estoy buscando una buena biblioteca que integre ODE rígidos en Python. El problema es que scipy's odeint me ofrece buenas soluciones a veces, pero el más mínimo cambio en las condiciones iniciales hace que se caiga y se rinda. El mismo problema lo resuelven felizmente los rígidos solucionadores de MATLAB (ode15s y ode23s), pero no puedo usarlo (incluso desde Python, porque ninguno de los enlaces de Python para la API de C de MATLAB implementa callbacks, y necesito pasar una función al solucionador de ODE). Estoy probando PyGSL, pero es horrendo y complejo. Cualquier sugerencia sería muy apreciada.

EDITAR: El problema específico que tengo con PyGSL es elegir la función de paso correcta. Hay varios de ellos, pero no análogos directos a ode15s u ode23s (fórmula bdf y Rosenbrock modificado si tiene sentido). Entonces, ¿cuál es una buena función de paso para elegir para un sistema rígido? Tengo que resolver este sistema por un tiempo realmente largo para asegurarme de que alcance el estado estacionario, y los solucionadores de GSL eligen un paso de tiempo minúsculo o uno demasiado grande.

Answer 1

Python puede llamar a C. El estándar de la industria es LSODE en ODEPACK. Es de dominio público. Puede descargar el C version. Estos solucionadores son extremadamente difíciles, por lo que es mejor usar un código bien probado.

Agregado: Asegúrese de que realmente tiene un sistema rígido, es decir, si las tasas (valores propios) difieren en más de 2 o 3 órdenes de magnitud. Además, si el sistema es rígido, pero solo está buscando una solución de estado estacionario, estos solucionadores le dan la opción de resolver algunas de las ecuaciones algebraicamente. De lo contrario, un buen solucionador de Runge-Kutta como DVERK será una solución buena, y mucho más simple.

adicional de esto, ya que no encajaría en un comentario: Esto es de la cabecera doc 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). 

Además, sí la cinética de Michaelis-Menten es no lineal. La aceleración Aitken funciona con eso, sin embargo. (Si desea una breve explicación, primero considere el caso simple de que Y sea un escalar. Ejecuta el sistema para obtener 3 puntos Y (T). Ajuste una curva exponencial a través de ellos (álgebra simple). Luego ajuste Y en la asíntota y repita. Ahora generalice que Y es un vector. Suponga que los 3 puntos están en un plano: está bien si no lo están). Además, a menos que tenga una función de forzamiento (como un goteo intravenoso constante), la eliminación de MM disminuirá de distancia y el sistema se acercará a la linealidad. Espero que ayude.

Answer 2

Actualmente estoy estudiando un poco de ODE y sus solucionadores, por lo que su pregunta es muy interesante para mí ...

Por lo que he oído y leído, para los problemas de rigidez en el camino correcto a seguir es elegir un método implícito como función de paso (corrígeme si estoy equivocado, todavía estoy aprendiendo las miserias de los solucionadores de ODE). No puedo citarlo donde lo leí, porque no recuerdo, pero aquí hay un thread de gsl-help donde se hizo una pregunta similar.

Así que, en resumen, parece que vale la pena tomar el método bsimp, aunque requiere una función jacobiana. Si no puede calcular el Jacobiano, lo intentaré con rk2imp, rk4imp, o cualquiera de los métodos del engranaje.

Answer 3

Si resuelve el problema con ode15s de Matlab, debería ser capaz de resolverlo con el vode solucionador de scipy. Para simular ode15s, utilizo los siguientes ajustes:

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

y luego se puede felizmente resolver su problema con ode15s.integrate(t_final). Debería funcionar bastante bien en un problema difícil.

(véase también http://www.scipy.org/NumPy_for_Matlab_Users)

Answer 4

PyDSTool envuelve el solucionador Radau, que es un excelente integrador rígido implícita. Esto tiene más configuración que odeint, pero mucho menos que PyGSL. El mayor beneficio es que su función RHS se especifica como una cadena (normalmente, aunque puede construir un sistema utilizando manipulaciones simbólicas) y se convierte en C, por lo que no hay devoluciones de llamada de pitón lento y todo será muy rápido.