En C57 introdujimos el exponente de Lyapunov para mapas unidimensionales. Allí una perturbación infinitesimal tenía una sola dirección posible y bastaba una tasa media para describir su crecimiento o contracción. C60 mostró cómo esa tasa puede traducirse, bajo hipótesis explícitas, en un horizonte idealizado de predicción. Desde entonces el problema se ha vuelto geométricamente más rico. Lorenz, Rössler, Hénon y la herradura de Smale viven en espacios de dimensión mayor que uno. Una perturbación ya no tiene una sola orientación, porque puede estirarse en una dirección, contraerse en otra y rotar mientras la trayectoria avanza.

La pregunta natural deja de ser «¿cuál es el exponente?» y pasa a ser cuáles son las tasas asintóticas independientes de expansión y contracción de perturbaciones infinitesimales. El conjunto ordenado de esas tasas es el espectro de Lyapunov.

Este capítulo desarrolla tres ideas. Primero, el espectro se deriva de la dinámica linealizada y no de la distancia entre dos trayectorias finitamente separadas. Segundo, sus exponentes contienen información geométrica sobre longitudes, áreas y volúmenes infinitesimales. Tercero, el mayor exponente puede utilizarse para construir un horizonte predictivo, aunque solo después de especificar la incertidumbre inicial, la tolerancia aceptable y qué otras fuentes de error se están ignorando.

El objetivo es aprender qué afirma el espectro, cómo se calcula y cuáles son los límites de esa afirmación, sin convertirlo en una etiqueta automática de «caos».

68.1. De una perturbación escalar a un vector tangente

Consideremos un sistema autónomo continuo

x˙=f(x),x∈Rd,\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x}), \qquad \mathbf{x}\in\mathbb{R}^d,

con f\mathbf{f} suficientemente regular. Sea x(t)\mathbf{x}(t) una trayectoria. Una perturbación infinitesimal δx(t)\delta\mathbf{x}(t) evoluciona, al primer orden, según la ecuación variacional

δx˙=Df(x(t)) δx,\dot{\delta\mathbf{x}} =D\mathbf{f}(\mathbf{x}(t))\,\delta\mathbf{x},

donde DfD\mathbf{f} es el jacobiano evaluado a lo largo de la trayectoria.

Si Φ(t)\Phi(t) es la matriz fundamental de esta ecuación, definida por

Φ˙(t)=Df(x(t))Φ(t),Φ(0)=I,\dot\Phi(t)=D\mathbf{f}(\mathbf{x}(t))\Phi(t), \qquad \Phi(0)=I,

entonces

δx(t)=Φ(t)δx(0).\delta\mathbf{x}(t)=\Phi(t)\delta\mathbf{x}(0).

La matriz Φ(t)\Phi(t) transporta simultáneamente todas las perturbaciones infinitesimales. En una dimensión era un número. En dd dimensiones es una transformación lineal que puede estirar ciertas direcciones y contraer otras.

Para un mapa diferenciable

xn+1=F(xn),\mathbf{x}_{n+1}=F(\mathbf{x}_n),

la ecuación tangente es

δxn+1=DF(xn)δxn,\delta\mathbf{x}_{n+1}=DF(\mathbf{x}_n)\delta\mathbf{x}_n,

y después de NN pasos

δxN=DF(xN−1)⋯DF(x0)δx0.\delta\mathbf{x}_N= DF(\mathbf{x}_{N-1})\cdots DF(\mathbf{x}_0)\delta\mathbf{x}_0.

El problema de Lyapunov consiste en estudiar las tasas asintóticas asociadas a estos productos de matrices.

68.2. Qué es un espectro de Lyapunov

En condiciones adecuadas, el teorema ergódico multiplicativo de Oseledec garantiza que, para casi toda trayectoria respecto de una medida invariante apropiada, existen tasas asintóticas

λ1≥λ2≥⋯≥λd\lambda_1\ge \lambda_2\ge\cdots\ge\lambda_d

que describen el crecimiento exponencial de perturbaciones en subespacios tangentes distinguidos por la dinámica.

De forma esquemática, para una perturbación con componente no nula en una dirección asociada a λi\lambda_i, su magnitud presenta un crecimiento asintótico del tipo

∥δx(t)∥∼Ceλit\|\delta\mathbf{x}(t)\| \sim C e^{\lambda_i t}

mientras se permanezca en el régimen infinitesimal. Para mapas se reemplaza tt por el número de iteraciones nn.

La notación anterior no debe leerse como si existiera una base fija global formada por vectores propios ordinarios. El jacobiano cambia a lo largo de la órbita. Las direcciones relevantes también cambian. El espectro no surge diagonalizando una única matriz, salvo en sistemas lineales o situaciones especiales.

Los exponentes tienen unidades. En un flujo se miden en inverso de tiempo. En un mapa se expresan por iteración. Comparar directamente un valor 0.90.9 de un flujo con un valor 0.40.4 de un mapa carece de significado si no se especifican las escalas temporales.

68.3. El espectro no es el espectro de autovalores del jacobiano

Esta distinción merece formularse de manera explícita. Los autovalores de

Df(x∗)D\mathbf{f}(\mathbf{x}^*)

en un equilibrio x∗\mathbf{x}^* describen la dinámica lineal local alrededor de ese equilibrio. Los exponentes de Lyapunov de una órbita no estacionaria, en cambio, integran la secuencia completa de jacobianos visitados por esa órbita.

En un sistema caótico, una trayectoria recorre regiones con diferentes tasas locales de expansión y contracción. Puede atravesar zonas donde todos los autovalores instantáneos parecen contractivos y, aun así, acumular una dirección inestable a largo plazo por la combinación de estiramiento, rotación y reinyección.

Por ello, calcular autovalores del jacobiano en una colección de puntos y promediarlos no produce en general el espectro de Lyapunov. Lo correcto es propagar vectores tangentes mediante la ecuación variacional o el producto de jacobianos y controlar numéricamente su crecimiento.

68.4. Longitudes, áreas y volúmenes infinitesimales

El primer exponente controla la tasa asintótica de crecimiento de una longitud infinitesimal genérica en la dirección más expansiva,

ℓ(t)∼ℓ(0)eλ1t.\ell(t)\sim \ell(0)e^{\lambda_1 t}.

El espectro completo contiene más información. El área de un paralelogramo infinitesimal generado por dos vectores independientes crece, bajo hipótesis regulares, a una tasa asociada a λ1+λ2.\lambda_1+\lambda_2.

De manera general, el volumen kk-dimensional máximo generado por kk vectores tangentes tiene una tasa asintótica

λ1+⋯+λk.\lambda_1+\cdots+\lambda_k.

Para k=dk=d aparece el volumen completo del espacio tangente. En un flujo suave,

ddtlog⁡∣det⁡Φ(t)∣=tr⁡Df(x(t))=∇⋅f(x(t)).\frac{d}{dt}\log |\det\Phi(t)| = \operatorname{tr}D\mathbf{f}(\mathbf{x}(t)) = \nabla\cdot\mathbf{f}(\mathbf{x}(t)).

Integrando,

1tlog⁡∣det⁡Φ(t)∣=1t∫0t∇⋅f(x(s)) ds.\frac1t\log|\det\Phi(t)| = \frac1t\int_0^t \nabla\cdot\mathbf{f}(\mathbf{x}(s))\,ds.

En el límite asintótico, cuando las hipótesis necesarias se cumplen,

λ1+⋯+λd=⟨∇⋅f⟩\lambda_1+\cdots+\lambda_d = \left\langle\nabla\cdot\mathbf{f}\right\rangle

respecto de la medida recorrida por la trayectoria.

Para un mapa invertible,

λ1+⋯+λd=⟨log⁡∣det⁡DF(x)∣⟩.\lambda_1+\cdots+\lambda_d = \left\langle \log |\det DF(\mathbf{x})|\right\rangle.

Estas identidades proporcionan controles de consistencia muy útiles para una simulación.

68.5. Qué significa un exponente cero en un flujo

En un flujo autónomo, una perturbación infinitesimal paralela a la dirección de la trayectoria corresponde esencialmente a desplazar el origen temporal. Si x(t)\mathbf{x}(t) es una solución, entonces x(t+ε)\mathbf{x}(t+\varepsilon) representa la misma órbita con una fase temporal ligeramente distinta.

Por ello, sobre una órbita regular no estacionaria y bajo condiciones habituales, aparece un exponente de Lyapunov igual a cero. En un atractor caótico tridimensional típico se espera entonces un patrón

λ1>0,λ2=0,λ3<0.\lambda_1>0, \qquad \lambda_2=0, \qquad \lambda_3<0.

Eso es lo que aparece en el sistema de Lorenz clásico.

La afirmación requiere matiz. Un equilibrio no tiene una dirección temporal tangente no nula del mismo modo, porque f(x∗)=0\mathbf{f}(\mathbf{x}^*)=0. En un punto fijo los exponentes están ligados a la parte real de los autovalores del jacobiano. Tampoco debe exigirse que un cálculo numérico finito entregue exactamente cero, ya que el segundo exponente suele converger lentamente y puede quedar en 10−310^{-3} o 10−210^{-2} durante una integración razonable.

68.6. Por qué no basta con propagar muchos vectores sin renormalizar

Supongamos que queremos calcular los tres exponentes de un sistema tridimensional. Una estrategia ingenua sería iniciar tres vectores independientes y propagarlos mucho tiempo.

El problema es que la dirección asociada a λ1\lambda_1 domina exponencialmente. Salvo casos excepcionales, los tres vectores terminan alineándose con la dirección más expansiva. Numéricamente se pierde la información sobre λ2\lambda_2 y λ3\lambda_3.

La solución estándar consiste en interrumpir periódicamente la propagación y ortogonalizar los vectores. Si QkQ_k contiene una base ortonormal al comienzo de un intervalo y MkM_k es la transformación tangente acumulada durante ese intervalo, formamos

MkQk=Qk+1RkM_k Q_k=Q_{k+1}R_k

mediante una descomposición QR. Los elementos diagonales de RkR_k contienen los factores de expansión que se acumulan antes de la nueva ortonormalización.

Después de muchos intervalos,

λi≈1T∑klog⁡∣(Rk)ii∣.\lambda_i \approx \frac{1}{T} \sum_k \log |(R_k)_{ii}|.

Este procedimiento está asociado a los métodos clásicos de Shimada y Nagashima y de Benettin y colaboradores.

Hay una precaución conceptual. Las columnas de QkQ_k son vectores de Gram–Schmidt producidos por el algoritmo. Aunque permiten recuperar las tasas correctas, no deben identificarse sin más con direcciones intrínsecas covariantes del sistema. El espectro tiene un significado dinámico estable; la base ortogonal utilizada para calcularlo contiene además decisiones numéricas.

68.6.1. Valores singulares, QR y tasas asintóticas

Existe otra forma natural de pensar la deformación tangente. Para un intervalo finito, la matriz de propagación Φ(t)\Phi(t) puede descomponerse mediante valores singulares. Si

Φ(t)=UΣVT,\Phi(t)=U\Sigma V^{\mathsf T},

los elementos diagonales

σ1(t)≥σ2(t)≥⋯≥σd(t)>0\sigma_1(t)\ge\sigma_2(t)\ge\cdots\ge\sigma_d(t)>0

son los factores principales de estiramiento acumulado durante ese intervalo. Las cantidades

1tlog⁡σi(t)\frac{1}{t}\log \sigma_i(t)

son exponentes singulares de tiempo finito. Bajo las hipótesis del teorema multiplicativo, sus límites asintóticos recuperan el espectro de Lyapunov.

Entonces, ¿por qué usar QR en vez de formar Φ(t)\Phi(t) durante un tiempo enorme y calcular una sola descomposición SVD al final? Porque el producto tangente se vuelve rápidamente mal condicionado. Si una dirección crece como e0.9te^{0.9t} y otra decrece como e−14te^{-14t}, después de un intervalo largo las escalas difieren por muchos órdenes de magnitud. La información de las direcciones contractivas puede perderse por redondeo antes de realizar la descomposición final.

La ortonormalización periódica evita esa catástrofe numérica. Cada intervalo conserva factores manejables y los logaritmos se acumulan en lugar de multiplicar números enormes o minúsculos. El resultado no depende de elegir un intervalo QR infinitesimal, aunque un intervalo excesivamente largo puede deteriorar el condicionamiento. Un intervalo demasiado corto aumenta el costo y puede acumular más operaciones de ortogonalización de las necesarias. Por ello el parámetro debe someterse a una prueba de convergencia.

También debe evitarse una confusión frecuente. Los valores singulares de una propagación finita describen el máximo estiramiento sobre ese intervalo específico. Los exponentes de Lyapunov son límites asintóticos ligados a una trayectoria o medida. Un sistema puede presentar amplificación transitoria muy grande aunque todos sus exponentes asintóticos sean negativos. Ese fenómeno aparece, por ejemplo, en sistemas lineales no normales. El espectro responde a la tasa de largo plazo y no a toda posible amplificación de corto plazo.

68.6.2. Normas, coordenadas y qué permanece invariante

La magnitud de una perturbación requiere elegir una norma. En un espacio físico heterogéneo, además, las variables pueden tener unidades distintas. Cambiar de metros a centímetros modifica numéricamente algunos componentes de un vector. Una transformación lineal de coordenadas puede inclinar o escalar las direcciones tangentes.

Los exponentes asintóticos son invariantes bajo cambios de coordenadas suaves y regulares cuyas derivadas y derivadas inversas permanezcan acotadas en la región relevante. La razón es que esas transformaciones solo introducen factores multiplicativos acotados en la norma de una perturbación. Al tomar logaritmos y dividir por un tiempo que tiende a infinito, una contribución acotada desaparece,

1tlog⁡C⟶0.\frac{1}{t}\log C\longrightarrow0.

Esto no significa que todo procedimiento numérico sea inmune a la elección de coordenadas. El condicionamiento de los cálculos, la convergencia de exponentes de tiempo finito y la interpretación de un umbral de error sí pueden cambiar. Si una variable se mide en una escala mil veces mayor que otra, una norma euclídea sin normalización puede hacer que el horizonte operacional quede dominado por esa variable.

Por tanto, conviene separar dos preguntas. El espectro asintótico es una propiedad dinámica robusta frente a transformaciones regulares. El horizonte aplicado a un observable necesita una métrica o tolerancia que tenga significado para la pregunta concreta.

68.7. Caso exacto de control con un flujo lineal diagonal

Antes de aplicar el algoritmo a un atractor caótico conviene probarlo en un sistema cuyo espectro conocemos exactamente,

x˙=0.2x,y˙=−0.5y,z˙=−z.\dot x=0.2x, \qquad \dot y=-0.5y, \qquad \dot z=-z.

La matriz fundamental es

Φ(t)=(e0.2t000e−0.5t000e−t).\Phi(t)= \begin{pmatrix} e^{0.2t}&0&0\\ 0&e^{-0.5t}&0\\ 0&0&e^{-t} \end{pmatrix}.

Por tanto,

λ1=0.2,λ2=−0.5,λ3=−1.\lambda_1=0.2, \qquad \lambda_2=-0.5, \qquad \lambda_3=-1.

La suma es −1.3,-1.3, que coincide con la divergencia constante del campo.

Esta prueba parece trivial, aunque cumple una función metodológica importante. Antes de interpretar un espectro calculado en Lorenz o Hénon, el código debe recuperar resultados conocidos en un problema donde no hay ambigüedad dinámica.

68.8. Espectro del mapa de Hénon

Retomemos el mapa clásico estudiado en C64,

H(x,y)=(1−1.4x2+y,  0.3x).H(x,y)= (1-1.4x^2+y,\;0.3x).

Su jacobiano es

DH(x,y)=(−2.8x10.30),DH(x,y)= \begin{pmatrix} -2.8x & 1\\ 0.3 & 0 \end{pmatrix},

y su determinante es constante,

det⁡DH=−0.3.\det DH=-0.3.

Por ello, si λ1\lambda_1 y λ2\lambda_2 son los exponentes asociados a una órbita típica del atractor,

λ1+λ2=log⁡0.3≈−1.20397.\lambda_1+\lambda_2=\log 0.3 \approx -1.20397.

El cálculo reproducible de este capítulo, después de descartar 5000 iteraciones y acumular 60000 iteraciones adicionales, obtiene aproximadamente

λ1≈0.418,λ2≈−1.622\lambda_1\approx0.418, \qquad \lambda_2\approx-1.622

por iteración.

La suma numérica es −1.2039728043,-1.2039728043, en acuerdo con log⁡0.3\log0.3 dentro del error numérico del cálculo.

Convergencia de los exponentes de Lyapunov del mapa de Hénon con el número de iteraciones acumuladas.
Figura 68.1. Convergencia de los exponentes de Lyapunov del mapa de Hénon con el número de iteraciones acumuladas. Ampliar figura ↗

La interpretación geométrica es coherente con C64. El mapa contrae áreas globalmente, porque la suma de exponentes es negativa, mientras que existe una dirección tangente media expansiva. La coexistencia de expansión direccional y contracción de área permite que un conjunto sea disipativo sin colapsar a una órbita periódica simple.

No hay contradicción entre λ1>0\lambda_1>0 y λ1+λ2<0.\lambda_1+\lambda_2<0.

La primera desigualdad describe expansión lineal en una dirección. La segunda describe contracción de áreas infinitesimales.

68.9. Espectro del sistema de Lorenz

Consideremos ahora

x˙=σ(y−x),\dot x=\sigma(y-x),
y˙=x(ρ−z)−y,\dot y=x(\rho-z)-y,
z˙=xy−βz,\dot z=xy-\beta z,

con los parámetros clásicos

σ=10,ρ=28,β=83.\sigma=10, \qquad \rho=28, \qquad \beta=\frac83.

Como vimos en C62, la divergencia es constante,

∇⋅f=−(σ+1+β)=−413≈−13.6667.\nabla\cdot\mathbf{f} =-(\sigma+1+\beta) =-\frac{41}{3} \approx-13.6667.

Por tanto, el espectro debe satisfacer

λ1+λ2+λ3=−413.\lambda_1+\lambda_2+\lambda_3 =-\frac{41}{3}.

Integrando simultáneamente la trayectoria y la ecuación variacional, con ortogonalización QR cada 0.10.1 unidades de tiempo, el cálculo reproducible obtiene después de 500 unidades de acumulación

λ1≈0.904,λ2≈0.002,λ3≈−14.573.\lambda_1\approx0.904, \qquad \lambda_2\approx0.002, \qquad \lambda_3\approx-14.573.

El segundo valor no debe interpretarse como un exponente genuinamente positivo. Su pequeña desviación respecto de cero es compatible con la convergencia finita del algoritmo. La suma numérica es −13.6666666666,-13.6666666666, prácticamente idéntica a la divergencia analítica.

Convergencia del espectro de Lyapunov del sistema de Lorenz para los parámetros clásicos.
Figura 68.2. Convergencia del espectro de Lyapunov del sistema de Lorenz para los parámetros clásicos. Ampliar figura ↗

El patrón es el esperado para un flujo caótico disipativo tridimensional, con una dirección expansiva, una dirección neutral asociada al avance temporal y una dirección fuertemente contractiva.

68.10. Exponentes de tiempo finito

El espectro asintótico resume una trayectoria larga. No afirma que la expansión sea constante en cada instante.

Podemos definir exponentes de Lyapunov de tiempo finito sobre una ventana de duración τ\tau. Para la dirección dominante, de forma esquemática,

λ1(τ)(t)=1τlog⁡∥δx(t+τ)∥∥δx(t)∥,\lambda_1^{(\tau)}(t) = \frac1\tau \log \frac{\|\delta\mathbf{x}(t+\tau)\|} {\|\delta\mathbf{x}(t)\|},

con la evolución tangente y la renormalización adecuadas.

En Lorenz, los valores calculados sobre bloques de cinco unidades de tiempo varían sustancialmente alrededor del promedio asintótico. En la simulación reproducible, la desviación estándar de esos valores es aproximadamente 0.335.0.335.

Exponentes de Lyapunov de tiempo finito del sistema de Lorenz sobre bloques de cinco unidades de tiempo.
Figura 68.3. Exponentes de Lyapunov de tiempo finito del sistema de Lorenz sobre bloques de cinco unidades de tiempo. Ampliar figura ↗

Esto significa que el sistema atraviesa regiones de expansión más intensa y regiones donde el crecimiento local es débil o incluso contractivo. El valor λ1≈0.904\lambda_1\approx0.904 es una tasa media a largo plazo y no una ley instantánea.

La diferencia importa cuando se habla de predicción. Dos estados con la misma incertidumbre inicial pueden tener horizontes efectivos distintos según la región del atractor donde se encuentren. Un horizonte basado en el exponente asintótico es una escala promedio, no una fecha exacta de pérdida de información para cada condición inicial.

68.11. Del mayor exponente al horizonte de predicción

Supongamos que una incertidumbre inicial genérica de magnitud δ0\delta_0 posee una componente no nula en la dirección más inestable. Mientras la perturbación permanezca suficientemente pequeña,

δ(t)≈δ0eλ1t.\delta(t) \approx \delta_0 e^{\lambda_1 t}.

Definamos una tolerancia Δ\Delta, de modo que a partir de esa separación consideramos que la predicción puntual ya no es suficientemente precisa para la pregunta planteada. Entonces el horizonte idealizado satisface

Δ≈δ0eλ1T∗,\Delta \approx \delta_0 e^{\lambda_1 T_*},

y por tanto

T∗≈1λ1log⁡(Δδ0)T_* \approx \frac{1}{\lambda_1} \log\left(\frac{\Delta}{\delta_0}\right)

para un flujo. En un mapa,

N∗≈1λ1log⁡(Δδ0)N_* \approx \frac{1}{\lambda_1} \log\left(\frac{\Delta}{\delta_0}\right)

se expresa en iteraciones.

Con δ0=10−9\delta_0=10^{-9} y Δ=10−1\Delta=10^{-1}, el espectro calculado da para Lorenz T∗≈20.4T_*\approx20.4 unidades de tiempo. Para Hénon, N∗≈44.1N_*\approx44.1 iteraciones.

Horizonte de predicción en función de la incertidumbre inicial y de la tolerancia aceptada.
Figura 68.4. Horizonte de predicción en función de la incertidumbre inicial y de la tolerancia aceptada. Ampliar figura ↗

Estos números no deben compararse como si fueran la misma unidad física. Sirven para mostrar cómo la tasa de expansión transforma una incertidumbre microscópica en una escala macroscópica.

Una consecuencia importante es logarítmica. Mejorar la precisión inicial por un factor de diez no multiplica por diez el horizonte. Solo añade

log⁡10λ1\frac{\log10}{\lambda_1}

unidades de tiempo o iteraciones. En un sistema con λ1>0\lambda_1>0, la predicción puntual enfrenta rendimientos decrecientes al aumentar la precisión de las condiciones iniciales.

68.12. El horizonte depende de qué incertidumbre se propaga

La fórmula anterior contiene varias hipótesis que deben declararse.

Primero, supone que la incertidumbre inicial tiene una componente en la dirección más expansiva. Si por construcción se restringe a otra dirección, el crecimiento inicial puede estar gobernado por otro exponente. En sistemas no normales, además, puede existir amplificación transitoria que no coincide de forma simple con el espectro asintótico.

Segundo, Δ\Delta depende de la pregunta. Un error de 0.10.1 puede ser enorme para una variable y despreciable para otra. El horizonte de predicción no es una propiedad exclusiva del sistema, porque depende también del criterio de resolución que imponemos.

Tercero, el espectro describe error de estado inicial dentro de un modelo dado. No incluye automáticamente incertidumbre de parámetros, ruido externo, error de observación ni error estructural del modelo. Si esos términos dominan antes que la amplificación caótica, la fórmula basada en λ1\lambda_1 sobreestima la utilidad predictiva real.

Cuarto, la ley exponencial es local en el espacio tangente. Una vez que la perturbación alcanza la escala del atractor, deja de tener sentido extrapolar δ(t)=δ0eλ1t\delta(t)=\delta_0e^{\lambda_1t} indefinidamente. La separación se satura, se pliega y puede incluso disminuir temporalmente.

68.12.1. Orientación de la incertidumbre y horizonte efectivo

La fórmula con λ1\lambda_1 describe una perturbación genérica porque, salvo una alineación especial, cualquier vector pequeño posee alguna componente sobre la dirección de mayor expansión. A medida que la evolución tangente actúa, esa componente termina dominando.

Sin embargo, durante tiempos finitos la orientación inicial puede importar mucho. Imaginemos un sistema tridimensional con espectro λ1>0>λ2>λ3.\lambda_1>0>\lambda_2>\lambda_3.

Una perturbación exactamente situada en el subespacio contractivo puede disminuir inicialmente. Si ese subespacio fuera estrictamente invariante, el crecimiento gobernado por λ1\lambda_1 podría no aparecer para esa perturbación idealizada. En un sistema real, pequeñas componentes debidas a ruido, redondeo o incertidumbre de orientación suelen reintroducir la dirección expansiva.

Esto explica por qué el horizonte no debe presentarse como una constante universal del atractor. Dos distribuciones de error con la misma norma δ0\delta_0 pueden producir distintos tiempos de pérdida de precisión si ocupan direcciones diferentes. Una descripción más completa debería propagar una matriz de covarianza, un conjunto de incertidumbre o un ensamble, en lugar de un radio escalar.

El mismo razonamiento se aplica al observable. Supongamos que el estado completo pierde precisión rápidamente, mientras que la cantidad que interesa es

y=h(x).y=h(\mathbf{x}).

La perturbación observable es aproximadamente

δy≈∇h(x)⋅δx.\delta y\approx \nabla h(\mathbf{x})\cdot\delta\mathbf{x}.

Si ∇h\nabla h es casi ortogonal a la dirección de máxima expansión durante un intervalo, el error en yy puede crecer más lentamente que el error del estado completo. Lo contrario también es posible si el observable amplifica la dirección inestable.

Así, un horizonte de predicción científicamente interpretable necesita especificar al menos tres elementos, a saber, qué incertidumbre inicial se considera, qué variable o función del estado se predice y qué error se acepta. El exponente máximo proporciona una escala dinámica fundamental, aunque no elimina esas decisiones de modelamiento.

68.12.2. Horizonte medio y horizonte local

Los exponentes de tiempo finito permiten refinar la idea anterior. Si sobre un intervalo la tasa dominante efectiva es λ1(τ)(t)\lambda_1^{(\tau)}(t), puede construirse un horizonte local aproximado usando esa tasa en lugar del promedio asintótico. Cuando la trayectoria entra en una región de fuerte expansión, el horizonte se acorta; cuando atraviesa una región contractiva, puede alargarse.

No debe abusarse de esta construcción. Un exponente medido sobre una ventana pasada no es necesariamente una tasa futura conocida. Para utilizarlo predictivamente habría que demostrar que el estado presente permite anticipar la expansión venidera o integrar la dinámica tangente hacia adelante con un modelo disponible. En un modelo conocido esto puede hacerse; en datos observados la situación es más difícil.

La figura 68.3 muestra esa variabilidad. El promedio asintótico estabiliza una tasa global, mientras los bloques de cinco unidades de tiempo conservan fluctuaciones amplias. La coexistencia de ambos niveles evita dos errores simétricos, el de creer que el exponente asintótico describe cada momento y el de concluir que la variabilidad local vuelve inútil toda tasa promedio.

68.13. Predictibilidad de estado y predictibilidad estadística

La pérdida de una predicción puntual no implica que todo pronóstico sea imposible.

Un sistema caótico puede impedir que se prediga con precisión el estado exacto dentro de muchas unidades de tiempo y, aun así, permitir que se estimen distribuciones invariantes, medias temporales, frecuencias de visita o probabilidades de eventos. El espectro de Lyapunov limita una forma concreta de predicción, la de seguir una trayectoria individual con precisión finita.

Esta distinción ya apareció implícitamente en C59 y C60. Un diagrama de bifurcación o una medida estadística puede converger aunque las trayectorias punto a punto de dos simulaciones diverjan. La pregunta correcta no es «¿el sistema es predecible?», sino qué observable se quiere predecir, con qué tolerancia y durante cuánto tiempo.

En aplicaciones empíricas esto es determinante. Un modelo atmosférico puede perder información sobre el estado puntual y conservar información probabilística útil. Un modelo poblacional puede fallar al predecir cada trayectoria individual pero conservar medias de largo plazo. Estas posibilidades no se deducen del espectro por sí solas, porque requieren estudiar los observables y la medida relevante.

68.14. Modelo conocido frente a datos observados

En Lorenz y Hénon conocemos las ecuaciones exactas. Podemos evaluar sus jacobianos, integrar las ecuaciones variacionales y aplicar QR. Ese escenario es mucho más favorable que estimar un espectro desde datos.

Si solo observamos una serie temporal escalar, no conocemos necesariamente el estado completo ni la dinámica tangente. Deben reconstruirse estados, elegir escalas de vecindad y separar crecimiento dinámico de ruido observacional. El método de Wolf y colaboradores fue históricamente importante porque abordó la estimación desde series experimentales, aunque ese problema requiere hipótesis adicionales.

Por ello, un espectro calculado a partir de un modelo y un espectro estimado desde datos no tienen el mismo estatus epistemológico. El primero verifica una propiedad de las ecuaciones implementadas. El segundo es una inferencia dependiente del procedimiento de reconstrucción y de la calidad de los datos.

El currículo retomará este problema en C336, cuando se estudie explícitamente la estimación de exponentes a partir de observaciones.

68.15. Verificación computacional

El paquete reproducible de este capítulo ejecuta tres niveles de control.

Primero, prueba el algoritmo contra el sistema lineal diagonal con espectro exacto (0.2,−0.5,−1)(0.2,-0.5,-1). Esto verifica la mecánica básica de propagación y cálculo de tasas.

Segundo, explota identidades estructurales. En Hénon se comprueba

λ1+λ2=log⁡∣b∣.\lambda_1+\lambda_2=\log|b|.

En Lorenz se comprueba

λ1+λ2+λ3=−(σ+1+β).\lambda_1+\lambda_2+\lambda_3 =-(\sigma+1+\beta).

Tercero, se examina convergencia temporal y variabilidad de tiempo finito. Un único número final sin trayectoria de convergencia puede ocultar que el cálculo aún no ha estabilizado.

Se ejecutaron 36 comprobaciones automáticas y todas fueron satisfactorias. Esto aumenta la confianza en la implementación. No constituye una demostración matemática de que los valores finitos obtenidos sean los exponentes exactos del atractor. Los exponentes asintóticos son límites; cualquier integración termina en tiempo finito.

68.16. Del espectro a la confianza en una trayectoria

C67 mostró que una herradura puede organizar infinitos itinerarios mediante dinámica simbólica. C68 añade una capa métrica, donde distintas direcciones tangentes pueden crecer o contraerse a tasas exponenciales diferentes. En Hénon observamos simultáneamente expansión lineal y contracción de área. En Lorenz aparece el patrón positivo, casi cero y negativo característico de un flujo disipativo caótico tridimensional.

El mayor exponente convierte la sensibilidad en una escala temporal. Sin embargo, el horizonte resultante sigue siendo una estimación basada en un modelo, una tolerancia y una hipótesis sobre la incertidumbre inicial. Después de ese horizonte, la pérdida de coincidencia punto a punto no implica necesariamente que la simulación sea inútil.

La pregunta siguiente es más delicada. Si una trayectoria numérica deja de permanecer cerca de la trayectoria exacta que comenzó en el mismo número representado, ¿puede seguir estando cerca de alguna trayectoria exacta del sistema? Esa pregunta motiva C69, dedicado al sombreado de trayectorias y a los límites de confianza numérica.

Referencias

Benettin, G., Galgani, L., Giorgilli, A., & Strelcyn, J.-M. (1980). Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory. Meccanica, 15(1), 9–20. https://doi.org/10.1007/BF02128236

Benettin, G., Galgani, L., Giorgilli, A., & Strelcyn, J.-M. (1980). Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 2: Numerical application. Meccanica, 15(1), 21–30. https://doi.org/10.1007/BF02128237

Hénon, M. (1976). A two-dimensional mapping with a strange attractor. Communications in Mathematical Physics, 50(1), 69–77. https://doi.org/10.1007/BF01608556

Lorenz, E. N. (1963). Deterministic nonperiodic flow. Journal of the Atmospheric Sciences, 20(2), 130–141. https://doi.org/10.1175/1520-0469(1963)020%3C0130:DNF%3E2.0.CO;2

Oseledec, V. I. (1968). A multiplicative ergodic theorem. Lyapunov characteristic numbers for dynamical systems. Transactions of the Moscow Mathematical Society, 19, 197–231.

Shimada, I., & Nagashima, T. (1979). A numerical approach to ergodic problem of dissipative dynamical systems. Progress of Theoretical Physics, 61(6), 1605–1616. https://doi.org/10.1143/PTP.61.1605

Wolf, A., Swift, J. B., Swinney, H. L., & Vastano, J. A. (1985). Determining Lyapunov exponents from a time series. Physica D, 16(3), 285–317. https://doi.org/10.1016/0167-2789(85)90011-9

Recursos reproducibles

Descargar recursos del capítulo (ZIP, 363 kB).

4 archivos CSV · 4 figuras · 1 programa Python.

El programa se ha ejecutado en una carpeta vacía y genera los archivos de datos y figuras incluidos en el paquete. Se corrigió la carpeta de salida para que el programa guarde los resultados junto al script.

Ver instrucciones y contenido del paquete.