Tutorial de solución de ecuaciones diferenciales ordinarias con forzamiento temporal determinista¶
Lourdes Martín Aguilar, César Flores López & Marco Herrera Valdez¶
Facultad de Ciencias, Universidad Nacional Autónoma de México¶
Supongamos que la concentración de un metabolito ($x$, mM) en el torrente sanguíneo es influenciada por dos hormonas de manera distinta. Si la concentración del metabolito tiende a un valor $a$ en la ausencia de influencias externas, los cambios en $x$ con respecto al tiempo se pueden determinar resolviendo una ecuación de la forma
\begin{equation}
\partial_t x = \frac{a - x}{\tau} + f(t),
\end{equation}
donde $f$ es una función de forzamiento que representa una influencia externa que comienza en un tiempo $t_0$, a partir de una condición inicial $x(0) = x_0$. Si las influencias de las dos hormonas mencionadas anteriormente se pueden representar con funciones $f_1(t)$ y $f_2(t)$, entonces nos interesa comparar el comportamiento de $x$ con respecto a $f_1$ y $f_2$ en el tiempo.
Las soluciones son obtenidas utilizando las rutinas de integración numérica del módulo scipy de python.
Primero debemos importar los módulos necesarios.
import scipy as sc
import matplotlib.pylab as gr
from scipy.integrate import odeint
%matplotlib inline
Después hacemos un diccionario con los parámetros y arreglos necesarios para la integración numérica.
Las funciones de forzamiento que utilizaremos como ejemplos son: \begin{eqnarray} f_1(t)=\alpha\left(t,t0,\tau{\alpha}, A_1\right)&=& A_1 H(t>t_0)\frac{t-t0}{\tau{\alpha}}\left(1- \frac{t-t0}{\tau{\alpha}} \right) \ f_2(t)=\sin\left(t,t_0,\omega, A_2\right)&=& A_2 H(t>t_0) \sin \left[\left(t-t_0\right)2\pi \right] \end{eqnarray} La función $f_1$ representa un pulso con comienzo súbito en un tiempo $t_0$ y decaimiento exponencial con una constante de tiempo $\tau_{\alpha}$. La función $f_2$ representa una fluctuación periódica de forma sinusoidal que comienza también en $t=t_0$.
Escribimos y graficamos las funciones alfa y seno:
alpha= lambda t,q: q["A1"]*sc.int32(t>q["t0"])*((t-q["t0"])/q["tauAlpha"])*sc.exp(1-(t-q["t0"])/q["tauAlpha"])
sine = lambda t,q: sc.int32(t>q["t0"])* q["A2"]*(1+sc.sin(t*q["freq"]*2*sc.pi))
y las graficamos para verificar que tienen el curso temporal correcto
# Parámetros para usar con las funciones alpha y sine
p=dict()
p["t0"]=5; p["ic"]=10.0;
p["a"]=5.0; p["A1"]=2; p["A2"]=0.3
p["tau"]=5.0; p["tauAlpha"]=0.5; p["height"]=0.3;
p["freq"]=0.3;
p["timeMax"]=30.0; p["timeStep"]=0.1;
p["sampTimes"]= sc.arange(0,p["timeMax"],p["timeStep"])
print(p.keys())
# Graficas de las funciones alpha y sine
fi=gr.figure(figsize=(11,4)); gr.ioff();
ax1=fi.add_subplot(1,2,1); ax2=fi.add_subplot(1,2,2)
ax1.plot(p["sampTimes"],alpha(p["sampTimes"],p),label=r"$f_1(t)$ ((mg/dl)/seg), $t_0$=%g"%(p["t0"]))
ax2.plot(p["sampTimes"],sine(p["sampTimes"],p),label=r"$f_2$ ((mg/dl)/seg), %g Hz, $t_0$=%g"%(p["freq"],p["t0"]))
ax1.set_xlabel("seg"); ax2.set_xlabel("seg")
ax1.legend(); ax2.legend(); gr.ion(); gr.draw()
Definimos una función que calcule el cambio en $x$ para que Python la integre numéricamente:
def linearForcedEq(x,t,par):
dx = (par["a"]-x)/par["tau"] + par['f'](t,par)
return dx #la función integrate.odeint de _scipy_
Ahora integramos la función usando la función integrate.odeint de scipy
p=dict()
p["t0"]=5; p["ic"]=10.0; p["a"]=3.0; p["A1"]=2; p["A2"]=0.3
p["tau"]=5.0; p["tauAlpha"]=0.5; p["freq"]=0.3;
p["timeMax"]=30.0; p["timeStep"]=0.1;
p["sampTimes"]= sc.arange(0,p["timeMax"],p["timeStep"])
# Integración numérica
p["f"]=sine
uSine= odeint(func=linearForcedEq, y0=p["ic"], t=p["sampTimes"],args=(p,)).transpose()
p["f"]=alpha
uAlpha= odeint(func=linearForcedEq, y0=p["ic"], t=p["sampTimes"],args=(p,)).transpose()
Graficamos la solución.
gr.figure(figsize=(11,4))
gr.plot(p["sampTimes"],uAlpha[0],label=r"$f(t)=\alpha(t$,%g,%g$)$"%(p["tauAlpha"],p["t0"]))
gr.plot(p["sampTimes"],uSine[0],label=r"$f(t)=\sin(t$,%g,%g$)$"%(p["freq"],p["t0"]))
gr.ylim(0,p["ic"])
gr.legend(loc="lower left",ncol=2)
