python - example - ¿Cómo leer un sistema de ecuaciones diferenciales desde un archivo de texto para resolver el sistema con scipy.odeint?

sympy python

Question

Tengo un sistema grande (> 2000 ecuaciones) de ODE que quiero resolver con python scipy''s odeint.

Tengo tres problemas que quiero resolver (¿quizás tendré que hacer 3 preguntas diferentes?). Para simplificar, los explicaré aquí con un modelo de juguete, pero tenga en cuenta que mi sistema es grande. Supongamos que tengo el siguiente sistema de ODE:

dS/dt = -beta*S
dI/dt = beta*S - gamma*I
dR/dt = gamma*I

con beta = c p I

donde c, p y gamma son parámetros que quiero pasar a odeint.

odeint espera un archivo como este:

def myODEs(y, t, params):
c,p, gamma = params
beta = c*p
S = y[0]
I = y[1]
R = y[2]
dydt = [-beta*S*I,
beta*S*I - gamma*I,
- gamma*I]
return dydt

que luego se puede pasar a odeint así:

myoutput = odeint(myODEs, [1000, 1, 0], np.linspace(0, 100, 50), args = ([c,p,gamma], ))

Genere un archivo de texto en Mathematica, digamos myOdes.txt, donde cada línea del archivo corresponde a la RHS de mi sistema de ODE, por lo que se ve así

#myODEs.txt
-beta*S*I
beta*S*I - gamma*I
- gamma*I

Mi archivo de texto se ve similar a lo que espera odeint, pero todavía no estoy allí. Tengo tres problemas principales:

1. ¿Cómo puedo pasar mi archivo de texto para que odeint entienda que este es el RHS de mi sistema?
2. ¿Cómo puedo definir mis variables de una manera inteligente, es decir, de forma sistemática? Como hay> 2000 de ellos, no puedo definirlos manualmente. Idealmente, los definiría en un archivo separado y leería eso también.
3. ¿Cómo puedo pasar los parámetros (hay muchos de ellos) como un archivo de texto también?

Leí esta pregunta que está cerca de mis problemas 1 y 2 e intenté copiarla (puse directamente valores para los parámetros para no tener que preocuparme por mi punto 3 anterior):

 systemOfEquations = []
with open("myODEs.txt", "r") as fp :
for line in fp :
systemOfEquations.append(line)
def dX_dt(X, t):
vals = dict(S=X[0], I=X[1], R=X[2], t=t)
return [eq for eq in systemOfEquations]
out = odeint(dX_dt, [1000,1,0], np.linspace(0, 1, 5))

pero obtuve el error:

 odepack.error: Result from function call is not a proper array of floats.
ValueError: could not convert string to float: -((12*0.01/1000)*I*S),

Editar: modifiqué mi código a:

 systemOfEquations = []
with open("SIREquationsMathematica2.txt", "r") as fp :
for line in fp :
pattern = regex.compile(r''.+?/s+=/s+(.+?)$'')
expressionString = regex.search(pattern, line)
systemOfEquations.append( sympy.sympify( expressionString) )
def dX_dt(X, t):
vals = dict(S=X[0], I=X[1], R=X[2], t=t)
return [eq for eq in systemOfEquations]
out = odeint(dX_dt, [1000,1,0], np.linspace(0, 100, 50), )

y esto funciona (no entiendo muy bien qué están haciendo las dos primeras líneas del bucle for). Sin embargo, me gustaría hacer el proceso de definición de las variables más automáticas, y todavía no sé cómo usar esta solución y pasar los parámetros en un archivo de texto. En la misma línea, ¿cómo puedo definir los parámetros (que dependerán de las variables) dentro de la función dX_dt?

¡Gracias por adelantado!

Answer 1

Esta no es una respuesta completa, sino más bien algunas observaciones / preguntas, pero son demasiado largas para los comentarios.

dX_dt es llamado muchas veces por odeint con una matriz 1d y tupla t . Usted proporciona t través del parámetro args . y es generado por odeint y varía con cada paso. dX_dt debe simplificarse para que se ejecute rápidamente.

Por lo general, una expresión como [eq for eq in systemOfEquations] se puede simplificar a systemOfEquations . [eq for eq...] no hace nada significativo. Pero puede haber algo sobre systemOfEquations que lo requiera.

Le sugiero que imprima systemOfEquations (para este pequeño caso de 3 líneas), tanto para su beneficio como para el nuestro. Está utilizando sympy para traducir las cadenas del archivo en ecuaciones. Necesitamos ver lo que produce

Tenga en cuenta que myODEs es una función, no un archivo. Se puede importar desde un módulo, que por supuesto es un archivo.

El punto para vals = dict(S=X[0], I=X[1], R=X[2], t=t) es producir un diccionario con el que las expresiones de sympy pueden funcionar. Una función dX_dt más directa (y creo que más rápida) se vería así:

def myODEs(y, t, params):
c,p, gamma = params
beta = c*p
dydt = [-beta*y[0]*y[1],
beta*y[0]*y[1] - gamma*y[1],
- gamma*y[1]]
return dydt

Sospecho que el dX_dt que ejecuta las expresiones generadas por Sympy será mucho más lento que uno "codificado" como este.

Voy a agregar la etiqueta sympy , porque, tal como está escrita, esa es la clave para traducir tu archivo de texto en una función que odeint puede usar.

Me inclinaría a poner la variabilidad de la ecuación en los parámetros t , en lugar de una lista de expresiones sympy.

Eso es reemplazar:

 dydt = [-beta*y[0]*y[1],
beta*y[0]*y[1] - gamma*y[1],
- gamma*y[1]] 

con algo así como

 arg12=np.array([-beta, beta, 0])
arg1 = np.array([0, -gamma, -gamma])
arg0 = np.array([0,0,0])
dydt = arg12*y[0]*y[1] + arg1*y[1] + arg0*y[0]

Una vez que esto sea correcto, las definiciones de argxx se pueden mover fuera de dX_dt y pasar a través de args . Ahora dX_dt es solo un cálculo simple y rápido.

Todo sympy enfoque sympy puede funcionar bien, pero me temo que en la práctica será lento. Pero alguien con más experiencia de sympy puede tener otras ideas.