Integrar EDOs rígidas con Python
scipy pygsl (4)
Estoy buscando una buena biblioteca que integre EDOs rígidas en Python. El problema es que el principio de scipy a veces me da buenas soluciones, pero el menor cambio en las condiciones iniciales hace que se caiga y se rinda. El mismo problema se resuelve bastante felizmente por los solucionadores rígidos de MATLAB (ode15s y ode23s), pero no puedo usarlo (incluso desde Python, porque ninguno de los enlaces de Python para la API de MATLAB C implementa devoluciones de llamada, y necesito pasar una función al solucionador de ODE). Estoy probando PyGSL, pero es terriblemente 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 hay análogos directos a ode15s u ode23s (fórmula bdf y Rosenbrock modificado si eso tiene sentido). Entonces, ¿cuál es la función de un buen paso para elegir un sistema rígido? Tengo que resolver este sistema durante un tiempo realmente largo para garantizar que alcance el estado estable, y los solucionadores de GSL eligen un paso de tiempo minúsculo o uno que sea demasiado grande.
Actualmente estoy estudiando un poco de EDO y sus solucionadores, por lo que su pregunta es muy interesante para mí ...
Por lo que he escuchado y leído, para problemas rígidos, el camino correcto es elegir un método implícito como una función escalonada (corríjame si me equivoco, todavía estoy aprendiendo los misterios de los solucionadores ODE). No puedo citarte donde leí esto, porque no lo recuerdo, pero aquí hay un thread de gsl-help donde se hizo una pregunta similar.
Entonces, en resumen, parece que el método bsimp vale la pena, aunque requiere una función jacobiana. Si no puedes calcular el jacobiano, lo intentaré con rk2imp , rk4imp , o cualquiera de los métodos de engranajes.
Python puede llamar a C. El estándar de la industria es LSODE en ODEPACK. Es de dominio público. Puedes descargar la versión C Estos solucionadores son extremadamente difíciles, por lo que es mejor usar un código bien probado.
Agregado: asegúrese de que realmente tenga 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 estable, estos solucionadores le dan la opción de resolver algebraicamente algunas de las ecuaciones. De lo contrario, un buen solucionador de Runge-Kutta como DVERK será una solución buena y mucho más simple.
Se agregó aquí porque no cabría en un comentario: Esto es del doc. Del encabezado de 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 no es lineal. Sin embargo, la aceleración de Aitken funciona con ella. (Si desea una breve explicación, primero considere el caso simple de que Y es un escalar. Ejecute el sistema para obtener 3 puntos Y (T). Ajuste una curva exponencial a través de ellos (álgebra simple). Luego, establezca Y en la asíntota y repita. Ahora simplemente generalice a Y como 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 forzado (como un goteo IV constante), la eliminación de MM decaerá De distancia y el sistema se acercará a la linealidad. Espero que ayude.
Si puedes resolver tu problema con los ode15s de Matlab, deberías poder resolverlo con el solucionador de vip scipy. Para simular ode15s , uso la siguiente configuración:
ode15s = scipy.integrate.ode(f)
ode15s.set_integrator(''vode'', method=''bdf'', order=15, nsteps=3000)
ode15s.set_initial_value(u0, t0)
y luego puede resolver su problema con ode15s.integrate(t_final) . Debería funcionar bastante bien en un problema rígido.
(Véase también http://www.scipy.org/NumPy_for_Matlab_Users )
PyDSTool envuelve el solucionador Radau, que es un excelente integrador rígido implícito. 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 (por lo general, aunque puede construir un sistema usando manipulaciones simbólicas) y se convierte en C, por lo que no hay devoluciones de llamada lentas de python y todo será muy rápido.