Año 2 · Módulo 11 · Nivel C

109.1. Cuando el cambio decisivo ocurre lejos del equilibrio

Las separatrices de un oscilador con dos pozos de potencial separan movimientos de distinta naturaleza. En el sistema sin forzamiento pueden coincidir con conexiones que salen de un punto de silla y regresan a él. Una perturbación periódica pequeña deforma esas conexiones y, según sus parámetros, las variedades estable e inestable pueden dejar de coincidir y cruzarse transversalmente. Ese cruce reorganiza las trayectorias en una región que la estabilidad lineal del equilibrio no describe.

El análisis local describe cómo se comporta una perturbación junto a un equilibrio o a una órbita periódica. Las conexiones homoclínicas, en cambio, exigen seguir trayectorias que recorren una zona amplia del espacio de fases antes de acercarse de nuevo al objeto invariante; su presencia depende de una relación global entre variedades. Por ello, los valores propios no deciden si la separatriz se ha escindido ni si existen intersecciones transversales.

El método de Melnikov es una herramienta perturbativa para ese problema. Parte de un sistema cuya separatriz conocemos, añade una perturbación periódica pequeña y calcula una función que aproxima, al primer orden, la separación orientada entre las variedades estable e inestable. Bajo las hipótesis apropiadas, un cero simple de esa función demuestra una intersección transversal para perturbaciones suficientemente pequeñas. Por esa vía se llega analíticamente a una estructura homoclínica compleja, aunque eso no certifica un atractor caótico visible en cualquier simulación.

Trabajaremos con un oscilador de Duffing de doble pozo porque permite hacer casi todos los cálculos a mano. Interesa entender de dónde surge la integral, qué información geométrica conserva y qué inferencias quedan por justificar después de calcularla.

109.2. Dos formas de acercarse a una silla

Sea un sistema autónomo suave en el plano, z˙=f(z)\dot z=f(z), con un equilibrio hiperbólico de tipo silla pp. Cerca de esa silla existe una dirección a lo largo de la cual las perturbaciones disminuyen al avanzar el tiempo y otra en la que crecen. La variedad estable Ws(p)W^s(p) reúne los estados cuyas trayectorias se aproximan a pp cuando t→+∞t\to+\infty; la variedad inestable Wu(p)W^u(p), los que se aproximan a pp cuando t→−∞t\to-\infty.

Una órbita zh(t)z_h(t) es homoclínica si se acerca al mismo objeto invariante en ambos límites temporales. Para una silla fija, esto significa

lim⁡t→−∞zh(t)=p,lim⁡t→+∞zh(t)=p.\lim_{t\to-\infty}z_h(t)=p, \qquad \lim_{t\to+\infty}z_h(t)=p.

No se trata de una órbita periódica, porque tarda un tiempo infinito en acercarse a la silla por cualquiera de los extremos; el lazo con que se representa una órbita homoclínica oculta esa aproximación asintótica en ambas direcciones temporales. Una conexión heteroclínica, en cambio, une dos objetos invariantes distintos. Esa diferencia importa al elegir el modelo y al definir la función de Melnikov.

En el plano autónomo, las trayectorias de un campo suave con unicidad no pueden cruzarse transversalmente entre sí, de modo que una conexión homoclínica planar todavía no tiene los cruces transversales que buscaremos. La situación cambia al añadir forzamiento periódico. Un sistema plano no autónomo equivale a un sistema autónomo de tres dimensiones si se incorpora la fase de la excitación. En una sección estroboscópica, las variedades de una órbita periódica de tipo silla pueden entonces encontrarse transversalmente.

109.3. Del campo dependiente del tiempo al mapa de retorno

Consideremos un sistema perturbado,

z˙=f0(z)+εg(z,t),g(z,t+T)=g(z,t),\dot z=f_0(z)+\varepsilon g(z,t), \qquad g(z,t+T)=g(z,t),

con 0<ε≪10<\varepsilon\ll1. Si T=2π/ωT=2\pi/\omega, añadimos una coordenada de fase θ\theta que satisface θ˙=ω\dot\theta=\omega, entendida módulo 2π2\pi. El sistema ampliado es autónomo en (z,θ)(z,\theta). La antigua silla puede dar origen a una órbita periódica hiperbólica del espacio ampliado; en nuestro ejemplo, esto ocurre para perturbaciones suficientemente pequeñas.

Otra descripción de la dinámica se obtiene muestreando en una fase fija del forzamiento. Si medimos zz en los instantes t0+nTt_0+nT, obtenemos un mapa de Poincaré, que aquí llamaremos mapa estroboscópico. Sus puntos fijos corresponden a movimientos de período TT del sistema forzado; las curvas estable e inestable de un punto fijo de tipo silla organizan los retornos cercanos.

Distinguir los espacios evita un error frecuente. En el mapa bidimensional puede haber intersecciones transversales de curvas invariantes porque esas curvas son conjuntos invariantes de una transformación discreta, a los que no se aplica la restricción de unicidad de las trayectorias del plano. El mismo fenómeno puede estudiarse con superficies invariantes en el espacio ampliado de tres dimensiones.

La pregunta perturbativa es si, al encender una fuerza débil, las ramas estable e inestable que coincidían sobre la antigua separatriz quedan separadas de manera permanente, se tocan tangencialmente o se intersectan transversalmente. Para responderla usaremos una cantidad conservada del sistema sin perturbar.

109.4. El oscilador de Duffing de doble pozo

Tomemos el sistema no perturbado

x˙=y,y˙=x−x3.\dot x=y,\qquad \dot y=x-x^3.

Equivale a la ecuación de segundo orden x¨−x+x3=0\ddot x-x+x^3=0. Su energía es

H(x,y)=y22−x22+x44.H(x,y)=\frac{y^2}{2}-\frac{x^2}{2}+\frac{x^4}{4}.

La derivada temporal satisface H˙=Hxx˙+Hyy˙\dot H=H_x\dot x+H_y\dot y. Como Hx=−x+x3H_x=-x+x^3 y Hy=yH_y=y, obtenemos H˙=(−x+x3)y+y(x−x3)=0\dot H=(-x+x^3)y+y(x-x^3)=0. La energía es constante sobre cada trayectoria del sistema sin forzar.

Los equilibrios son (0,0)(0,0), (1,0)(1,0) y (−1,0)(-1,0). La matriz jacobiana es

Df0(x,y)=(011−3x20).Df_0(x,y)=\begin{pmatrix}0&1\\1-3x^2&0\end{pmatrix}.

En el origen, los valores propios son +1+1 y −1-1, de modo que es una silla hiperbólica. En x=±1x=\pm1 son ±i2\pm i\sqrt2, que en este sistema hamiltoniano de un grado de libertad corresponden a centros. La energía en el origen es cero y en cada mínimo del potencial vale −1/4-1/4. El potencial V(x)=−x2/2+x4/4V(x)=-x^2/2+x^4/4 tiene dos pozos, separados por la barrera situada en el origen.

La curva de energía H=0H=0 contiene la silla y dos lazos homoclínicos, uno hacia x>0x>0 y otro hacia x<0x<0. Estas curvas son las separatrices, que distinguen las oscilaciones confinadas a un pozo de los movimientos que pueden recorrer ambos. La forma de la separatriz se obtiene resolviendo y2=x2−x4/2y^2=x^2-x^4/2, pero aquí tenemos además una parametrización explícita, que permitirá calcular la perturbación de sus ramas invariantes.

109.5. Una órbita homoclínica que podemos verificar

Para el lazo situado a la derecha, consideremos

xh(t)=2sech⁡t,yh(t)=−2sech⁡ttanh⁡t.x_h(t)=\sqrt2\operatorname{sech}t, \qquad y_h(t)=-\sqrt2\operatorname{sech}t\tanh t.

Como sech⁡t=1/cosh⁡t\operatorname{sech}t=1/\cosh t, su derivada es −sech⁡ttanh⁡t-\operatorname{sech}t\tanh t, de modo que x˙h=yh\dot x_h=y_h, como exige la primera ecuación. Para comprobar la segunda derivamos otra vez y obtenemos

x¨h=2sech⁡t(2tanh⁡2t−1).\ddot x_h=\sqrt2\operatorname{sech}t\left(2\tanh^2t-1\right).

Con tanh⁡2t=1−sech⁡2t\tanh^2t=1-\operatorname{sech}^2t queda x¨h=2sech⁡t−22sech⁡3t=xh−xh3\ddot x_h=\sqrt2\operatorname{sech}t-2\sqrt2\operatorname{sech}^3t=x_h-x_h^3. Hemos verificado la ecuación sin recurrir a un gráfico. Cuando t→±∞t\to\pm\infty, xh→0x_h\to0 e yh→0y_h\to0; en t=0t=0 el punto de retorno es (2,0)(\sqrt2,0).

La rama izquierda se obtiene cambiando el signo de ambas coordenadas. La elección de rama afecta al signo de algunos términos de la función de Melnikov, aunque no al umbral de existencia de ceros simples bajo forzamiento armónico simétrico. El origen temporal también interviene, porque desplazar tt cambia la fase sin alterar la geometría de la conexión.

En este ejemplo no hace falta aproximar numéricamente la separatriz antes de integrar. Tener una solución cerrada es una comodidad pedagógica que el método no exige; cuando la órbita homoclínica se conoce solo de forma numérica, las integrales de Melnikov pueden evaluarse controlando con cuidado el error y el comportamiento de las colas temporales.

C109 Figura 01 Duffing homoclinica
Figura 109.1. C109 Figura 01 Duffing homoclinica Ampliar figura ↗

Figura 109.1. El sistema conservativo x˙=y\dot x=y, y˙=x−x3\dot y=x-x^3 posee dos lazos homoclínicos sobre el nivel de energía H=0H=0. Se representan el potencial, las curvas de energía y el recorrido de la solución exacta xh(t)=2 sech(t)x_h(t)=\sqrt{2}\,\mathrm{sech}(t), yh(t)=−2 sech(t)tanh⁡(t)y_h(t)=-\sqrt{2}\,\mathrm{sech}(t)\tanh(t). Los centros y la silla se distinguen de las órbitas periódicas de los pozos.

109.6. Encender una perturbación pequeña

Con amortiguamiento lineal y una fuerza periódica, ambos de orden ε\varepsilon, la ecuación queda

x¨−x+x3=ε(Fcos⁡(ωt)−γx˙).\boxed{\ddot x-x+x^3=\varepsilon\bigl(F\cos(\omega t)-\gamma\dot x\bigr).}

Usaremos F≥0F\ge0, γ≥0\gamma\ge0, ω>0\omega>0 y 0<ε≪10<\varepsilon\ll1. El parámetro pequeño multiplica a toda la perturbación; FF y γ\gamma son coeficientes de orden uno en esta parametrización. Si se escriben directamente un coeficiente de amortiguamiento físico dd y una amplitud física aa, la identificación es d=εγd=\varepsilon\gamma y a=εFa=\varepsilon F; mezclar las dos convenciones conduce a umbrales incorrectos.

El sistema forzado satisface

x˙=y,y˙=x−x3+ε(Fcos⁡(ωt)−γy).\dot x=y,\qquad \dot y=x-x^3+\varepsilon\bigl(F\cos(\omega t)-\gamma y\bigr).

La energía anterior deja de conservarse; al derivarla sobre la trayectoria forzada se obtiene la identidad exacta

dHdt=εyFcos⁡(ωt)−εγy2.\frac{dH}{dt}=\varepsilon yF\cos(\omega t)-\varepsilon\gamma y^2.

El primer término es el trabajo instantáneo de la fuerza externa y puede ser positivo o negativo según la fase. El segundo nunca es positivo, porque el amortiguamiento extrae energía. Que a lo largo de una excursión completa se recupere tanta energía como se disipa depende, entonces, de la fase de la excitación. La función de Melnikov expresa esa dependencia.

109.7. Construir la función sin perder su significado

Imaginemos que la órbita forzada permanece próxima al antiguo lazo homoclínico durante una excursión completa. Como primera aproximación, evaluamos la variación de energía sobre la solución no perturbada (xh(t),yh(t))(x_h(t),y_h(t)). Si la fase de la fuerza se desplaza en t0t_0, el balance al primer orden queda

M(t0)=∫−∞∞yh(t)[Fcos⁡(ω(t+t0))−γyh(t)]dt.\boxed{M(t_0)=\int_{-\infty}^{\infty}y_h(t) \bigl[F\cos\bigl(\omega(t+t_0)\bigr)-\gamma y_h(t)\bigr]dt.}

La integral converge porque la órbita homoclínica se aproxima exponencialmente a la silla, de modo que yh(t)y_h(t) decae cuando ∣t∣→∞|t|\to\infty. En problemas más generales, donde las tasas asintóticas pueden ser distintas, esta convergencia debe comprobarse. El desplazamiento t0t_0 es la fase relativa entre el recorrido homoclínico y el forzamiento; durante la integración permanece fijo.

En una sección transversal adecuada, la separación con signo entre las variedades perturbadas admite una expansión proporcional a εM(t0)+O(ε2)\varepsilon M(t_0)+O(\varepsilon^2). El factor de proporcionalidad depende de la sección y de la forma de medir la separación. Por eso MM no es literalmente una distancia euclídea con unidades universales. Lo que cuenta son sus ceros simples, que se conservan al multiplicar por factores suaves no nulos.

El cálculo combina una trayectoria conocida y una perturbación, sin integrar la dinámica completa para todas las condiciones iniciales. Así extrae información global de una expansión perturbativa local en el parámetro ε\varepsilon, evaluada a lo largo de una órbita que recorre una región amplia del espacio de fases.

109.8. Dos integrales, dos efectos físicos

Podemos separar el balance en

M(t0)=FIforz(t0)−γIdis,M(t_0)=F I_{\mathrm{forz}}(t_0)-\gamma I_{\mathrm{dis}},

con Iforz(t0)=∫yh(t)cos⁡(ω(t+t0))dtI_{\mathrm{forz}}(t_0)=\int y_h(t)\cos(\omega(t+t_0))dt e Idis=∫yh2(t)dtI_{\mathrm{dis}}=\int y_h^2(t)dt, integrando siempre en toda la recta real. La integral disipativa es positiva, mientras que la contribución forzada oscila con la fase. Si el forzamiento es demasiado débil, el balance de primer orden queda con un signo constante; si supera un umbral, pueden aparecer fases con balances de signos opuestos.

El lenguaje energético ayuda a interpretar el cálculo, pero hay que usarlo con precisión. Lo que se calcula es la corrección de primer orden a lo largo de una trayectoria de referencia; ni la trayectoria de energía cero sobrevive intacta al forzamiento ni la energía se conserva sobre la órbita perturbada. El teorema geométrico relaciona después esa corrección con las variedades invariantes perturbadas.

En el caso de una perturbación hamiltoniana, donde no hay amortiguamiento, la función de Melnikov puede cambiar de signo sin un término disipativo constante. En sistemas que no tienen una energía natural se usa una formulación más general basada en determinantes o productos exteriores de campos, con el factor de integración adecuado. La derivación mediante H˙\dot H es especialmente transparente para sistemas hamiltonianos planos perturbados como el que hemos elegido.

109.9. Calcular la disipación exactamente

Para nuestra separatriz,

yh2(t)=2sech⁡2ttanh⁡2t.y_h^2(t)=2\operatorname{sech}^2t\tanh^2t.

Introduzcamos u=tanh⁡tu=\tanh t. Entonces du=sech⁡2t dtdu=\operatorname{sech}^2t\,dt, con límites −1-1 y +1+1, por lo que

Idis=2∫−11u2du=2[u33]−11=43.I_{\mathrm{dis}}=2\int_{-1}^{1}u^2du =2\left[\frac{u^3}{3}\right]_{-1}^{1}=\frac43.

Para la separatriz normalizada que hemos usado, la contribución del amortiguamiento a la función de Melnikov es exactamente −4γ/3-4\gamma/3. El coeficiente cambia al modificar la escala temporal, el término cúbico o la amplitud de la órbita homoclínica. Cuando se comparan expresiones publicadas para versiones distintas de Duffing, hay que reconstruir primero las transformaciones de variables y parámetros que conectan los modelos.

El término disipativo no depende de t0t_0, porque al recorrer el mismo lazo no perturbado el amortiguamiento extrae la misma cantidad integrada de energía, sea cual sea la fase de la fuerza externa. Así, MM será un seno desplazado verticalmente por una constante negativa; la existencia de ceros dependerá de cuál de los dos términos domine.

109.10. Calcular la contribución armónica

Con la identidad trigonométrica del coseno de una suma,

Iforz(t0)=cos⁡(ωt0)∫yh(t)cos⁡(ωt)dt−sin⁡(ωt0)∫yh(t)sin⁡(ωt)dt.I_{\mathrm{forz}}(t_0)=\cos(\omega t_0)\int y_h(t)\cos(\omega t)dt -\sin(\omega t_0)\int y_h(t)\sin(\omega t)dt.

Como yh(t)y_h(t) es impar y cos⁡(ωt)\cos(\omega t) es par, la primera integral vale cero. Para la segunda, observemos que yh(t)=2 d(sech⁡t)/dty_h(t)=\sqrt2\,d(\operatorname{sech}t)/dt; una integración por partes, cuyo término de borde desaparece por el decaimiento exponencial, da

∫−∞∞yh(t)sin⁡(ωt)dt=−2 ω∫−∞∞sech⁡tcos⁡(ωt)dt.\int_{-\infty}^{\infty}y_h(t)\sin(\omega t)dt =-\sqrt2\,\omega\int_{-\infty}^{\infty}\operatorname{sech}t\cos(\omega t)dt.

La transformada de Fourier conocida de sech⁡t\operatorname{sech}t da

∫−∞∞sech⁡tcos⁡(ωt)dt=πcosh⁡(πω/2).\int_{-\infty}^{\infty}\operatorname{sech}t\cos(\omega t)dt =\frac{\pi}{\cosh(\pi\omega/2)}.

En consecuencia,

Iforz(t0)=2πωcosh⁡(πω/2)sin⁡(ωt0).I_{\mathrm{forz}}(t_0)= \frac{\sqrt2\pi\omega}{\cosh(\pi\omega/2)}\sin(\omega t_0).

Con la otra rama homoclínica cambiaría el signo de la contribución forzada, aunque el umbral absoluto sería el mismo por la simetría del potencial y el carácter aditivo del forzamiento. En un potencial asimétrico, las dos conexiones pueden tener integrales y umbrales distintos.

109.11. El umbral de intersección y el sentido de un cero simple

Al juntar ambos términos se obtiene

M(t0)=2πFωcosh⁡(πω/2)sin⁡(ωt0)−4γ3.\boxed{M(t_0)=\frac{\sqrt2\pi F\omega}{\cosh(\pi\omega/2)}\sin(\omega t_0)-\frac{4\gamma}{3}.}

Definamos A=2πFω/cosh⁡(πω/2)A=\sqrt2\pi F\omega/\cosh(\pi\omega/2) y B=4γ/3B=4\gamma/3. Para F,γ≥0F,\gamma\ge0, tenemos M=Asin⁡(ωt0)−BM=A\sin(\omega t_0)-B. Si A>BA>B, hay dos ceros simples por período de fase y en cada uno la derivada M′(t0)=Aωcos⁡(ωt0)M'(t_0)=A\omega\cos(\omega t_0) no se anula. El criterio queda

F>Fcrit=4γcosh⁡(πω/2)32πω.\boxed{F>F_{\mathrm{crit}}= \frac{4\gamma\cosh(\pi\omega/2)}{3\sqrt2\pi\omega}.}

Bajo las condiciones habituales de suavidad, hiperbolicidad y perturbación suficientemente pequeña, estos ceros simples garantizan intersecciones transversales de las variedades estable e inestable de la órbita periódica hiperbólica perturbada. La conclusión es perturbativa y no indica hasta qué valor de ε\varepsilon sigue siendo válida.

Si A=BA=B, la función toca cero en el máximo del seno, donde la derivada también se anula. El primer orden no decide entonces si la intersección real es tangencial, transversal o inexistente; podrían decidirlo términos de orden superior. Si A<BA<B, la aproximación de primer orden no tiene ceros y no garantiza transversalidad; otras estructuras dinámicas o perturbaciones finitas requieren un análisis aparte. El criterio da una condición suficiente local, sin clasificar toda la dinámica perturbada.

109.12. Resolver un caso con números

Elijamos ω=1\omega=1 y γ=0,1\gamma=0,1. El factor que multiplica a FF es

2πcosh⁡(π/2)≈1,7706524.\frac{\sqrt2\pi}{\cosh(\pi/2)}\approx1,7706524.

El término disipativo vale 4γ/3=0,13333334\gamma/3=0,1333333 y la amplitud umbral resulta Fcrit≈0,075302F_{\mathrm{crit}}\approx0,075302. Con F=0,05F=0,05, la contribución sinusoidal tiene amplitud aproximada 0,088530,08853, menor que la pérdida por fricción; por tanto, M(t0)<0M(t_0)<0 para toda fase y no se cumple el criterio de Melnikov de transversalidad. Con F=0,10F=0,10, la amplitud aumenta a 0,177070,17707, excede 0,133330,13333 y aparecen ceros simples.

En el segundo caso, el seno en un cero debe satisfacer sin⁡t0=0,1333333/0,1770652≈0,753\sin t_0=0,1333333/0,1770652\approx0,753. Hay dos soluciones por período, una con coseno positivo y otra con coseno negativo, de modo que sus pendientes tienen signos opuestos, como corresponde al intercambio de posición relativa de las ramas invariantes al variar la fase.

La desigualdad final no contiene explícitamente el parámetro de perturbación ε\varepsilon, porque ese factor común se ha extraído al calcular el primer orden, pero el forzamiento y el amortiguamiento físicos siguen siendo proporcionales a ε\varepsilon. Que el símbolo no aparezca en el umbral no implica que este valga para cualquier intensidad; fuera del régimen perturbativo, la frontera debe contrastarse con el sistema completo.

C109 Figura 02 Umbral Melnikov
Figura 109.2. C109 Figura 02 Umbral Melnikov Ampliar figura ↗

Figura 109.2. Función de Melnikov M(ϕ)=2πFωsin⁡(ϕ)/cosh⁡(πω/2)−4γ/3M(\phi)=\sqrt{2}\pi F\omega\sin(\phi)/\cosh(\pi\omega/2)-4\gamma/3 para γ=0,1\gamma=0{,}1 y ω=1\omega=1, con ϕ=ωt0\phi=\omega t_0. Para F=0,05F=0{,}05 no hay ceros; en F=Fcrit≈0,07530F=F_{\mathrm{crit}}\approx0{,}07530 aparece un cero no simple, insuficiente para concluir transversalidad al primer orden; con F=0,10F=0{,}10 hay dos ceros simples por período. La curva de la derecha muestra cómo cambia el umbral con la frecuencia. La conclusión de Melnikov requiere una perturbación suficientemente pequeña.

109.13. La frecuencia también selecciona qué perturbaciones importan

La amplitud de la contribución forzada contiene el factor

K(ω)=2πωcosh⁡(πω/2).K(\omega)=\frac{\sqrt2\pi\omega}{\cosh(\pi\omega/2)}.

Cuando ω→0+\omega\to0^+, el numerador tiende a cero y el denominador a uno, por lo que K(ω)∼2πωK(\omega)\sim\sqrt2\pi\omega. Una fuerza casi constante durante el tránsito por el lazo hace contribuciones positivas y negativas que casi se cancelan en el balance de primer orden. En el extremo ω→∞\omega\to\infty, el denominador crece exponencialmente y KK vuelve a ser pequeño; la oscilación demasiado rápida se cancela por promediado a lo largo del recorrido homoclínico.

Entre ambos extremos hay una franja de frecuencias con mayor contribución integrada. Para γ\gamma fijo, el umbral FcritF_{\mathrm{crit}} es inversamente proporcional a K(ω)K(\omega), de modo que es alto para frecuencias muy pequeñas y muy grandes y menor en una región intermedia. Este efecto de integración es distinto de la resonancia lineal ordinaria de un oscilador alrededor de un mínimo de potencial. La escala temporal que interviene aquí es la de la excursión próxima a la separatriz y su permanencia cerca de la silla.

Así, una perturbación no se califica de «pequeña» o «grande» solo por su amplitud instantánea, porque también cuenta cómo se acopla con el recorrido temporal del sistema. El cálculo de Melnikov cuantifica parte de esa dependencia mientras el sistema conserve la estructura perturbativa supuesta.

109.14. Consecuencias de la intersección transversal

Una intersección homoclínica transversal de las variedades invariantes de un punto fijo hiperbólico de un difeomorfismo bidimensional conduce, bajo las hipótesis correspondientes, a una estructura de herradura de Smale para alguna iteración del mapa. Una herradura contiene un conjunto invariante con dinámica simbólica; tiene órbitas periódicas de infinitos períodos distintos e infinitos itinerarios no periódicos. La topología de ese conjunto permite formular una conclusión matemática mucho más fuerte que «la gráfica parece irregular».

Esa conclusión tiene límites. El conjunto con herradura puede no ser atrayente, ocupar una región muy pequeña y ser difícil de observar desde condiciones iniciales genéricas. En un sistema disipativo, las órbitas típicas pueden terminar en un atractor periódico aunque exista un conjunto hiperbólico caótico no atrayente que organice transitorios largos. Tampoco se deduce que toda medida observada tenga una dimensión fractal dada ni que el espectro de Lyapunov calculado sobre una trayectoria finita identifique automáticamente ese conjunto.

El criterio de Melnikov predice una intersección transversal; la transición de un experimento hacia un régimen caótico estable es otra cuestión. Que exista una estructura invariante complicada y que esta tenga peso estadístico para observadores reales son problemas distintos. El segundo se estudia con técnicas numéricas para localizar variedades, cuencas y medidas invariantes.

109.15. Un cálculo reproducible sin confiar en una gráfica

La integral puede comprobarse numéricamente sin una simulación prolongada del Duffing forzado. El código siguiente usa NumPy y SciPy para evaluar la función de Melnikov por cuadratura y compararla con la expresión analítica. La integral sobre toda la recta se aproxima en [−30,30][-30,30]; en este ejemplo, las colas son pequeñas por el decaimiento exponencial de la separatriz. Lo que se comprueba es la etapa algebraica y numérica de la predicción; la persistencia de las variedades perturbadas queda fuera de esta prueba.

import numpy as np
from scipy.integrate import quad

F, gamma, omega = 0.10, 0.10, 1.0

def yh(t):
    return -np.sqrt(2) * np.tanh(t) / np.cosh(t)

def melnikov_integral(t0):
    fun = lambda t: yh(t) * (
        F * np.cos(omega * (t + t0)) - gamma * yh(t)
    )
    return quad(fun, -30, 30, epsabs=1e-11)[0]

def melnikov_formula(t0):
    a = np.sqrt(2) * np.pi * F * omega
    a /= np.cosh(np.pi * omega / 2)
    return a * np.sin(omega * t0) - 4 * gamma / 3

for phase in [0.0, np.pi / 6, np.pi / 2, np.pi]:
    numeric = melnikov_integral(phase)
    analytic = melnikov_formula(phase)
    print(f'{phase:.5f}: {numeric:.9f}, {analytic:.9f}, '
          f'error={abs(numeric-analytic):.2e}')

En t0=0t_0=0 ambas expresiones dan aproximadamente −0,133333-0,133333 y en t0=π/2t_0=\pi/2, cerca de 0,0437320,043732. El cambio de signo entre estas fases indica al menos un cero; la periodicidad permite localizar el segundo. Para aplicar el criterio hay que demostrar además que la derivada no se anula en esos ceros, lo que aquí se hace analíticamente con la derivada del seno.

El código permite detectar un error común, que es olvidar el desplazamiento temporal ωt0\omega t_0 en el argumento de la fuerza. Sin él, la integral toma un valor fijo y se pierde la variable que permite estudiar la intersección de las variedades. Otro error consiste en truncar en un intervalo muy estrecho y atribuir el sesgo de cola a una ruptura real de la fórmula.

109.16. Cómo buscar el fenómeno en una simulación completa

Comprobada la integral, una simulación del Duffing forzado puede integrar el sistema en varias fases del forzamiento y construir el mapa estroboscópico, de período T=2π/ωT=2\pi/\omega. Luego se localiza, por ejemplo con disparo y Newton, el punto fijo hiperbólico del mapa que continúa a la silla. Su derivada se obtiene integrando las ecuaciones variacionales. Los vectores propios de esa derivada dan las direcciones locales desde las que se inicia el cálculo de las ramas estable e inestable.

Para la variedad inestable, se parte de perturbaciones muy pequeñas del punto fijo a lo largo de su vector propio expansivo y se itera el mapa hacia adelante. Para la estable se usa el mapa inverso o un procedimiento de continuación adecuado, ya que una trayectoria integrada hacia adelante desde un estado cualquiera no la reconstruye. Tampoco basta reunir muchos puntos de una única simulación para afirmar que se ha observado un cruce transversal, porque curvas visualmente próximas pueden pertenecer a partes distintas de un conjunto o deberse a una resolución insuficiente.

La comprobación numérica exige determinar la orientación de las tangentes en el punto de encuentro. En una intersección transversal las tangentes son linealmente independientes; una tangencia puede ser muy difícil de distinguir en una imagen sin un cálculo adicional. También hay que examinar la sensibilidad a las tolerancias, la precisión del período de muestreo y el número de iteraciones. La aproximación de primer orden de Melnikov indica en qué región de parámetros conviene hacer estas pruebas.

Por último, un exponente de Lyapunov positivo estimado sobre una trayectoria finita tampoco equivale a la intersección de variedades. Cuando está bien estimado para una medida invariante, informa del crecimiento de las perturbaciones a lo largo de esa dinámica, mientras que el resultado homoclínico describe una estructura geométrica. Ambos resultados pueden apoyarse mutuamente.

109.17. Extensiones y situaciones donde el método no basta

El cálculo clásico presupone una separatriz explícita en un sistema hamiltoniano plano y una perturbación periódica suficientemente pequeña; otros problemas pueden requerir una versión distinta de la función de Melnikov. En sistemas de más grados de libertad aparecen múltiples direcciones estables e inestables y, a veces, hay que calcular funciones vectoriales o considerar variedades de objetos invariantes de dimensión superior. Poder escribir una integral no garantiza que una condición escalar baste para decidir esas intersecciones.

Los sistemas lento-rápidos, estudiados en los capítulos 106 y 107, plantean una dificultad diferente. Cerca de una pérdida de hiperbolicidad normal, las distancias entre variedades pueden ser muy pequeñas y la elección de escala perturbativa es delicada. Un método asintótico que retiene términos algebraicos en ε\varepsilon podría no detectar diferencias exponencialmente pequeñas. El primer paso consiste en identificar el límite perturbativo y el objeto geométrico cuya separación se calcula.

Con perturbaciones aleatorias, la fuerza ya no repite exactamente su fase cada período. Hay herramientas de sistemas dinámicos aleatorios y extensiones estocásticas, pero la función periódica M(t0)M(t_0) no se puede reutilizar mecánicamente. Con datos empíricos, además, suelen desconocerse el campo vectorial y la ubicación de las variedades; una serie temporal irregular que alterna episodios regulares y caóticos no autoriza a inferir una intersección homoclínica.

109.18. Cinco ejercicios resueltos

Ejercicio 1. Energía y tipo de equilibrio. Compruebe que (0,0)(0,0) es una silla del Duffing x˙=y\dot x=y, y˙=x−x3\dot y=x-x^3. La matriz jacobiana en el origen es (0110)\left(\begin{smallmatrix}0&1\\1&0\end{smallmatrix}\right). Sus valores propios, +1+1 y −1-1, tienen signos opuestos, de modo que el origen es hiperbólico de tipo silla. La función de energía vale cero en ese equilibrio; por tanto, su separatriz pertenece al nivel H=0H=0.

Ejercicio 2. Recalcular la disipación. Use yh=−2sech⁡ttanh⁡ty_h=-\sqrt2\operatorname{sech}t\tanh t para evaluar ∫−∞∞yh2dt\int_{-\infty}^{\infty}y_h^2dt. Con u=tanh⁡tu=\tanh t resulta 2∫−11u2du=4/32\int_{-1}^{1}u^2du=4/3. Si el amortiguamiento aparece en la perturbación con coeficiente γ\gamma, su contribución al balance es −4γ/3-4\gamma/3. El signo negativo se conserva sea cual sea el lazo escogido.

Ejercicio 3. Obtener el umbral con frecuencia distinta. Mantenga γ=0,1\gamma=0,1, pero adopte ω=2\omega=2. El factor armónico es aproximadamente 0,766550,76655, luego Fcrit=(0,133333)/0,76655≈0,17394F_{\mathrm{crit}}=(0,133333)/0,76655\approx0,17394. Una fuerza F=0,10F=0,10, suficiente a frecuencia uno, ya no supera el umbral a frecuencia dos. Lo que disminuye es la contribución integrada del forzamiento a lo largo de la separatriz; el sistema no se ha vuelto linealmente más estable.

Ejercicio 4. Distinguir cero simple y contacto. Sea M(ϕ)=Asin⁡ϕ−BM(\phi)=A\sin\phi-B, con A,B>0A,B>0. Si A>BA>B, en los ceros vale ∣sin⁡ϕ∣<1|\sin\phi|<1, por lo que cos⁡ϕ≠0\cos\phi\ne0 y M′≠0M'\ne0. Si A=BA=B, el único cero módulo 2π2\pi está en ϕ=π/2\phi=\pi/2, donde M′=0M'=0. El primer caso satisface el criterio estándar; el segundo requiere términos superiores u otro análisis.

Ejercicio 5. Interpretar una trayectoria irregular. Una simulación del Duffing forzado muestra transitorios irregulares y termina después en una oscilación periódica. ¿Contradice esto la existencia de una intersección homoclínica transversal? No. La geometría transversal puede sostener un conjunto invariante caótico no atrayente que origine transitorios prolongados, mientras otra región contiene el atractor periódico. Para distinguir ambas afirmaciones deben localizarse variedades, cuencas y regiones recurrentes del mapa estroboscópico.

109.19. Una lectura crítica de las aplicaciones

El oscilador de Duffing aparece como aproximación de sistemas elásticos, resonadores, circuitos y estructuras mecánicas con no linealidad cúbica. Esa utilidad no autoriza a traducir automáticamente cualquier transición abrupta de una variable clínica, económica o ecológica a una «intersección de Melnikov». En ámbitos donde no se conoce el mecanismo real, la geometría homoclínica es una hipótesis que exige suficientes variables de estado, una regla dinámica plausible, una excitación temporal caracterizable y evidencia específica de retornos cerca de una estructura hiperbólica.

Podemos imaginar una serie de actividad neuronal que alterna ritmos relativamente regulares y episodios irregulares. Antes de hablar de homoclinicidad, habría que excluir artefactos, ruido no estacionario, forzamientos no medidos y cambios de estado independientes. Incluso si un modelo de oscilador reprodujera la serie, podría hacerlo por mecanismos distintos, porque un ajuste visual no establece unicidad causal. Para sostener una conexión homoclínica convendría encontrar una región de parámetros donde el modelo reconstruido muestre una estructura de silla y validar predicciones adicionales, por ejemplo cambios de los retornos al modificar una frecuencia externa controlada.

La ventaja del método es que convierte una idea geométrica difícil de observar directamente en una desigualdad comprobable. Esa conversión depende, sin embargo, de un modelo reducido, una aproximación asintótica y supuestos geométricos concretos. Conviene distinguir entre el resultado demostrado dentro del modelo, la robustez numérica del cálculo y la evidencia de que la realidad estudiada pertenece a esa clase de modelos.

109.20. Qué prueba el criterio de Melnikov

En el oscilador de Duffing de doble pozo, la separatriz no perturbada da una trayectoria homoclínica conocida. La función de Melnikov calcula al primer orden cómo se separan las variedades estable e inestable cuando añadimos amortiguamiento y forzamiento periódico. Su forma compara la contribución del forzamiento, que oscila con la fase, con la pérdida constante por amortiguamiento.

Con las hipótesis de regularidad, un cero simple de esa función implica intersecciones transversales para perturbaciones suficientemente pequeñas. La herradura de Smale que de ahí se sigue, bajo condiciones apropiadas, no garantiza que una simulación finita muestre un atractor caótico físicamente relevante; la estabilidad local, la dinámica transitoria y la atracción de conjuntos son preguntas adicionales.

Si la perturbación deja de ser pequeña, si la órbita homoclínica no se conoce o si la función de Melnikov no tiene ceros simples, que el método no llegue a una conclusión no demuestra que falten intersecciones. El cálculo orienta la búsqueda; la integración completa y el análisis geométrico la contrastan.

Referencias

  • Guckenheimer, J., & Holmes, P. (1983). Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. Springer.
  • Wiggins, S. (2003). Introduction to Applied Nonlinear Dynamical Systems and Chaos (2.ª ed.). Springer.
  • Ott, E. (2002). Chaos in Dynamical Systems (2.ª ed.). Cambridge University Press.
  • Granados, A. (2016). Notes on the computation of periodic orbits using Newton and Melnikov's method. arXiv:1610.03363. https://arxiv.org/abs/1610.03363
  • Motonaga, S., & Yagasaki, K. (2021). Persistence of periodic and homoclinic orbits, first integrals and commutative vector fields in dynamical systems. arXiv:2108.08026. https://arxiv.org/abs/2108.08026