Exponentes de Lyapunov y sensibilidad

Definición del exponente de Lyapunov

La sensibilidad a condiciones iniciales no significa que dos trayectorias cercanas se separen para siempre a una distancia ilimitada. En una región acotada, la distancia termina por saturarse. El exponente de Lyapunov describe la tasa promedio de crecimiento o contracción de perturbaciones infinitesimales, antes de que esa saturación destruya el régimen lineal.

Para un flujo

$$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x}),$$

consideremos una trayectoria de referencia $\mathbf{x}(t)$ y otra cercana $\mathbf{x}(t)+\delta\mathbf{x}(t)$. Al expandir,

$$\mathbf{f}(\mathbf{x}+\delta\mathbf{x}) =\mathbf{f}(\mathbf{x}) +D\mathbf{f}(\mathbf{x})\delta\mathbf{x} +O(\|\delta\mathbf{x}\|^2).$$

Al restar la ecuación de referencia y retener el término lineal se obtiene la ecuación variacional

$$\dot{\delta\mathbf{x}} =D\mathbf{f}(\mathbf{x}(t))\,\delta\mathbf{x}.$$

El jacobiano cambia a lo largo de la órbita. Por ello los exponentes no se obtienen, en general, calculando autovalores en un punto cualquiera del atractor. Los autovalores de un equilibrio describen su vecindad; el espectro de Lyapunov promedia la evolución tangente a lo largo del régimen visitado.

Para una dirección inicial $\mathbf{v}_0$, una tasa a tiempo finito es

$$\lambda_T(\mathbf{v}_0) =\frac1T\log \frac{\|\delta\mathbf{x}(T)\|}{\|\delta\mathbf{x}(0)\|}.$$

El límite $T\to\infty$, cuando existe bajo las hipótesis adecuadas, produce un exponente característico. En una simulación solo se observan valores finitos, dependientes del transitorio, la duración, la condición inicial y el método numérico. Alligood, Sauer y Yorke separan con cuidado la definición asintótica del algoritmo que la aproxima Alligood, 1996 Kathleen T. Alligood, Tim D. Sauer y James A. Yorke (1996) Chaos: An Introduction to Dynamical Systems Springer Ubicación consultada: §§3.1, 5.1–5.2 y 9.6 .

Exponente de Lyapunov del mapa logístico

Para un mapa unidimensional $x_{n+1}=F(x_n)$, una perturbación pequeña cumple

$$\delta x_{n+1}=F'(x_n)\delta x_n.$$

Después de $N$ iteraciones,

$$\delta x_N= \left[\prod_{n=0}^{N-1}F'(x_n)\right]\delta x_0.$$

Tomar logaritmos convierte el producto en suma:

$$\lambda_N =\frac1N\log\left|\frac{\delta x_N}{\delta x_0}\right| =\frac1N\sum_{n=0}^{N-1}\log|F'(x_n)|.$$

Para el mapa logístico, $F_r'(x)=r(1-2x)$, así que

$$\lambda_N =\frac1N\sum_{n=0}^{N-1}\log|r(1-2x_n)|.$$

Si la órbita converge a un punto fijo $x^*$, la suma converge a $\log|F'(x^*)|$. Un ciclo estable de periodo $p$ produce

$$\lambda=\frac1p\log\left| \prod_{j=0}^{p-1}F'(x_j) \right|<0.$$

En $r=4$, para órbitas típicas respecto de la medida natural, $\lambda=\log2$. La interpretación es una duplicación promedio de una perturbación infinitesimal por iteración. No significa que la distancia observable supere el tamaño de $[0,1]$; al alcanzar ese tamaño, la aproximación lineal deja de describir la separación.

Espectro completo de Lyapunov

En dimensión $d$ existen hasta $d$ exponentes, ordenados como

$$\lambda_1\geq\lambda_2\geq\cdots\geq\lambda_d.$$

Una esfera infinitesimal de perturbaciones se transforma en un elipsoide. Sus ejes principales crecen, en promedio, con tasas $\lambda_i$. Un exponente positivo señala expansión en alguna dirección; uno negativo, contracción. En un flujo autónomo no estacionario suele aparecer un exponente cero asociado a desplazarse a lo largo de la propia trayectoria. “Cero” significa cercano a cero dentro de la convergencia alcanzada, no exactamente cero en cualquier ejecución finita.

La suma del espectro está relacionada con la contracción promedio de volumen:

$$\sum_{i=1}^d\lambda_i =\lim_{T\to\infty}\frac1T \int_0^T\nabla\cdot\mathbf{f}(\mathbf{x}(t))\,dt.$$

En Lorenz clásico, cuya divergencia es $-(\sigma+1+\beta)$, la suma esperada es $-(10+1+8/3)=-41/3\approx-13.667$. Esta igualdad ofrece un control potente: si el espectro calculado tiene una suma muy diferente y no mejora al aumentar $T$ o refinar el paso, se revisan la ecuación variacional, la acumulación QR y la integración.

La presencia de $\lambda_1>0$ en una órbita acotada es evidencia característica de sensibilidad. Aun así, una estimación positiva breve puede originarse en un transitorio o en una región localmente expansiva. Se exige estabilidad del signo y del orden de magnitud al prolongar la medición.

Ortonormalización QR del sistema tangente

Si se propagan varios vectores tangentes sin intervención, todos terminan alineándose con la dirección de mayor crecimiento. Los vectores contractivos se vuelven demasiado pequeños para la precisión disponible y se pierde el resto del espectro. El algoritmo de Benettin resuelve este problema mediante ortogonalizaciones periódicas.

Sea $Q_0$ una matriz ortogonal cuyas columnas son perturbaciones iniciales. Se integra simultáneamente la trayectoria y la matriz variacional durante un intervalo $\tau$:

$$\dot\Phi=J(t)\Phi,\qquad \Phi(0)=Q_0.$$

Al final se factoriza

$$\Phi(\tau)=Q_1R_1,$$

con $Q_1$ ortogonal y $R_1$ triangular superior. Los logaritmos de $|R_{1,ii}|$ registran el crecimiento de cada dirección durante el intervalo. Se reinicia la matriz tangente con $Q_1$, se repite el proceso y se acumula:

$$\lambda_i^{(K)} =\frac{1}{K\tau}\sum_{k=1}^K\log|R_{k,ii}|.$$

La ortogonalización no altera las tasas asintóticas; impide desbordamiento y mantiene separadas las direcciones. Sprott resume este procedimiento y sus controles para flujos simples Sprott, 2010 Julien Clinton Sprott (2010) Elegant Chaos: Algebraically Simple Chaotic Flows World Scientific Ubicación consultada: §§1.9–1.10 .

Un intervalo QR demasiado largo permite que los vectores se alineen y puede desbordar. Uno extremadamente corto aumenta costo y puede amplificar efectos de redondeo. La comparación se hace manteniendo $\tau$ en unidades físicas al cambiar el paso de integración.

Método de dos trayectorias y dinámica tangente

Cuando no se dispone del jacobiano, puede estimarse el exponente máximo con una trayectoria de referencia y otra perturbada. Se parte de una separación $\varepsilon$, se integran ambas durante un intervalo $\tau$, se mide $d_k$, se acumula $\log(d_k/\varepsilon)$ y se devuelve la segunda trayectoria a distancia $\varepsilon$ en la dirección alcanzada. Después de $K$ ciclos,

$$\widehat\lambda_1= \frac{1}{K\tau} \sum_{k=1}^K\log\frac{d_k}{\varepsilon}.$$

La reinyección evita la saturación, pero introduce dos escalas que deben comprobarse. Si $\varepsilon$ es demasiado grande, la diferencia sale del régimen lineal; si es demasiado pequeña, resta números casi iguales y amplifica redondeo. Si $\tau$ es largo, la separación puede saturarse antes de medirla; si es muy corto, el cociente queda dominado por fluctuaciones y error de integración. Se repite el cálculo con varios $\varepsilon$ y $\tau$ y se busca una región estable.

Este método aproxima principalmente la dirección de mayor crecimiento. No recupera por sí solo todo el espectro, porque las perturbaciones terminan alineándose. Para las demás tasas se necesita una base tangente y ortogonalización, o algoritmos equivalentes. Tampoco debe confundirse la distancia euclidiana elegida con una magnitud física universal: un cambio de unidades altera separaciones finitas, aunque los exponentes asintóticos se preserven bajo cambios suaves y acotados de coordenadas.

Los cocientes de cada intervalo definen tasas locales o de tiempo finito. Su histograma muestra que una órbita puede atravesar regiones contractivas y expansivas aunque el promedio final sea positivo. Esas fluctuaciones ayudan a explicar episodios de predictibilidad desigual, pero dependen de $\tau$ y no reemplazan la curva acumulada. El informe presenta ambos objetos con nombres distintos.

Verificación en un sistema lineal

Consideremos el sistema lineal

$$\dot x=2x,\qquad \dot y=-3y.$$

Las soluciones son $x(t)=x_0e^{2t}$ y $y(t)=y_0e^{-3t}$. La matriz fundamental es

$$\Phi(t)= \begin{pmatrix} e^{2t}&0\\ 0&e^{-3t} \end{pmatrix}.$$

Por tanto, los exponentes son $\lambda_1=2$ y $\lambda_2=-3$. Su suma, $-1$, coincide con la divergencia $2-3=-1$. Este sistema sirve como prueba unitaria del código QR: si el algoritmo no recupera ambas tasas, el error no puede atribuirse a caos ni a un transitorio complicado.

Para Lorenz, el jacobiano que se integra junto con el estado es

$$J(x,y,z)= \begin{pmatrix} -\sigma&\sigma&0\\ \rho-z&-1&-x\\ y&x&-\beta \end{pmatrix}.$$

Se valida primero en un parámetro donde el sistema converge a un equilibrio. Allí los exponentes deben ser compatibles con las partes reales de los autovalores del equilibrio alcanzado. Después se pasa al régimen clásico. Esta calibración detecta signos cambiados, orden incorrecto de variables y errores en la reinicialización.

Dimensión de Kaplan–Yorke y otros estimadores

Los exponentes también describen cómo cambia un elemento de volumen en distintas direcciones. Si

$$\lambda_1\ge\lambda_2\ge\cdots\ge\lambda_d$$

y $j$ es el mayor índice para el cual $\lambda_1+\cdots+\lambda_j\ge0$, la dimensión de Kaplan–Yorke se define como

$$D_{KY}=j+ \frac{\lambda_1+\cdots+\lambda_j}{|\lambda_{j+1}|}.$$

La parte entera cuenta direcciones cuya expansión acumulada todavía no es negativa; la fracción interpola cuánto compensa la siguiente dirección contractiva. Ott presenta exponentes y dimensiones fractales como diagnósticos relacionados, pero distingue las definiciones geométricas —por cajas, información o medida puntual— de las relaciones que infieren una dimensión mediante el espectro. Ott, 2002 Edward Ott (2002) Chaos in Dynamical Systems 2.ª ed. · Cambridge University Press Ubicación consultada: cap. 3, §§3.1–3.9, pp. 69–101; cap. 4, §4.4, pp. 137–144

Por ejemplo, con un espectro aproximado $(0.91,0,-14.57)$ se tiene $j=2$ y

$$D_{KY}=2+\frac{0.91}{14.57}\approx2.062.$$

El resultado expresa una predicción espectral: el conjunto ocupa más que una superficie lisa y menos que un volumen. No es una medición directa de dimensión de cajas. Para contrastarlo se estima esa otra dimensión en varias escalas, se identifica una región de pendiente estable y se informa la diferencia. Una coincidencia razonable fortalece la interpretación; una discrepancia obliga a revisar convergencia, escala de ajuste e hipótesis.

Exponentes de Lyapunov de tiempo finito

Sea $\Phi(t_0+T,t_0)$ la matriz que transporta perturbaciones tangentes durante una ventana $T$. Los exponentes de tiempo finito pueden definirse mediante sus valores singulares $\sigma_i$:

$$\lambda_i^{(T)}(t_0) =\frac1T\log\sigma_i\!\left(\Phi(t_0+T,t_0)\right).$$

Mover $t_0$ a lo largo de la misma trayectoria produce una distribución, no un solo número. Una región casi tangente a una dirección estable puede dar una tasa breve negativa y una zona de fuerte estiramiento una positiva, aunque el promedio largo sea positivo. Vallejo y Sanjuán usan precisamente estas distribuciones para separar regímenes y escalas de predictibilidad; también incorporan los errores de modelo y de integración al alcance del pronóstico. Vallejo, 2017 Juan C. Vallejo y Miguel A. F. Sanjuán (2017) Predictability of Chaotic Dynamics: A Finite-time Lyapunov Exponents Approach Springer International Publishing Ubicación consultada: cap. 2, §§2.4–2.8, pp. 33–54; cap. 4, §§4.1–4.2, pp. 91–123 Abrir fuente

La figura semilogarítmica de esta unidad muestra una sola ventana de crecimiento antes de la saturación. Repetir el ajuste en ventanas deslizantes convierte esa pendiente en una serie $\lambda^{(T)}(t_0)$. Su histograma debe acompañarse con $T$: una ventana corta conserva fluctuaciones y una larga las promedia. Cambiar de norma o de coordenadas puede modificar los valores finitos, mientras que los exponentes asintóticos se conservan bajo cambios suaves y acotados en las condiciones habituales.

El efecto sobre el horizonte se puede calcular. Con $\delta_0=10^{-6}$ y una tolerancia $\Delta=10^{-2}$,

$$\log(\Delta/\delta_0)=\log(10^4)\approx9.21.$$

Una ventana con tasa $0.2$ da $T_p\approx46.1$ unidades; otra con tasa $1$ da $T_p\approx9.21$. Un único exponente promedio ocultaría esa diferencia de cinco veces. Estos valores tampoco garantizan que una trayectoria numérica siga a una trayectoria física: redondeo, discretización y error estructural del modelo entran como perturbaciones adicionales, y la propiedad de ensombrecimiento requiere un análisis independiente.

Convergencia y horizonte de predicción

Una tabla mínima incluye paso $h$, tolerancias si el método es adaptativo, transitorio, tiempo de medición $T$, intervalo QR $\tau$, condición inicial, los $d$ exponentes y su suma. Se repite al menos con $T$, $2T$ y con dos pasos. Además de los valores finales se guardan las curvas $\lambda_i(t)$; una meseta visible permite distinguir una cifra estable de un promedio todavía errante.

Tél y Gruiz sitúan el cálculo de parámetros de caos dentro de un conjunto de procedimientos numéricos, no como un botón aislado Tél, 2006 Tamás Tél y Márton Gruiz (2006) Chaotic Dynamics: An Introduction Based on Classical Mechanics Cambridge University Press Ubicación consultada: apéndice A.5 . Los controles más importantes son:

• descartar el transitorio antes de acumular;
• conservar el intervalo QR físico al refinar $h$;
• evitar una perturbación de tamaño finito si se usa dinámica tangente;
• comprobar la suma contra la divergencia media;
• repetir desde otra condición inicial de la misma cuenca;
• informar cuántas cifras permanecen al cambiar $T$ y $h$.

Si $\lambda_1>0$, una incertidumbre inicial $\delta_0$ alcanza una tolerancia $\Delta$ aproximadamente en

$$T_p\approx\frac{1}{\lambda_1} \log\left(\frac{\Delta}{\delta_0}\right).$$

Este horizonte es una estimación local/promedio. La tasa finita fluctúa y los errores de modelo no se reducen a la condición inicial. Por eso se presenta como escala orientativa, con $\delta_0$ y $\Delta$ declarados.

Ejercicios de cálculo y convergencia del espectro

1. Deriva $\lambda=\log|a|$ para $x_{n+1}=ax_n$ y clasifica los casos $|a|<1$, $|a|=1$ y $|a|>1$.
2. Calcula el exponente de un ciclo de periodo dos como promedio de los logaritmos de sus multiplicadores.
3. Verifica el jacobiano de Lorenz componente por componente.
4. Explica qué se pierde si dos vectores tangentes se propagan sin ortogonalización.
5. Construye un criterio para elegir la ventana lineal de una curva de separación semilogarítmica.
6. Con $\lambda_1=0.8$, $\delta_0=10^{-8}$ y $\Delta=10^{-2}$, estima el horizonte de predicción y enumera dos motivos por los que el valor real puede diferir.

El informe del espectro de Lyapunov debe incluir el algoritmo, la curva de convergencia, el control de paso, la suma de exponentes y el horizonte temporal.