Saltos, ciclos e histéresis en un caso resuelto

El problema

Los capítulos 41 a 49 presentaron por separado osciladores forzados, ciclos límite, bifurcaciones de equilibrios, la bifurcación de Hopf, sistemas lentos y rápidos, retardos, resonancia y sincronización. Este capítulo reúne varias de esas ideas en un solo problema, en el que aparecen juntas sin significar lo mismo. La pregunta es cómo un sistema determinista continuo puede mostrar saltos entre estados, una respuesta periódica y un bucle de histéresis sin que esos tres fenómenos sean equivalentes.

Usaremos el modelo escalar

Con $\mu$ fijo, la ecuación es autónoma y unidimensional. Puede tener uno o tres equilibrios, pero ningún ciclo límite. Si $\mu$ varía lentamente, los equilibrios estables pueden desaparecer en bifurcaciones de nodo-silla y la trayectoria realiza transiciones rápidas. Si $\mu$ es periódico, el sistema se vuelve no autónomo y su respuesta, pasado el transitorio, es periódica. Al representar esa respuesta frente a la entrada aparece un bucle de histéresis. El interés del caso está en separar lo que el bucle debe a la geometría de los equilibrios de lo que debe a la velocidad del barrido.

Estructura de gradiente

Con $f(x;\mu)=\mu+x-x^3$, el modelo es un flujo gradiente, $\dot x=-\partial V/\partial x$, con el potencial

Para $\mu$ fijo se cumple

La función $V$ decrece estrictamente fuera de los equilibrios, por lo que el sistema autónomo no puede sostener una oscilación. Después de una vuelta completa $V$ tendría que recuperar su valor inicial, aunque habría disminuido durante el recorrido. Cuando $\mu(t)$ sea periódico aparecerán respuestas periódicas, sin que el sistema a parámetro fijo esconda un ciclo límite.

Equilibrios y pliegues

Los equilibrios cumplen $\mu=x^3-x$, una forma que parametriza el diagrama de bifurcación por $x$. La estabilidad depende de $f_x=1-3x^2$. Los equilibrios con $|x|>1/\sqrt3$ son estables y los que tienen $|x|<1/\sqrt3$ son inestables. Las ramas cambian de estructura donde se cumplen a la vez $f=0$ y $f_x=0$, es decir, en $x=\pm1/\sqrt3$. Sustituyendo en $\mu=x^3-x$ se obtienen los pliegues

Para $|\mu|<\mu_c$ hay tres equilibrios, dos exteriores estables y uno central inestable. Para $|\mu|>\mu_c$ queda un solo equilibrio, estable (figura 50.1a). El modelo contiene dos bifurcaciones de nodo-silla como las del capítulo 44. Aquí la geometría global une ambos pliegues en una curva en forma de S y crea una región intermedia de biestabilidad (Strogatz, 2015, sec. 3.6).

Biestabilidad sin histéresis

Con $\mu=0$ los equilibrios son $-1$, $0$ y $1$. Los dos exteriores son estables y el central es inestable. Una trayectoria que parte de $x_0<0$ converge a $-1$ y una que parte de $x_0>0$ converge a $+1$, de modo que el equilibrio inestable separa las dos cuencas.

Eso demuestra biestabilidad, todavía no histéresis. Para hablar de histéresis hay que especificar una entrada que cambia, un protocolo de variación y una respuesta que dependa del recorrido de la entrada. La coexistencia de dos atractores favorece ciertos mecanismos histeréticos sin definir la histéresis. El capítulo 9 ya advirtió que la memoria, la dependencia de la trayectoria, la multiestabilidad y la histéresis se superponen en parte sin ser intercambiables. En este modelo la relación puede mostrarse de manera explícita, porque $\mu$ hace de entrada controlada.

El paisaje de potencial

La forma de $V$ muestra la misma estructura. Con $\mu=0$ el potencial tiene dos mínimos en $x=\pm1$, separados por un máximo en $x=0$. Las barreras de ambos pozos valen $1/4$. Al aumentar $\mu$, el término $-\mu x$ inclina el potencial, el pozo derecho se hace más profundo y el izquierdo más somero. Con $\mu=0{,}2$ los equilibrios están en $-0{,}879$, $-0{,}209$ y $1{,}088$. Las barreras valen entonces $0{,}082$ para el pozo izquierdo y $0{,}480$ para el derecho. El mínimo izquierdo deja de ser el global mucho antes de desaparecer. Persiste como mínimo local hasta $\mu=\mu_c$, donde choca con el máximo y ambos dejan de existir (figura 50.1b). Al disminuir $\mu$ ocurre lo mismo, de forma simétrica, en $-\mu_c$.

Mientras exista un mínimo local, una trayectoria determinista situada en su cuenca converge hacia él aunque otro pozo sea más profundo. La comparación de profundidades no sustituye a la dinámica ni elimina las barreras entre cuencas. En un sistema con ruido la situación cambia, porque las fluctuaciones pueden cruzar la barrera más baja antes de que la rama desaparezca.

Barridos lentos en ambos sentidos

Supongamos que $\mu$ empieza por debajo de $-\mu_c$, donde solo existe la rama izquierda, para luego aumentar lentamente. Al entrar en la región $-\mu_c<\mu<\mu_c$ aparecen otros dos equilibrios, pero la trayectoria ya está cerca de la rama izquierda. Si el parámetro varía despacio y el estado está lejos del pliegue, la trayectoria sigue esa rama. En $\mu=0$ el sistema permanece cerca de $x=-1$ aunque también exista el equilibrio estable $x=+1$. La entrada actual no determina qué estado se observa, porque importa desde qué rama llegó la trayectoria. La rama izquierda persiste hasta $\mu=\mu_c$, donde choca con la inestable. Para $\mu>\mu_c$ no queda ningún equilibrio a la izquierda y la trayectoria se dirige rápidamente a la rama derecha.

El salto no es una discontinuidad de $x(t)$. El campo es suave y la solución es continua. La palabra salto indica que $x$ cambia en una escala de tiempo corta comparada con la variación lenta de $\mu$, como en los sistemas lentos y rápidos del capítulo 47.

Al disminuir $\mu$ desde valores positivos grandes, la trayectoria sigue la rama derecha y no vuelve a la izquierda al pasar por $\mu=0$. La rama derecha desaparece recién en $-\mu_c$. Los umbrales cuasiestáticos son, por tanto, distintos, $\mu_\uparrow=\mu_c$ para el barrido ascendente y $\mu_\downarrow=-\mu_c$ para el descendente. Entre ambos valores, una misma entrada es compatible con dos respuestas estables y el recorrido previo decide cuál se observa. Esa es la histéresis en sentido dinámico. La relación entre entrada y salida no puede describirse con una función $x=h(\mu)$ independiente de la historia, porque hay dos ramas estables accesibles y los cambios entre ellas ocurren en umbrales que dependen de la dirección del barrido.

El área cuasiestática del bucle

En el límite de un barrido arbitrariamente lento, el bucle está limitado por las dos ramas estables entre $-\mu_c$ y $\mu_c$. Su área es

donde $x_R$ y $x_L$ son las ramas derecha e izquierda. Como $d\mu=(3x^2-1)\,dx$ y las ramas son simétricas, basta integrar la rama derecha y duplicar. Cuando $\mu$ va de $-\mu_c$ a $\mu_c$, esa rama va de $x=1/\sqrt3$ a $x=2/\sqrt3$, de modo que

El valor depende de la normalización del modelo, como muestra el ejercicio 50.3. El procedimiento sí se generaliza. Antes de interpretar un bucle hay que identificar las ramas estables y los umbrales que corresponden al protocolo aplicado.

Forzamiento periódico y respuesta forzada

Con $\mu(t)=A\sin(\Omega t)$ y $A>\mu_c$, la entrada supera ambos pliegues y obliga al sistema a abandonar alternativamente cada rama. Tomaremos $A=0{,}8$. La ecuación ya no es autónoma, porque dos instantes con el mismo $x$ pueden tener derivadas distintas si la fase de la entrada difiere. Con la fase $\phi=\Omega t\pmod{2\pi}$ se obtiene el sistema autónomo ampliado

cuyo espacio de estados es un cilindro. Pasado el transitorio, la respuesta es periódica con el período de la entrada. En la formulación ampliada corresponde a una órbita cerrada, que no es el ciclo límite autosostenido del capítulo 42, porque sin forzamiento el sistema vuelve a un equilibrio.

Con $\Omega$ pequeño, la trayectoria sigue durante la mayor parte del ciclo una de las ramas estables instantáneas (figura 50.2a). Sigue lentamente la rama izquierda mientras $\mu$ aumenta, salta hacia la derecha cerca de $\mu_c$, sigue lentamente la rama derecha mientras $\mu$ disminuye y salta de vuelta cerca de $-\mu_c$. El patrón recuerda una oscilación de relajación, aunque el reloj que obliga a recorrer el ciclo está fuera del subsistema escalar. La fase $\phi$ incorpora ese reloj al estado sin convertirlo en una variable interna generada por $x$.

Retraso del salto a velocidad finita

En el límite cuasiestático el salto ascendente ocurriría en $\mu_c$. A velocidad finita, la trayectoria permanece un tiempo cerca de la región donde desapareció el equilibrio, el retraso dinámico estudiado en el capítulo 43 (Kuehn, 2015). Cerca del pliegue ascendente, con $x=-1/\sqrt3+y$ y $\mu=\mu_c+\delta\mu$, el campo es

El término lineal en $y$ desaparece porque el punto es un pliegue. Si el parámetro cruza el umbral con velocidad $r$, de modo que $\delta\mu\approx rt$, la ecuación dominante es $\dot y=rt+\sqrt3\,y^2$. Con las escalas $t=(\sqrt3\,r)^{-1/3}\,T$ e $y=(\sqrt3\,r)^{1/3}\,Y/\sqrt3$ se reduce a $dY/dT=T+Y^2$, una ecuación de Riccati cuya solución escapa a infinito cuando $T$ alcanza el primer cero de la función de Airy, $a_1\approx2{,}338$. El retraso predicho es

Con el forzamiento sinusoidal, la velocidad en el pliegue es $r=\Omega\sqrt{A^2-\mu_c^2}$. Ya en régimen periódico, el bloque siguiente mide el valor de $\mu$ en el que $x$ cruza el cero durante el salto ascendente y lo compara con la predicción. Calcula también el área del bucle como $-\oint x\,d\mu$ sobre un período. La integración usa el método DOP853 (Hairer et al., 1993) de SciPy (Virtanen et al., 2020).

import numpy as np
from scipy.integrate import solve_ivp
from scipy.special import ai_zeros

A = 0.8
mu_c = 2 / (3 * np.sqrt(3))
a1 = -ai_zeros(1)[0][0]                      # primer cero de Ai, 2.33811

cruce = lambda t, x, Om: x[0]                # x = 0 recorrido hacia arriba
cruce.direction = 1

print("    Omega   mu del salto   retraso   predicción   área del bucle")
for Om in (0.0005, 0.002, 0.02, 0.1):
    T = 2 * np.pi / Om
    campo = lambda t, x, Om: [A * np.sin(Om * t) + x[0] - x[0]**3]
    sol = solve_ivp(campo, (0, 3 * T), [-1.2], args=(Om,), events=cruce,
                    method="DOP853", rtol=1e-10, atol=1e-12, dense_output=True)
    t_salto = sol.t_events[0][-1]             # último salto ascendente
    mu_salto = A * np.sin(Om * t_salto)
    r = Om * np.sqrt(A**2 - mu_c**2)          # velocidad de mu en el pliegue
    prediccion = a1 * 3**(-1 / 6) * r**(2 / 3)
    t = np.linspace(2 * T, 3 * T, 200001)     # tercer período, ya estacionario
    x = sol.sol(t)[0]
    area = -np.trapezoid(x * A * Om * np.cos(Om * t), t)   # -(integral de x dmu)
    print(f"{Om:9.4f}   {mu_salto:10.4f}   {mu_salto - mu_c:7.4f}"
          f"   {prediccion:9.4f}   {area:10.4f}")

El programa imprime lo siguiente.

    Omega   mu del salto   retraso   predicción   área del bucle
   0.0005       0.3944    0.0095      0.0097       1.5335
   0.0020       0.4083    0.0234      0.0244       1.5838
   0.0200       0.4838    0.0989      0.1132       1.8736
   0.1000       0.6298    0.2449      0.3311       2.4819

La predicción es asintótica. Con $\Omega=0{,}0005$ difiere del retraso medido en un 2 %. Con $\Omega=0{,}0002$ el cociente entre ambos llega a $0{,}99$. Con frecuencias mayores la predicción sobrestima el retraso, porque la velocidad del barrido deja de ser constante durante el paso y el término cúbico se hace apreciable (figura 50.3a). Los valores medidos no son nuevos puntos de bifurcación. Los pliegues del sistema congelado siguen en $\pm0{,}3849$. Lo que cambia es el momento en que una trayectoria sometida a un parámetro variable abandona la vecindad de la rama desaparecida.

El área también crece con la frecuencia. El exceso sobre $3/2$ proviene casi por completo de los retrasos. Cada salto se desplaza en $\delta\mu$ y la distancia vertical entre las ramas en el pliegue es $2/\sqrt3+1/\sqrt3=\sqrt3$, por lo que el exceso es aproximadamente $2\sqrt3\,\delta\mu$. Con $\Omega=0{,}002$ esa estimación da $0{,}081$, frente al exceso medido de $0{,}084$.

Balance de disipación

Cuando $\mu$ depende del tiempo, $V(x;\mu(t))$ ya no es monótona. Su derivada total es

El primer término es disipativo y el segundo es el intercambio debido al cambio externo del parámetro. En una respuesta periódica $V$ vuelve a su valor después de un período, así que

El área orientada del bucle mide la disipación acumulada en cada ciclo. La igualdad es propia de este modelo de gradiente y no define la histéresis en general. Con $\Omega=0{,}02$ ambos miembros valen $1{,}873638$ y coinciden hasta el duodécimo decimal, lo que comprueba la integración de manera independiente.

Histéresis independiente y dependiente de la velocidad

El modelo separa dos componentes que suelen aparecer mezcladas. Aun cuando la frecuencia tiende a cero, los umbrales ascendente y descendente siguen separados y el área tiende a $3/2$. Esa parte del bucle se debe a la biestabilidad y a la desaparición de las ramas en pliegues distintos. A frecuencia finita se suma el retraso respecto de los pliegues, que agranda el bucle.

Un sistema monoestable con relajación lenta también produce un bucle a frecuencia finita, por el desfase entre entrada y salida, aunque ese bucle colapsa cuando el barrido se hace muy lento. Con $\dot x=\mu(t)-x$ y el mismo forzamiento, el área es $\pi A^2\Omega/(1+\Omega^2)$, que vale $0{,}040$ con $\Omega=0{,}02$ y tiende a cero con $\Omega$ (figura 50.3b). Un área no nula a una sola frecuencia no distingue los dos mecanismos. Hace falta estudiar cómo cambia el área al reducir la velocidad del barrido.

Tres estructuras cerradas que no son un ciclo límite

El caso contiene tres curvas cerradas que podrían llamarse ciclos. El bucle entrada-salida en el plano $(\mu,x)$ representa la histéresis y sus ejes combinan una entrada externa con una variable del sistema. La respuesta temporal periódica cumple $x(t+T)=x(t)$ con $T=2\pi/\Omega$ después del transitorio. La órbita periódica del sistema ampliado en el cilindro $(x,\phi)$ permite usar el lenguaje de los sistemas autónomos, aunque la coordenada $\phi$ sea la fase de la entrada. Ninguna demuestra que el subsistema $\dot x=\mu+x-x^3$ tenga un ciclo límite, algo imposible en un sistema escalar continuo que, además, tiene el potencial decreciente $V$ cuando $\mu$ es fijo.

Dónde reside la dependencia de la historia

En el sistema no autónomo, el valor actual de $x$ junto con la fase $\phi$ determina el futuro. No hace falta ninguna variable de memoria adicional ni un retardo explícito. La histéresis aparece porque, dentro de la región biestable, historias distintas de la entrada conducen a estados actuales distintos. Quien observe solo $\mu$ tiene un estado subespecificado, ya que el mismo valor de $\mu$ puede acompañarse de un $x$ en cualquiera de las dos ramas. En $\mu=0$, el barrido ascendente lento encuentra al sistema cerca de $-1$ y el descendente cerca de $+1$. La información histórica está en el valor actual de $x$. El caso contrasta con el del capítulo 48, donde el valor puntual $x(t)$ podía ser insuficiente aun conociendo el parámetro y el tiempo, porque la ecuación requería un segmento de historia.

Salto por pérdida de estabilidad y salto por perturbación

El capítulo 44 distinguió dos maneras de cambiar de estado. Ambas caben en este diagrama. Con $\mu=0$ existen las dos ramas estables. Una perturbación que lleve $x$ de un valor negativo a uno positivo cruza la frontera en $x=0$ y el sistema termina en el otro atractor sin que cambie ningún parámetro ni desaparezca ningún equilibrio. Durante el barrido ascendente, en cambio, la rama izquierda desaparece en $\mu_c$ mediante una nodo-silla y el cambio pertenece a la familia parametrizada.

Las dos trayectorias muestran una transición rápida de valores negativos a positivos. Una serie temporal de $x$ no basta para saber cuál de los dos mecanismos ocurrió. Distinguirlos requiere información sobre la entrada, los equilibrios disponibles, las cuencas y el protocolo aplicado. Este caso cierra el recorrido iniciado en el capítulo 41.

Ejercicios

Ejercicio 50.1. Con $\mu=0{,}2$, (a) encuentre los tres equilibrios y clasifíquelos. (b) Si $\mu$ llegó a $0{,}2$ aumentando lentamente desde $-1$, ¿cerca de qué equilibrio está el sistema? (c) ¿Y si llegó disminuyendo desde $1$?

Ejercicio 50.2. Con $\mu=0$ y con $\mu=0{,}2$, calcule la altura de la barrera de potencial de cada pozo, es decir, la diferencia entre $V$ en el equilibrio inestable y en cada mínimo. Explique por qué el sistema determinista permanece en el pozo más somero.

Ejercicio 50.3. Considere $\dot x=\mu+ax-x^3$ con $a>0$. (a) Encuentre los pliegues. (b) Con el cambio $x=\sqrt a\,X$, demuestre que el área cuasiestática del bucle es $\tfrac32a^2$. (c) ¿Qué indica este resultado sobre el significado del número $3/2$?

Ejercicio 50.4. (a) Derive la expresión $dV/dt=-\dot x^2-x\dot\mu$. (b) Deduzca la igualdad entre la disipación por ciclo y $-\oint x\,d\mu$. (c) Explique por qué el bucle de la figura 50.2b se recorre en sentido antihorario y por qué eso hace positiva el área $-\oint x\,d\mu$.

Ejercicio 50.5. Para el sistema monoestable $\dot x=\mu(t)-x$ con $\mu=A\sin(\Omega t)$, (a) obtenga la respuesta estacionaria. (b) Demuestre que el área del bucle es $\pi A^2\Omega/(1+\Omega^2)$. (c) Compare con el modelo biestable cuando $\Omega\to0$.

Ejercicio 50.6. (a) Calcule el retraso predicho del salto ascendente con $A=0{,}8$ y $\Omega=0{,}002$. (b) Repita el cálculo con $A=1{,}2$ y la misma frecuencia. (c) Explique por qué una amplitud mayor produce un retraso mayor aunque el período no cambie.

Referencias

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

Kuehn, C. (2015). Multiple time scale dynamics. Springer. https://doi.org/10.1007/978-3-319-12316-5

Strogatz, S. H. (2015). Nonlinear dynamics and chaos: With applications to physics, biology, chemistry, and engineering (2.ª ed.). Westview Press.

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