Retardos y dependencia de la historia

Del estado presente a la historia reciente

Una variable lenta puede conservar información sobre la trayectoria reciente, como mostró el capítulo 47. Esa memoria es aparente, porque una vez conocido el estado completo no hace falta consultar valores anteriores para calcular la derivada actual. En este capítulo la ley de evolución depende explícitamente de estados pasados. La forma más sencilla es

El término $x(t-\tau)$ cambia qué información forma el estado y la manera de especificar las condiciones iniciales. Incluso una ecuación escalar adquiere un espectro infinito de modos y puede perder estabilidad mediante oscilaciones, algo imposible en una ecuación diferencial ordinaria escalar autónoma.

En una ecuación ordinaria $\dot x=f(x)$, conocer $x(t_0)$ basta, con las hipótesis habituales, para determinar el futuro. Con un retardo fijo, conocer $x(0)$ no permite calcular $\dot x(0)$, porque también hace falta $x(-\tau)$. Tampoco bastan esos dos valores aislados. Durante el intervalo $0<t<\tau$, la derivada usa $x(t-\tau)$, que recorre todos los tiempos entre $-\tau$ y $0$. La condición inicial adecuada es una función de historia,

y el espacio de estados natural, en el caso escalar con historias continuas, es el espacio de funciones $C([-\tau,0],\mathbb R)$, que tiene dimensión infinita (Hale y Verduyn Lunel, 1993). La variable observada sigue siendo una sola magnitud, aunque el estado matemático en el instante $t$ es el segmento de trayectoria

La evolución transforma un segmento de historia en el siguiente y el presente $x(t)=x_t(0)$ es solo su extremo derecho. El capítulo 2 ya advirtió que una variable medida no tiene por qué ser un estado suficiente. Aquí el motivo es explícito, porque la propia ecuación consulta el pasado.

Dos historias con el mismo presente

Consideremos $\dot x(t)=-x(t)-2x(t-1)$ con dos historias en $[-1,0]$, $\phi_1(s)=1$ y $\phi_2(s)=1-0{,}3s$. Ambas valen 1 en $s=0$, aunque $\phi_1(-1)=1$ y $\phi_2(-1)=1{,}3$. Las derivadas iniciales son, por tanto, $\dot x_1(0)=-1-2=-3$ y $\dot x_2(0)=-1-2{,}6=-3{,}6$. El mismo valor presente produce velocidades distintas porque las historias difieren. Las trayectorias siguen separadas después (figura 48.1a). En $t=1$ valen $-0{,}896$ y $-1{,}055$. En $t=5$ valen $-0{,}198$ y $-0{,}121$.

El ejemplo ilustra una consecuencia de la estructura del problema, cuya dimensión infinita proviene de su formulación. Tampoco se sigue de él que toda dependencia de la historia requiera una ecuación con retardo. En el capítulo 47, una variable lenta no observada producía una dependencia aparente del pasado dentro de una ecuación ordinaria. Elegir entre las dos representaciones es una decisión de modelamiento que, en una aplicación real, requiere evidencia.

Retroalimentación inmediata y retardada

Estudiaremos el sistema lineal

Si $\tau=0$ se obtiene $\dot x=-(a+b)x$, cuya solución $x_0e^{-(a+b)t}$ converge exponencialmente a cero. Con $\tau>0$ el equilibrio sigue siendo $x^*=0$, porque un estado constante cumple $x(t)=x(t-\tau)$. El retardo no mueve el equilibrio y puede, en cambio, cambiar su estabilidad. En los capítulos anteriores vimos bifurcaciones en las que los equilibrios aparecían, desaparecían o intercambiaban estabilidad. Aquí el parámetro $\tau$ altera el espectro de las perturbaciones alrededor del mismo equilibrio. La intuición es que una corrección basada en información antigua puede actuar cuando el sistema ya cambió y producir una sobrecorrección, aunque decidir si eso ocurre exige resolver el problema espectral.

Una ecuación característica trascendente

Al buscar soluciones $x(t)=e^{\lambda t}$ se obtiene $x(t-\tau)=e^{\lambda t}e^{-\lambda\tau}$. Sustituyendo y dividiendo por $e^{\lambda t}$ resulta la ecuación característica

A diferencia de lo que ocurre en los sistemas lineales de dimensión finita, esta ecuación no es un polinomio. La exponencial contiene a la incógnita y genera infinitas raíces complejas (Smith, 2011). La estabilidad asintótica lineal exige que todas tengan parte real negativa. Una sola raíz con parte real positiva produce un modo creciente. Un par complejo que alcanza el eje imaginario abre, a su vez, la posibilidad de una pérdida oscilatoria de estabilidad.

En este problema las raíces pueden escribirse con la función $W$ de Lambert, definida como la inversa de $w\mapsto we^w$ (Corless et al., 1996). La ecuación equivale a $(\lambda+a)e^{\lambda\tau}=-b$. Con $y=(\lambda+a)\tau$ queda $ye^y=-b\tau e^{a\tau}$, de modo que

Cada rama $W_k$ produce una raíz. El espectro es infinito y numerable, a diferencia del de una matriz de dimensión fija.

El retardo crítico y la transversalidad

Una raíz cruza el eje imaginario cuando $\lambda=i\omega$ con $\omega>0$. Al separar las partes real e imaginaria de la ecuación característica se obtienen

La ecuación de la parte real exige $\cos(\omega\tau)=-a/b$, lo que requiere $b>a$. Al elevar ambas al cuadrado y sumar resulta $a^2+\omega^2=b^2$, por lo que la frecuencia de cruce es $\omega_c=\sqrt{b^2-a^2}$. El primer cruce ocurre cuando $\omega_c\tau_c=\arccos(-a/b)$, es decir,

Hay cruces posteriores, que se obtienen sumando múltiplos de $2\pi$ a la fase, aunque el primero es el que decide la pérdida inicial de estabilidad. Si $a\ge b>0$, el equilibrio es estable para todo retardo no negativo (Erneux, 2009). Introducir un retardo no obliga, entonces, a desestabilizar el sistema.

Encontrar raíces imaginarias tampoco basta para afirmar que la estabilidad cambia. Hay que saber hacia dónde se desplazan al aumentar $\tau$. Derivando implícitamente la ecuación característica respecto de $\tau$ y usando en el cruce la relación $be^{-\lambda\tau}=-(\lambda+a)$, se obtiene

El par cruza transversalmente del semiplano izquierdo al derecho cuando $\tau$ aumenta, así que el primer cruce es una pérdida de estabilidad. El cálculo es local en el parámetro y en el estado. No dice qué ocurre lejos del equilibrio ni garantiza que una versión no lineal desarrolle un ciclo estable.

Un ejemplo con a = 1 y b = 2

Con $\dot x(t)=-x(t)-2x(t-\tau)$ se cumple $b>a$. La frecuencia de cruce es $\omega_c=\sqrt3\approx1{,}7321$ y, como $\arccos(-1/2)=2\pi/3$, el retardo crítico vale

El período del modo crítico es $2\pi/\omega_c\approx3{,}628$. La derivada del cruce vale $0{,}324$.

El espectro revela una historia más rica que la simple pérdida de estabilidad (figura 48.2a). En este ejemplo, la rama principal $W_0$ da la raíz de mayor parte real, como se comprueba evaluando varias ramas. Para retardos pequeños esa raíz es real y más negativa que $-3$, por lo que el retardo acelera al comienzo el decaimiento. Con $\tau=0{,}1$ la raíz dominante vale $-3{,}977$. Cuando $2\tau e^{\tau}=1/e$, es decir, en $\tau\approx0{,}157$, las dos raíces reales dominantes se juntan y forman un par complejo. Desde ese valor la respuesta estable es oscilatoria y su decaimiento se vuelve más lento a medida que $\tau$ aumenta. Con $\tau=0{,}70\,\tau_c$ la raíz dominante es $-0{,}210\pm2{,}254\,i$ y las oscilaciones se amortiguan. En $\tau_c$ el par llega al eje imaginario. Con $\tau=1{,}05\,\tau_c$ la parte real vale $0{,}018$ y la amplitud crece (figura 48.1b).

La solución supercrítica del modelo lineal crece sin cota. Esa divergencia indica que la aproximación lineal ya no contiene un mecanismo que limite la amplitud, sin predecir una divergencia física. En un modelo no lineal, la saturación u otros términos pueden contener el crecimiento y producir una oscilación acotada.

Qué muestra el espectro que no muestra una trayectoria

El espectro responde a una pregunta distinta de la que responde una serie temporal. Cerca del umbral, la tasa de crecimiento puede ser muy pequeña, como el $0{,}018$ del ejemplo. Una simulación corta puede no revelar ese crecimiento. La ausencia de crecimiento en una ventana finita no demuestra estabilidad. El signo de la raíz dominante, calculado con la ecuación característica, decide la estabilidad lineal con independencia de cuánto se prolongue la simulación.

La aproximación de retardo pequeño

Cuando $\tau$ es pequeño respecto de la escala en que cambia $x$, parece natural desarrollar $x(t-\tau)\approx x(t)-\tau\dot x(t)$. En el modelo lineal esto da $(1-b\tau)\dot x\approx-(a+b)x$ y, con $a=1$ y $b=2$, la tasa aproximada

La aproximación reproduce la aceleración inicial del decaimiento, aunque con error apreciable. Con $\tau=0{,}1$ predice $-3{,}75$, frente a $-3{,}977$. Con $\tau=0{,}3$ predice $-7{,}5$, mientras que la raíz dominante exacta es el par complejo $-2{,}548\pm4{,}007\,i$, con un decaimiento tres veces más lento y oscilaciones que la aproximación no puede describir. La fórmula se vuelve singular en $\tau=1/b=0{,}5$, pero el sistema exacto pierde estabilidad recién en $\tau_c\approx1{,}209$.

La discrepancia no es un error de álgebra. El desarrollo supone que $\tau$ es lo bastante pequeño para controlar los términos omitidos, hipótesis que falla mucho antes de la singularidad artificial. La ecuación con retardo tiene, además, infinitas ramas espectrales que la aproximación de primer orden reduce a una sola ecuación ordinaria. En este ejemplo, la aproximación describe razonablemente la corrección de la tasa para retardos menores que 0,1 y falla por completo al localizar la pérdida de estabilidad.

La logística retardada

Un retardo puede participar en el nacimiento de una oscilación acotada cuando hay no linealidad. Un ejemplo histórico es la ecuación logística de Hutchinson (1948),

Con $y=N/K$, $s=rt$ y $\rho=r\tau$ queda

con equilibrio positivo $y^*=1$. Al escribir $y=1+u$ y conservar el primer orden se obtiene $du/ds=-u(s-\rho)$, cuya ecuación característica es $\lambda+e^{-\lambda\rho}=0$. En un cruce imaginario $\lambda=i\omega$ se necesita $\cos(\omega\rho)=0$ y $\omega=\sin(\omega\rho)$, por lo que el primer cruce ocurre con $\omega=1$ y $\rho_c=\pi/2$. En unidades originales, $r\tau_c=\pi/2$. El período del modo crítico, $2\pi$ en tiempo adimensional, es cuatro veces el retardo crítico.

La combinación decisiva es el producto adimensional $r\tau$. Un mismo retardo físico puede ser pequeño para una dinámica lenta y grande para una rápida, de modo que la duración de un retardo, expresada en segundos o en días, no mide por sí misma su importancia dinámica.

La linealización indica cuándo pierde estabilidad el equilibrio. El término no lineal organiza la dinámica de amplitud finita (figura 48.3). Con $\rho=1{,}4$ la población vuelve al equilibrio. Con $\rho=1{,}7$ oscila entre $0{,}334$ y $1{,}943$ con período $6{,}96$, cerca de cuatro veces el retardo. Con $\rho=2$ oscila entre $0{,}069$ y $2{,}90$, con período $8{,}80$. La amplitud crece rápidamente al alejarse del umbral y los mínimos se acercan a cero. Mackey y Glass (1977) mostraron que una ecuación escalar con retroalimentación retardada puede producir regímenes todavía más complejos, incluido el caos. Una ecuación con retardo escalar no es un sistema dinámico de dimensión uno. Las restricciones que valen para los flujos de dimensión finita baja no se le aplican.

Retardos explícitos, variables ocultas y cadenas de compartimentos

Supongamos que una señal tarda en propagarse a través de varios compartimentos intermedios. El proceso puede modelarse con un retardo fijo o con variables adicionales que describan las etapas de tránsito. Las dos construcciones son distintas. Un retardo discreto concentra toda la influencia en el valor exacto $x(t-\tau)$. Una cadena de compartimentos distribuye los tiempos de tránsito. Solo en ciertos límites una sucesión de etapas aproxima un retardo casi fijo (Smith, 2011).

Un sistema puede contener además varios retardos, depender de un promedio ponderado del pasado,

o tener un retardo que depende del propio estado. Estas extensiones comparten la idea de que la evolución usa información distribuida en el tiempo, aunque no comparten la ecuación característica ni las condiciones de estabilidad. Las conclusiones obtenidas para un retardo discreto constante no se trasladan sin demostración.

Por eso, observar una demora entre dos señales no identifica la representación correcta. La demora aparente puede provenir del transporte, de un procesamiento intermedio, de una variable lenta no observada o del filtrado del instrumento. Una correlación desplazada en el tiempo no establece por sí sola una ecuación con retardo ni una relación causal. Para interpretar $\tau$ como un parámetro mecanístico hace falta evidencia independiente sobre el proceso que genera la demora.

El método de pasos

Con un retardo constante existe una idea numérica natural. Durante el primer intervalo, $0\le t\le\tau$, el valor retardado pertenece a la historia conocida $\phi$, por lo que la ecuación con retardo se convierte en una ecuación ordinaria no autónoma. Una vez resuelto ese tramo, en $\tau\le t\le2\tau$ el valor retardado pertenece al segmento recién calculado. El procedimiento se repite así intervalo por intervalo (Hale y Verduyn Lunel, 1993).

El bloque siguiente aplica el método al ejemplo con $a=1$ y $b=2$. Cada tramo se integra con el método DOP853 (Hairer et al., 1993) de SciPy (Virtanen et al., 2020). Su salida densa se usa como historia en el tramo siguiente. El programa estima la tasa de crecimiento a partir de los máximos de $|x|$ y la compara con la raíz dominante que da la función de Lambert.

import numpy as np
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks
from scipy.special import lambertw

a, b = 1.0, 2.0
tau_c = np.arccos(-a / b) / np.sqrt(b**2 - a**2)

def metodo_de_pasos(tau, historia, t_final):
    """Integra x'(t) = -a x(t) - b x(t - tau) de a un intervalo de largo tau."""
    t_todo, x_todo = [], []
    pasado, x0 = historia, historia(0.0)
    for k in range(int(t_final / tau)):
        t0, t1 = k * tau, (k + 1) * tau
        campo = lambda t, x, p=pasado: [-a * x[0] - b * p(t - tau)]
        sol = solve_ivp(campo, (t0, t1), [x0], method="DOP853",
                        rtol=1e-10, atol=1e-12, dense_output=True)
        pasado = lambda s, f=sol.sol: f(s)[0]   # historia para el tramo siguiente
        x0 = sol.y[0, -1]
        t = np.linspace(t0, t1, 200)
        t_todo.append(t)
        x_todo.append(sol.sol(t)[0])
    return np.concatenate(t_todo), np.concatenate(x_todo)

print(f"retardo crítico tau_c = {tau_c:.6f}")
print(" tau/tau_c   tasa simulada   Re(lambda) por W de Lambert")
for f in (0.70, 1.05):
    tau = f * tau_c
    t, x = metodo_de_pasos(tau, lambda s: 1.0, 60.0)
    picos, _ = find_peaks(np.abs(x))
    tardios = picos[t[picos] > 20]             # descarta el transitorio
    tasa = np.polyfit(t[tardios], np.log(np.abs(x[tardios])), 1)[0]
    raiz = lambertw(-b * tau * np.exp(a * tau), 0) / tau - a
    print(f"{f:9.2f}   {tasa:13.5f}   {raiz.real:13.5f}")

La ejecución produce la salida siguiente.

retardo crítico tau_c = 1.209200
 tau/tau_c   tasa simulada   Re(lambda) por W de Lambert
     0.70        -0.20978        -0.20978
     1.05         0.01805         0.01805

Las tasas simuladas coinciden con las espectrales hasta el quinto decimal. El valor de $\tau_c$ proviene de la ecuación característica y la simulación confirma que la implementación es coherente con él.

El procedimiento funciona bien en este ejemplo, con un único retardo constante y una historia suave, pero no es un integrador general. La derivada de la solución es discontinua en $t=0$ cuando la historia no empalma con la ecuación. Esa discontinuidad se propaga a derivadas de orden mayor en $t=\tau,2\tau,\ldots$ Los métodos especializados manejan esas discontinuidades, los retardos múltiples y la interpolación con control de error (Shampine y Thompson, 2001). Un integrador de ecuaciones ordinarias que no incorpora bien el término retardado puede producir una trayectoria suave y equivocada.

Estabilidad en un espacio de historias

En una ecuación con retardo se perturba una función inicial. La cercanía requiere una norma sobre historias, por ejemplo

Que el equilibrio $x^*=0$ sea estable significa que las historias suficientemente próximas a la historia nula generan trayectorias que permanecen próximas. La estabilidad asintótica añade la convergencia a cero. Por eso el retrato de fase de una ecuación escalar con retardo no es una recta. Pueden dibujarse proyecciones en coordenadas como $(x(t),x(t-\tau))$, como en la figura 48.3b, o reconstruirse atractores a partir de observaciones retardadas, aunque esas imágenes proyectan un estado funcional más rico. El teorema de Poincaré–Bendixson del capítulo 39 no se aplica por el hecho de dibujar dos coordenadas retardadas, porque la ecuación no se ha convertido en un flujo plano autónomo.

Tres intuiciones que el retardo contradice

Un retardo no implica inestabilidad. En $\dot x=-2x(t)-x(t-\tau)$ se cumple $a>b>0$ y el equilibrio es asintóticamente estable para todo $\tau\ge0$. El retardo puede reducir el margen de amortiguamiento o introducir oscilaciones transitorias sin obligar a cruzar el eje imaginario.

Una oscilación tampoco implica una fuerza externa. El modelo lineal retardado es autónomo. Cerca del umbral sus modos son oscilatorios porque las raíces dominantes son complejas, sin ninguna señal periódica que imponga la frecuencia. El capítulo 49 estudiará el caso distinto en que una entrada periódica actúa sobre el sistema.

Una memoria aparente, por último, no implica un retardo explícito. El sistema lento-rápido del capítulo 47 responde de manera distinta al mismo impulso según el valor de su variable lenta. Quien observe solo la variable rápida podría describir el efecto como memoria, aunque el modelo completo es una ecuación ordinaria bidimensional y su futuro queda determinado por el estado presente. Confundir ambos mecanismos puede llevar a postular un retardo que no existe. El capítulo 49 añadirá una tercera fuente de dependencia temporal, una entrada externa cuya fase evoluciona con independencia del sistema.

Ejercicios

Ejercicio 48.1. Considere $\dot x(t)=-0{,}5\,x(t)-x(t-\tau)$. (a) Calcule la frecuencia de cruce y el retardo crítico. (b) Calcule el período del modo crítico. (c) ¿Qué ocurre si el coeficiente del término retardado se reduce a $0{,}4$?

Ejercicio 48.2. Demuestre que, si $a\ge b>0$, la ecuación $\lambda+a+be^{-\lambda\tau}=0$ no tiene raíces sobre el eje imaginario para ningún $\tau\ge0$. Considere por separado el caso $\lambda=0$.

Ejercicio 48.3. Para $a=1$ y $b=2$, (a) calcule $\operatorname{Re}(d\lambda/d\tau)$ en el primer cruce. (b) Estime la parte real de la raíz dominante en $\tau=1{,}05\,\tau_c$ con una aproximación lineal en $\tau$ y compárela con el valor exacto, $0{,}01805$.

Ejercicio 48.4. Con $a=1$ y $b=2$, (a) calcule el retardo a partir del cual la raíz dominante deja de ser real resolviendo $2\tau e^{\tau}=1/e$. (b) Explique por qué la aproximación de retardo pequeño no puede describir el comportamiento para $\tau$ mayor que ese valor.

Ejercicio 48.5. En la logística de Hutchinson, (a) calcule $\rho=r\tau$ para una población con $r=0{,}5$ por día y retardos de 2 y de 4 días. Decida en cada caso si el equilibrio es estable. (b) Estime, a partir del cruce, el período de las oscilaciones cerca del umbral en días. (c) Explique por qué el mismo retardo de 4 días puede ser irrelevante para una población con $r=0{,}05$ por día.

Ejercicio 48.6. Aplique a mano el método de pasos a $\dot x(t)=-x(t-1)$ con la historia $\phi(s)=1$ en $[-1,0]$. (a) Obtenga $x(t)$ en $[0,1]$. (b) Obtenga $x(t)$ en $[1,2]$ y calcule $x(2)$. (c) Compare las derivadas laterales de $x$ en $t=0$ y las segundas derivadas laterales en $t=1$.

Referencias

Corless, R. M., Gonnet, G. H., Hare, D. E. G., Jeffrey, D. J. y Knuth, D. E. (1996). On the Lambert W function. Advances in Computational Mathematics, 5(1), 329–359. https://doi.org/10.1007/BF02124750

Erneux, T. (2009). Applied delay differential equations. Springer. https://doi.org/10.1007/978-0-387-74372-1

Hairer, E., Nørsett, S. P. y Wanner, G. (1993). Solving ordinary differential equations I: Nonstiff problems (2.ª ed. rev.). Springer. https://doi.org/10.1007/978-3-540-78862-1

Hale, J. K. y Verduyn Lunel, S. M. (1993). Introduction to functional differential equations. Springer. https://doi.org/10.1007/978-1-4612-4342-7

Hutchinson, G. E. (1948). Circular causal systems in ecology. Annals of the New York Academy of Sciences, 50(4), 221–246. https://doi.org/10.1111/j.1749-6632.1948.tb39854.x

Mackey, M. C. y Glass, L. (1977). Oscillation and chaos in physiological control systems. Science, 197(4300), 287–289. https://doi.org/10.1126/science.267326

Shampine, L. F. y Thompson, S. (2001). Solving DDEs in MATLAB. Applied Numerical Mathematics, 37(4), 441–458. https://doi.org/10.1016/S0168-9274(00)00055-6

Smith, H. (2011). An introduction to delay differential equations with applications to the life sciences. Springer. https://doi.org/10.1007/978-1-4419-7646-8

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., … SciPy 1.0 Contributors. (2020). SciPy 1.0: Fundamental algorithms for scientific computing in Python. Nature Methods, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2