LECCION 22 | ~30 min

Perturbacion-renormalizacion

Como calcular λ sin algebra lineal

El algoritmo

En la leccion anterior definimos el exponente de Lyapunov y comprendimos su significado profundo: es la tasa exponencial a la que trayectorias cercanas divergen en el espacio de fases. Pero la definicion formal involucra un limite y, en principio, requeriria resolver las ecuaciones variacionales del sistema, es decir, calcular y multiplicar la Jacobiana a lo largo de toda la trayectoria. Para sistemas de alta dimension o ecuaciones complicadas, esto puede ser prohibitivamente costoso.

El algoritmo de perturbacion-renormalizacion, propuesto por Benettin, Galgani, Giorgilli y Strelcyn en 1980, ofrece una alternativa elegante que no requiere ninguna derivada, ninguna Jacobiana, ninguna descomposicion matricial. Solo necesitamos poder integrar las ecuaciones del sistema dos veces en paralelo y medir distancias. Es un metodo puramente geometrico.

El procedimiento consta de seis pasos que se repiten ciclicamente, cada uno con un proposito preciso en la estimacion del exponente de Lyapunov:

Paso 1 — Inicializar

Partimos de dos trayectorias separadas por una distancia $\varepsilon$. La referencia arranca en $\mathbf{x}_0$, la perturbada en $\mathbf{x}_0 + \varepsilon \hat{\mathbf{e}}$, donde $\hat{\mathbf{e}}$ es un vector unitario en cualquier direccion.

Paso 2 — Integrar

Ambas trayectorias evolucionan bajo las mismas ecuaciones durante un intervalo de tiempo $\tau$, el periodo de renormalizacion. Usamos el mismo integrador (RK4) con el mismo paso temporal.

Paso 3 — Medir

Calculamos la nueva separacion $d_i = \|\mathbf{x}'(t_i) - \mathbf{x}(t_i)\|$. Si el sistema es caotico, $d_i \gg \varepsilon$. Si es estable, $d_i \lesssim \varepsilon$.

Paso 4 — Acumular

Sumamos el logaritmo del factor de amplificacion: $\ln(d_i / \varepsilon)$. Cada renormalizacion aporta una muestra de la tasa de divergencia local en esa region del atractor.

Paso 5 — Renormalizar

Devolvemos la perturbada a distancia $\varepsilon$ de la referencia, conservando la direccion: $\mathbf{x}' \leftarrow \mathbf{x} + \varepsilon \cdot (\mathbf{x}' - \mathbf{x}) / d_i$. Esto es el paso critico del algoritmo.

Paso 6 — Repetir

Volvemos al Paso 2 y repetimos $N$ veces. El exponente de Lyapunov se estima como el promedio acumulado de todas las contribuciones logaritmicas dividido por el tiempo total.

Despues de $N$ ciclos de renormalizacion, cada uno de duracion $\tau$, el exponente de Lyapunov maximo se calcula como:

$$\lambda_1 = \frac{1}{N\tau} \sum_{i=1}^{N} \ln \frac{d_i}{\varepsilon}$$

Esta formula es simplemente el promedio temporal de las tasas logaritmicas de expansion, exactamente como prescribe la definicion formal del exponente de Lyapunov. La belleza del metodo es su universalidad: funciona igual para cualquier sistema de ecuaciones diferenciales, mapas iterados, o incluso datos experimentales.

¿Por que renormalizar?

La renormalizacion es el ingrediente critico que distingue este algoritmo de un calculo ingenuo. Sin ella, el metodo fracasa inevitablemente, y la razon es geometrica. Consideremos que ocurre cuando dejamos divergir las dos trayectorias sin intervencion.

En un sistema caotico, la separacion crece exponencialmente: $d(t) \sim \varepsilon \cdot e^{\lambda t}$. Pero el espacio de fases no es infinito. El atractor tiene un diametro finito $D$ (para el Lorenz, $D \approx 40$ unidades). Una vez que la separacion alcanza el orden de $D$, las dos trayectorias estan en regiones completamente distintas del atractor. Su distancia fluctua erraticamente entre 0 y $D$, sin relacion alguna con la tasa local de divergencia exponencial que queremos medir.

El tiempo hasta la saturacion es sorprendentemente corto:

$$t_{\text{sat}} \approx \frac{1}{\lambda} \ln \frac{D}{\varepsilon}$$

Para el Lorenz con $\varepsilon = 10^{-6}$ y $D \approx 40$: $t_{\text{sat}} \approx \frac{1}{0.906} \ln(4 \times 10^7) \approx 19$ unidades de tiempo. Despues de apenas 19 segundos de simulacion, toda la informacion sobre la divergencia local se ha perdido irremediablemente.

La renormalizacion resuelve esto con elegancia. Al devolver periodicamente la perturbacion al regimen lineal (separaciones $\sim \varepsilon \ll D$), aseguramos que cada ciclo mida genuinamente la tasa de expansion exponencial local. Conservar la direccion de la perturbacion es igualmente importante: permite que el vector de separacion se alinee naturalmente con la direccion de maxima expansion, convergiendo al primer vector de Lyapunov.

Eleccion de parametros

La perturbacion ϵ

La eleccion de $\varepsilon$ es un equilibrio entre dos limites. Si $\varepsilon$ es demasiado grande, la perturbacion no es infinitesimal y la hipotesis de linealidad se viola: la separacion ya no crece como una exponencial pura, sino que esta contaminada por efectos no lineales. El resultado sera un $\lambda$ sistematicamente subestimado.

Si $\varepsilon$ es demasiado pequeno, los errores de redondeo del punto flotante de 64 bits (precision de $\sim 10^{-16}$) dominan la diferencia entre las trayectorias. El resultado sera ruidoso e inestable. Para aritmetica de doble precision, el rango optimo es:

$$10^{-8} \lesssim \varepsilon \lesssim 10^{-4}$$

En la practica, $\varepsilon \approx 10^{-6}$ es un valor robusto para la mayoria de los sistemas. Una verificacion esencial: si el $\lambda$ calculado cambia significativamente al variar $\varepsilon$ dentro de este rango, hay un problema con el integrador o con la escala del sistema.

El periodo de renormalizacion τ

El intervalo $\tau$ entre renormalizaciones tiene sus propias restricciones. Si $\tau$ es demasiado corto, el factor $d_i / \varepsilon$ estara muy proximo a 1, y su logaritmo sera un numero pequeno dominado por ruido numerico. Si $\tau$ es demasiado largo, la perturbacion puede saturarse antes de la renormalizacion, perdiendo la informacion que necesitamos.

Una regla practica robusta: $\tau$ debe ser del orden de una unidad de tiempo natural del sistema. Para el Lorenz, una orbita tipica alrededor de un ala del atractor tarda aproximadamente $T \sim 0.7$ a $1.5$ unidades de tiempo, asi que $\tau \approx 0.5$ a $2.0$ funciona bien. Para mapas iterados, $\tau = 1$ (un paso de iteracion) es la eleccion natural.

Convergencia

El exponente de Lyapunov converge lentamente. Cada ciclo de renormalizacion aporta una muestra de $\ln(d_i / \varepsilon)$, pero estas muestras estan correlacionadas temporalmente y tienen varianza considerable. La convergencia sigue la ley del promedio: la fluctuacion tipica del promedio parcial decrece como $1/\sqrt{N}$, lo que significa que necesitamos del orden de $10^4$ renormalizaciones para obtener 2 digitos significativos de $\lambda$, y $10^6$ para 3 digitos.

El promedio movil (running average) es la herramienta clave para monitorear la convergencia. Lo graficamos como funcion del numero de renormalizaciones y buscamos la estabilizacion:

$$\bar{\lambda}_N = \frac{1}{N\tau} \sum_{i=1}^{N} \ln \frac{d_i}{\varepsilon}$$

Al principio, $\bar{\lambda}_N$ oscila ampliamente: las primeras iteraciones estan contaminadas por el transitorio y la orientacion arbitraria de la perturbacion inicial. Gradualmente, el promedio se estabiliza alrededor del valor verdadero. El exponente ha convergido cuando las oscilaciones residuales son menores que la precision deseada.

"La precision del exponente de Lyapunov esta limitada por tres factores independientes: el numero de renormalizaciones N (convergencia estadistica), la precision del integrador numerico (errores de truncacion del RK4), y la magnitud de ϵ (regimen lineal vs. errores de redondeo). Duplicar los digitos significativos requiere cuadruplicar N."

Una verificacion practica indispensable es la independencia respecto a los parametros. Si ejecutamos el algoritmo con $\varepsilon = 10^{-5}$, $\varepsilon = 10^{-6}$, y $\varepsilon = 10^{-7}$, los tres deben converger al mismo $\lambda$ (dentro del error estadistico). Lo mismo con $\tau$: variaciones moderadas no deben alterar el resultado. Si lo alteran, algo esta mal.

El metodo de perturbacion-renormalizacion calcula solo el exponente de Lyapunov maximo $\lambda_1$. Para obtener el espectro completo ($\lambda_1, \lambda_2, \ldots, \lambda_n$), se necesita la generalizacion de Benettin que usa multiples vectores ortogonales y re-ortogonalizacion via Gram-Schmidt o descomposicion QR en cada renormalizacion. Pero para determinar si un sistema es caotico, el exponente maximo es suficiente: si $\lambda_1 > 0$, hay caos.

Ejercicios

Implementacion para el mapa logistico
Implementa el algoritmo de perturbacion-renormalizacion para el mapa logistico $x_{n+1} = r x_n(1-x_n)$ con $r=4$. En este caso $\tau = 1$ iteracion. Calcula $\lambda$ y compara con el valor exacto $\lambda = \ln 2 \approx 0.693$. ¿Cuantas iteraciones necesitas para 3 digitos correctos?
Sensibilidad a ϵ
Usa el LAB de abajo. Mueve el slider de $\varepsilon$ desde $10^{-8}$ hasta $10^{-2}$ y observa como cambia la estimacion de $\lambda$. ¿A partir de que valor la estimacion se degrada? ¿Se degrada por exceso o por defecto? Explica por que, distinguiendo los dos mecanismos de error (no linealidad vs. redondeo).
Efecto de τ
Con $\varepsilon$ fijo en $10^{-6}$, varia $\tau$ en el LAB. ¿La convergencia es mas rapida con $\tau$ grande o pequeno? ¿Hay un valor de $\tau$ a partir del cual el resultado se vuelve incorrecto? Argumenta con la formula $t_{\text{sat}} \approx \frac{1}{\lambda} \ln(D/\varepsilon)$.

LAB: Perturbacion-renormalizacion en vivo

Visualizacion del atractor de Lorenz (proyeccion x-z) con dos trayectorias evolucionando en paralelo. La trayectoria de referencia se muestra en gris, la perturbada en cyan. Cada $\tau$ unidades de tiempo, se produce una renormalizacion: una flecha amarilla indica el desplazamiento de la perturbada de vuelta a distancia $\varepsilon$. El promedio movil de $\lambda_1$ se muestra prominentemente. Ajusta $\varepsilon$ y $\tau$ con los sliders para explorar su efecto.

perturbacion-renorm.js

ϵ = 1e-6 τ = 1.00