Actualizar la condición inicial en ODE solver cada paso de tiempo

Estoy deseando resolver un sistema de EDO donde, durante los primeros 30,000 segundos, quiero que una de mis variables de estado comience desde el mismo valor inicial. Después de esos 30,000 segundos, quiero cambiar el valor inicial de esa variable de estado a algo diferente y simular el sistema por el rest del tiempo. Aquí está mi código:

def ode_rhs(y, t): ydot[0] = -p[7]*y[0]*y[1] + p[8]*y[8] + p[9]*y[8] ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8] ydot[2] = -p[10]*y[2]*y[3] + p[11]*y[9] + p[12]*y[9] ydot[3] = -p[13]*y[3]*y[6] + p[14]*y[10] + p[15]*y[10] - p[10]*y[2]*y[3] + p[11]*y[9] + p[9]*y[8] - p[21]*y[3] ydot[4] = -p[19]*y[4]*y[5] - p[16]*y[4]*y[5] + p[17]*y[11] - p[23]*y[4] + y[7]*p[20] ydot[5] = -p[19]*y[4]*y[5] + p[15]*y[10] - p[16]*y[4]*y[5] + p[17]*y[11] + p[18]*y[11] + p[12]*y[9] - p[22]*y[5] ydot[6] = -p[13]*y[3]*y[6] + p[14]*y[10] - p[22]*y[6] - p[25]*y[6] - p[23]*y[6] ydot[7] = 0 ydot[8] = p[7]*y[0]*y[1] - p[8]*y[8] - p[9]*y[8] ydot[9] = p[10]*y[2]*y[3] - p[11]*y[9] - p[12]*y[9] - p[21]*y[9] ydot[10] = p[13]*y[3]*y[6] - p[14]*y[10] - p[15]*y[10] - p[22]*y[10] - p[21]*y[10] - p[23]*y[10] ydot[11] = p[19]*y[4]*y[5] + p[16]*y[4]*y[5] - p[17]*y[11] - p[18]*y[11] - p[22]*y[11] - p[23]*y[11] ydot[12] = p[22]*y[10] + p[22]*y[11] + p[22]*y[5] + p[22]*y[6] + p[21]*y[10] + p[21]*y[3] + p[21]*y[9] + p[24]*y[13] + p[25]*y[6] + p[23]*y[10] + p[23]*y[11] + p[23]*y[4] + p[23]*y[6] ydot[13] = p[15]*y[10] + p[18]*y[11] - p[24]*y[13] return ydot pysb.bng.generate_equations(model) alias_model_components() p = np.array([k.value for k in model.parameters]) ydot = np.zeros(len(model.odes)) y0 = np.zeros(len(model.odes)) y0[0:7] = p[0:7] t = np.linspace(0.0,1000000.0,100000) r = odeint(ode_rhs,y0,t) 

Entonces, en otras palabras, quiero establecer y0 [1] en el mismo valor (100) cada vez que se llame a odeint durante los primeros 30,000 segundos. Estoy tratando efectivamente de dejar que el sistema se equilibre por un período de tiempo antes de ingresar una señal en el sistema. Pensé en hacer algo como if t < 30000: y0[1] = 100 como la primera línea de mi función ode_rhs() , pero no estoy muy seguro de que funcione.

Parece que desea que y1 (t) permanezca constante (con el valor 100) para la fase de equilibrio. Puede hacer esto asegurándose de que dy1 (t) / dt = 0 durante esta fase. Hay (al menos) dos formas de lograrlo. El primero es modificar el cálculo de ydot[1] en ode_rhs siguiente manera:

 if t < 30000: ydot[1] = 0.0 else: ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8] 

y usa la condición inicial 100 para y[1] .

Tenga en cuenta que esto introduce una discontinuidad en el lado derecho del sistema, pero el solucionador adaptable utilizado por odeint (el código de Fortran LSODA) suele ser lo suficientemente robusto para manejarlo.

Aquí hay un ejemplo autocontenido. He hecho argumentos p y t1 a ode_rhs . t1 es la duración de la fase de equilibrio.

 import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def ode_rhs(y, t, p, t1): ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2] if t < t1: ydot[1] = 0.0 else: ydot[1] = -p[0]*y[0]*y[1] + p[1]*y[2] ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2] return ydot ydot = np.zeros(3) p = np.array([0.01, 0.25, 0.1]) y0 = [20.0, 100.0, 0.0] t = np.linspace(0, 200, 2001) t1 = 20.0 sol = odeint(ode_rhs, y0, t, args=(p, t1)) plt.figure(1) plt.clf() plt.subplot(3, 1, 1) plt.plot(t, sol[:, 0]) plt.axvline(t1, color='r') plt.grid(True) plt.ylabel('y[0]') plt.subplot(3, 1, 2) plt.plot(t, sol[:, 1]) plt.axvline(t1, color='r') plt.grid(True) plt.ylabel('y[1]') plt.ylim(0, 110) plt.subplot(3, 1, 3) plt.plot(t, sol[:, 2]) plt.axvline(t1, color='r') plt.grid(True) plt.ylabel('y[2]') plt.xlabel('t') plt.show() 

Una pequeña variación del método anterior es modificar el sistema agregando un parámetro que sea 0 o 1. Cuando el parámetro es 0, el sistema de eclibración se resuelve y cuando el parámetro es 1, se resuelve el sistema completo. El código para ydot[1] (en mi ejemplo más pequeño) es entonces

 ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2]) 

donde full es el parametro

Para manejar la fase de equilibrio, el sistema se resuelve una vez en 0 <= t full=0 . Luego, el valor final de la solución de equilibrio se usa como condición inicial para la segunda solución, se ejecuta con full=1 . La ventaja de este método es que no está obligando al solucionador a lidiar con la discontinuidad.

Así es como se ve en el código.

 import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def ode_rhs(y, t, p, full): ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2] ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2]) ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2] return ydot ydot = np.zeros(3) p = np.array([0.01, 0.25, 0.1]) y0 = [20.0, 100.0, 0.0] t1 = 20.0 # Equilibration time tf = 200.0 # Final time # Solve the equilibration phase. teq = np.linspace(0, t1, 100) full = 0 soleq = odeint(ode_rhs, y0, teq, args=(p, full)) # Solve the full system, using the final point of the # equilibration phase as the initial condition. y0 = soleq[-1] # Note: the system is autonomous, so we could just as well start # at t0=0. But starting at t1 makes the plots (below) align without # any additional shifting of the time arrays. t = np.linspace(t1, tf, 2000) full = 1 sol = odeint(ode_rhs, y0, t, args=(p, full)) plt.figure(2) plt.clf() plt.subplot(3, 1, 1) plt.plot(teq, soleq[:, 0], t, sol[:, 0]) plt.axvline(t1, color='r') plt.grid(True) plt.ylabel('y[0]') plt.subplot(3, 1, 2) plt.plot(teq, soleq[:, 1], t, sol[:, 1]) plt.axvline(t1, color='r') plt.grid(True) plt.ylabel('y[1]') plt.ylim(0, 110) plt.subplot(3, 1, 3) plt.plot(teq, soleq[:, 2], t, sol[:, 2]) plt.axvline(t1, color='r') plt.grid(True) plt.ylabel('y[2]') plt.xlabel('t') plt.show() 

Y aquí está la ttwig que genera (la ttwig del primer ejemplo es casi exactamente la misma): Parcela de equilibrio y soluciones completas.

Esta respuesta no me funcionó, ya que necesitaba modificar mis condiciones iniciales periódicamente . Por lo tanto, me gustaría proponer una solución alternativa, que es alternar las condiciones para la ecuación diferencial dentro de la función en sí:

Miramos el valor de t y lo ajustamos:

 if int(t) > 20: full = 1 else: full = 0 

Aquí está dentro de una función:

 import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt def ode_rhs(y, t, p, full): if int(t) > 20: full = 1 else: full = 0 ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2] ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2]) ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2] return ydot ydot = np.zeros(3) # intial conditions p = np.array([0.01, 0.25, 0.1]) y0 = [20.0, 100.0, 0.0] t = np.linspace(0, 200, 100) full = 0 # solve equation solution = odeint(ode_rhs, y0, t, args=(p, full)) plt.figure() plt.clf() plt.subplot(3, 1, 1) plt.plot(t, solution[:, 0]) plt.axvline(20, color='r') # vertical line plt.grid(True) plt.ylabel('y[0]') plt.subplot(3, 1, 2) plt.plot(t, solution[:, 1]) plt.axvline(20, color='r') # vertical line plt.grid(True) plt.ylabel('y[1]') plt.ylim(0, 110) plt.subplot(3, 1, 3) plt.plot(t, solution[:, 2]) plt.axvline(20, color='r') # x=20 vertical line plt.grid(True) plt.ylabel('y[2]') plt.xlabel('t') plt.show() 

Eso permite llamar a la función para resolver la ecuación exactamente una vez.

  • No tienes que meterte con las condiciones iniciales del paso anterior.
  • Más fácil de trazar
  • El código es más limpio y fácil de manejar

Y lo que es más importante, ahora puede ajustar los parámetros periódicamente dentro de la ecuación. Por ejemplo, digamos que tiene t = [0: 200] y desea cambiar el valor de lleno cada 20 pasos, puede hacerlo de esta manera:

 if int(t/20) % 2 == 0: full = 1 else: full = 0 

valor alternativo de la variable * completo *