Lyapunov Exponents

MATLAB script for solving augmented IVP with variational equation, calculating Lyapunov exponents and Kaplan—Yorke dimension.

Table of Contents

• Notes on Nonautonomous Systems
• The Variational Equation
• The Lyapunov Exponents
• The Kaplan—Yorke Dimension
• Example

Notes on Nonautonomous Systems

Consider a nonautonomous IVP

$$\dot{\mathbf{z}}\left(t\right) = \mathbf{f}\left(t, \mathbf{z}\right), \quad \left.\mathbf{z}\left(t\right)\right|_{t=t_0}=\mathbf{z}_0,$$

where $\mathbf{z}\left(t\right):\mathbb{R} \mapsto \mathbb{R}^m, ~\mathbf{f}\left(t, \mathbf{z}\right): \mathbb{R}\times\mathbb{R}^m \to \mathbb{R}^m$.

In order to compute the Lyapunov exponents, we first need to convert the nonautonomous system to an autonomous system as follows:

$$z_{m+1} = t, \\\ \dot{\mathbf{y}}\left(t\right) = \mathbf{g}\left(\mathbf{y}\right), \quad \left.\mathbf{y}\left(t\right)\right|_{t=t_0}=\mathbf{y}_0$$

where $\mathbf{y}:\mathbb{R}\mapsto \mathbb{R}^{m+1}, ~\mathbf{g}\left(\mathbf{y}\right): \mathbb{R}^{m+1} \to \mathbb{R}^{m+1}$:

$$\mathbf{y} = \begin{bmatrix} \mathbf{z}\\\ z_{m+1} \end{bmatrix}, \quad \mathbf{g}\left(\mathbf{y}\right) = \begin{bmatrix} \mathbf{f}\left(z_{m+1} , \mathbf{z}\right)\\\ 1 \end{bmatrix}, \quad \mathbf{y}_0 = \begin{bmatrix} \mathbf{z}_0\\\ t_0 \end{bmatrix}.$$

The Variational Equation

Consider an autonomous IVP

$$\dot{\mathbf{z}}\left(t\right) = \mathbf{f}\left(\mathbf{z}\right), \quad \left.\mathbf{z}\left(t\right)\right|_{t=t_0}=\mathbf{z}_0,$$

where $\mathbf{z}\left(t\right):\mathbb{R} \mapsto \mathbb{R}^m, ~\mathbf{f}\left(\mathbf{z}\right): \mathbb{R}^m \to \mathbb{R}^m$. For this system the variational equation has the following form:

$$\dot{\boldsymbol{\Phi}}\left(t\right) = \mathbf{J}_\mathbf{f}\left(\mathbf{z}\right) \boldsymbol{\Phi}\left(t\right), \quad \left.\boldsymbol{\Phi}\left(t\right)\right|_{t=t_0} = \mathbf{I}$$

where $\mathbf{J}_\mathbf{f}\left(\mathbf{z}\right): \mathbb{R} \mapsto \mathbb{R}^{m\times m}$ is Jacobian of the function $\mathbf{f}\left(\mathbf{z}\right)$

$$\mathbf{J}_\mathbf{f}\left(\mathbf{z}\right) = \frac{\text{d}\mathbf{f}}{\text{d}\mathbf{z}} = \begin{bmatrix} \boldsymbol{\nabla}^\top f_1 \\\ \vdots \\\ \boldsymbol{\nabla}^\top f_m \end{bmatrix}= \begin{bmatrix} \dfrac{\partial f_1}{\partial z_1} &\cdots& \dfrac{\partial f_1}{\partial z_m}\\\ \vdots & \ddots & \vdots \\\ \dfrac{\partial f_m}{\partial z_1} &\cdots& \dfrac{\partial f_m}{\partial z_m} \end{bmatrix},$$

and $\boldsymbol{\Phi}\left(t\right): \mathbb{R} \mapsto \mathbb{R}^{m\times m}$ is variational matrix. To find out what happens to the variations, you need to solve the variational equation and the system equation simultaneously. To do this, you work with a new augmented state vector $\mathbf{z}_\ast\left(t\right):\mathbb{R}\mapsto\mathbb{R}^{m + m^2}$ of length $m + m^2$:

$$\dot{\mathbf{z}}_\ast\left(t\right) = \begin{bmatrix} \dot{\mathbf{z}}\left(t\right)\\\ \text{vec}\left(\dot{\boldsymbol{\Phi}}\left(t\right)\right) \end{bmatrix} = \begin{bmatrix} \mathbf{f}\left(\mathbf{z}\right)\\\ \text{vec}\left(\mathbf{J}_\mathbf{f}\left(\mathbf{z}\right)\boldsymbol{\Phi}\left(t\right)\right) \end{bmatrix}$$

The augmented IVP is solved using General Algorithm for the Explicit Runge—Kutta Method.

The Lyapunov Exponents

The calculation of the Lyapunov exponent was based on the QR decomposition method, the application of which can be viewed via the script odeExplicitSolversLyapunovExponents.m.

Syntax

[t, zsol, lyap_exp] = odeExplicitSolversLyapunovExponents(odefun, tspan, tau, incond)
[t, zsol, lyap_exp, dzdt_eval] = odeExplicitSolversLyapunovExponents(..., "Method", method_name, "SafeRegime", flag)

Input Arguments

• odefun: function handle that defines the augmented system of ODEs to be integrated. This function must compute the derivatives for both the original system state and the flattened variational matrix;
• tspan: interval of integration, specified as a two-element vector;
• tau: fixed time discretization step;
• incond: vector of initial conditions. This vector must consist of the initial state vector concatenated with the elements of the initial orthogonal matrix (usually the flattened identity matrix);
• Method (Name-Value Pair): string specifying the explicit Runge-Kutta method to be used. Available options:
  • "RK3", "RK4" (default), "RKB5", "RKN5", "RKB6", "RKB7", "RKCV8", "RKF10", "RKF12", "RKF14";

Output Arguments

• t: column vector of evaluation points used to perform the integration;
• zsol: matrix in which each row corresponds to the full solution vector (state + variations) at the value returned in the corresponding row of t;
• lyap_exp: matrix of Lyapunov exponents in which each row contains the spectrum $[\lambda_1, \dots, \lambda_m]$ calculated at the value returned in the corresponding row of t;
• dzdt_eval: (optional) matrix of the evaluated derivative vectors at each time step.

The Kaplan—Yorke Dimension

Let the Lyapunov exponents be sorted in descending order $\lambda _{1}\geq \lambda _{2}\geq \dots \geq \lambda _{m}$, then

$$D_\text{KY}=k+{\frac {\sum _{i=1}^{k}\lambda_{i}}{|\lambda_{k+1}|}},$$

where for $k$

$$\sum _{i=1}^{k}\lambda _{i}\geq 0, \quad \sum _{i=1}^{k + 1}\lambda _{i}<0.$$

Example

The Rössler Attractor in the ExampleOfUse.mlx was chosen as an example:

$$\begin{gather} \begin{cases} \dot{x} =-y-z,\\\ \dot{y} = x+\alpha y, \\\ \dot{z} = \beta+z\left(x-\varsigma\right), \end{cases} \\\ \begin{bmatrix} \alpha\\\ \beta\\\ \varsigma \end{bmatrix}=\begin{bmatrix} 0.2\\\ 0.2\\\ 5.7 \end{bmatrix}. \end{gather}$$