Cálculo numérico de los exponentes de Lyapunov sin Jacobiano

Question

Tengo un modelo numérico para el que me gustaría calcular el espectro de Lyapunov. El número de dimensiones es de cientos, y no tengo una expresión analítica para el jacobiano disponible. El sistema es caótico (exponente de Lyapunov mayor positivo) y la dimensión fractal no es enorme (3-5). La integración numérica se realiza con un método de Euler explícito.

¿Cuál sería la forma mejor/simple de calcular el espectro completo de los exponentes de Lyapunov? ¿Alguien sabe de scripts de Matlab que pueda utilizar y/o adaptar? Mi objetivo es calcular la dimensión de Kaplan-Yorke y determinar si el sistema es hipercaótico.

Answer 1

Es casi seguro que no quieres todo el espectro de Lyapunov, sino sólo los dos, tres o quizá cuatro exponentes de Lyapunov más grandes. Estrictamente hablando, dos son suficientes para mostrar el hipercaos, pero es posible que quieras ir a por el primer cero uno como comprobación de cordura. Obtener, por ejemplo, el centésimo exponente de Lyapunov requiere mucho tiempo de cálculo y su resultado será probablemente muy inexacto. Lo más importante es que no se gana nada con ello. Tampoco es necesario para la dimensión de Kaplan-Yorke.

Dicho esto, si no quieres trabajar con el jacobiano, la separación de órbitas es el camino a seguir. Más concretamente, se integran muchas versiones ligeramente diferentes del sistema y se observa cómo evoluciona su separación (los vectores de Lyapunov). Para los exponentes de Lyapunov más allá del primero, tienes que tener cuidado de eliminar regularmente cualquier componente en la dirección de los vectores de Lyapunov anteriores.

Técnicamente, el procedimiento para los dos primeros exponentes de Lyapunov es:

1. Seleccione algunos parámetros $ε$ que representa el tamaño de una pequeña perturbación, y $τ$ que representa el intervalo de reescalado. Para la mayoría de los sistemas, $ε=10^{-12}$ es una buena opción. Utilice un puñado de oscilaciones de su dinámica para  $τ$ pero más adelante se hablará de esta elección.

2. Cree tres instancias de su sistema, y llame para denotar sus estados en el momento $t$ por $y_0(t)$ , $y_1(t)$ y $y_2(t)$ .

3. Establecer la condición inicial $y_0(0)$ a algo en el atractor. Establecer $y_1(0)=y_0(0) +ε·r_1$ y $y_2(0)=y_0(0) +ε·r_2$ , donde $r_1$ y $r_2$ son vectores aleatorios.

4. Integrar para $τ$ unidades de tiempo.

5. Calcula $$v_1(t) = \frac{y_1(t)-y_0(t)}{|y_1(t)-y_0(t)|}; \qquad λ_1(t) = \frac{1}{τ}\log \left( \frac{|y_1(t)-y_0(t)|}{ε} \right).$$

   $v_1$ es su primer vector de Lyapunov. $λ_1$ es su primer exponente local de Lyapunov. Si $|y_1(t)-y_0(t)|$ no es mucho menor que el diámetro de su atractor, su elección de $τ$ es demasiado alto.

6. Calcula $$v_2(t) = \frac{y_2(t)-y_0(t)-\langle y_2, v_1 \rangle v_1}{\left | y_2(t)-y_0(t)-\langle y_2, v_1 \rangle v_1 \right|}; \qquad λ_2(t) = \frac{1}{τ} \log \left(\frac{\left| y_2(t)-y_0(t)-\langle y_2, v_1 \rangle v_1 \right|}{ε} \right).$$

   $v_2$ es su segundo vector de Lyapunov. $λ_2$ es su segundo exponente local de Lyapunov. Obsérvese que esto es análogo al paso 5, salvo que se elimina la proyección sobre $v_1$ .

7. Set $y_1(t) = y_0(t) + εv_1(t)$ y $y_2(t) = y_0(t) + εv_2(t)$ . Este reescalado es asegurarse de que sus separaciones sean pequeñas.

8. Vaya al paso 4.

Después de esto se puede promediar el $λ_1$ y $λ_2$ (excepto los primeros) para obtener los respectivos exponentes de Lyapunov.

También se podría estimar el jacobiano a partir de diferencias finitas, pero eso es sólo sustituir una  $ε$ con otro.