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

Tengo un sistema grande (> 2000 ecuaciones) de EDO que quiero resolver con el método de Python Scipy.

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’s:

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 está esperando 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], )) 

Generé un archivo de texto en Mathematica, por ejemplo, myOdes.txt, donde cada línea del archivo corresponde al 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 es similar a lo que odeint espera, pero todavía no estoy del todo 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 una manera sistemática? Como hay más de 2000, no puedo definirlos manualmente. Lo ideal sería definirlos en un archivo separado y leer eso también.
  3. ¿Cómo puedo pasar los parámetros (también hay muchos) como un archivo de texto?

Leí esta pregunta que está cerca de mis problemas 1 y 2 y traté de 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 me salió 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), 

Edición: modifiqué mi código para:

  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 bien lo que están haciendo las primeras dos líneas del bucle for). Sin embargo, me gustaría hacer el proceso de definir las variables de forma más automática, y todavía no sé cómo usar esta solución y pasar parámetros en un archivo de texto. En la misma línea, ¿cómo puedo definir parámetros (que dependerán de las variables) dentro de la función dX_dt?

¡Gracias por adelantado!

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

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

Generalmente, 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 las systemOfEquations que lo requiera.

Le sugiero que imprima systemOfEquations (para este caso pequeño de 3 líneas), tanto para su beneficio como para el nuestro. Está utilizando sympy para traducir las cadenas del archivo a 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 a vals = dict(S=X[0], I=X[1], R=X[2], t=t) es producir un diccionario con el que puedan trabajar las expresiones sympy . 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 sympy etiqueta sympy , porque, como está escrito, esa es la clave para traducir su archivo de texto en una función que odeint pueda usar.

Me sentiría inclinado a poner la variabilidad de la ecuación en los parámetros t , más bien 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 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, entonces las definiciones de argxx se pueden mover fuera de dX_dt y se pueden pasar a través de args . Ahora dX_dt es un cálculo simple y rápido.

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