Numerical analysis of dynamical systems: unstable periodic orbits, hidden transient chaotic sets, hidden attractors, and finite-time Lyapunov dimension

TL;DR

This paper discusses the challenges of numerically analyzing chaotic dynamical systems, focusing on hidden attractors, transient chaos, and finite-time Lyapunov dimensions using classic models like Lorenz, Rossler, and Vallis.

Contribution

It introduces methods for detecting hidden attractors and accurately computing finite-time Lyapunov exponents and dimensions in low-order chaotic models.

Findings

Identification of hidden chaotic attractors in Lorenz system.

Demonstration of finite-time Lyapunov dimension computation.

Analytical localization of attractors in Vallis system.

Abstract

In this article, on the example of the known low-order dynamical models, namely Lorenz, Rossler and Vallis systems, the difficulties of reliable numerical analysis of chaotic dynamical systems are discussed. For the Lorenz system, the problems of existence of hidden chaotic attractors and hidden transient chaotic sets and their numerical investigation are considered. The problems of the numerical characterization of a chaotic attractor by calculating finite-time time Lyapunov exponents and finite-time Lyapunov dimension along one trajectory are demonstrated using the example of computing unstable periodic orbits in the Rossler system. Using the example of the Vallis system describing the El Nino-Southern Oscillation it is demonstrated an analytical approach for localization of self-excited and hidden attractors, which allows to obtain the exact formulas or estimates of their Lyapunov dimensions.

Figures

Equations

$$\{\begin{array}{l}\dot{x}=-\sigma (x-y),\\ \dot{y}=rx-y-xz,\\ \dot{z}=-bz+xy,\end{array}$$

$$\{\begin{array}{l}\dot{x}=-y-z,\\ \dot{y}=x+0.2y,\\ \dot{z}=0.2-5.7z+xz.\end{array}$$

$$\{\begin{array}{l}\dot{x}=By-C(x+p),\\ \dot{y}=xz-y,\\ \dot{z}=-xy-z+1,\end{array}$$

$$D=\frac{Q^2}{4}+\frac{P^3}{27},Q=\frac{p}{3}\left(\frac{2p^2}{9}+\frac{B}{C}+2\right),P=-\left(\frac{p^2}{3}+\frac{B}{C}-1\right)$$

$$x^3+p{x}^2+\left(1-\frac{B}{C}\right)x+p=0.$$

$$\begin{gathered} x_j=-\frac{p}{3}+{\xi }^{j-1}\sqrt[3]{-\frac{Q}{2}+\sqrt{D}} \\ +{\xi }^{2(j-1)}\sqrt[3]{-\frac{Q}{2}-\sqrt{D}},j=1,2,3,\text{where}\xi =-\frac{1}{2}+\frac{\sqrt{3}}{2}i. \end{gathered}$$

$$\mathcal{B}\!=\!\left\{\!(x,y,z)\!\in\!\mathbb{R}^{3}\,\big{|}\,\tfrac{1}{B}x^{2}\!+\!y^{2}\!+\!(z\!+\!1)^{2}\!\leq\!4\!+\!\tfrac{(Cp)^{2}}{(2C-1)B}\!\right\}\!.$$

$${\dim }_{\mathrm{L}}{\mathcal{A}}_{\mathrm{V}}\le 3-\frac{2(C+2)}{C+1+\sqrt{(C-1{)}^2+\frac{25}{4}B}}.$$

$${\dim }_{\mathrm{L}}{O}_1=2.58834,{\dim }_{\mathrm{L}}{O}_2=2.00368,{\dim }_{\mathrm{L}}{O}_3=\mathrm{2.07175.}$$

$${\dim }_{\mathrm{L}}{O}_2<{\dim }_{\mathrm{L}}{O}_3<{\dim }_{\mathrm{L}}{O}_1\le {\dim }_{\mathrm{L}}\mathcal{A}\le \mathrm{2.65903.}$$

Full text

I Introduction

History of the turbulence phenomena study is associated with the consideration of various models, which include the Navier-Stokes equations, their Galerkin approximations, and the development of the theory of chaos Landau-1944 ; Hopf-1948 ; RuelleT-1971 ; Smale-1967 . Here let us note the significant results by D. Ruelle, F. Takens RuelleT-1971 , and S. Smale Smale-1967 , who proposed a chaotic attractor as a mathematical prototype describing the onset of turbulence, and by O. Ladyzhenskaya, who studied the case when the two-dimensional Navier-Stokes equation generates a dynamical system and proved the finite dimensionality of its attractor Ladyzhenskaya-1982 . The first vivid example of chaotic attractor in a hydrodynamic system was obtained by E. Lorenz Lorenz-1963 . Using the Galerkin method he derived a crude three-dimensional mathematical model for Rayleigh-Bénard convective flow, which has the following form:

[TABLE]

where $r$, $\sigma$, $b$ are positive parameters. For $0<r<1$, there is one globally asymptotically stable equilibrium $S_{0}=(0,0,0)$. For $r>1$, equilibrium $S_{0}$ is a saddle, and a pair of symmetric equilibria $S_{\pm}=(\pm\sqrt{b(r-1)},\pm\sqrt{b(r-1)},r-1)$ appears.

For the parameters $r=28$, $\sigma=10$, $b=8/3$ in system (1) E. Lorenz numerically found a chaotic attractor in the model. In general, for numerical localization of attractor, it is necessary to explore its basin of attraction and choose an initial point in it. If for a particular attractor its basin of attraction is connected with the unstable manifold of unstable equilibrium, then the localization procedure is quite simple. From this perspective, the following classification of attractors is suggested LeonovKV-2011-PLA ; LeonovK-2013-IJBC ; LeonovKM-2015-EPJST ; Kuznetsov-2016 : an attractor is called a self-excited attractor if its basin of attraction intersects an arbitrarily small open neighborhood of an equilibrium; otherwise, it is called a hidden attractor. Numerical localization of hidden attractors is much more challenging and requires the development of special methods. The classical Lorenz attractor is a self-excited one with respect to all equilibria $S_{0}$, $S_{\pm}$, and it is an open question (Kuznetsov-2016, , p. 14) whether for some parameters there exists a hidden Lorenz attractor. This question is related to the ’’chaotic’’ generalization LeonovK-2015-AMC of the second part of Hilbert’s 16th problem on the number and mutual disposition of attractors and repellers in the chaotic multidimensional dynamical systems and, in particular, their dependence on the degree of polynomials in the model; see corresponding discussion, e.g. in SprottJKK-2017 ; ZhangChen-2017 . There are a number of physical dynamical models which possess hidden chaotic attractors, e.g, the Rabinovich system (describes the interaction of plasma waves) KuznetsovLMPS-2018 ; ChenKLM-2017-IJBC , the Glukhovsky-Dolghansky system (describes convective fluid motion in a rotating cavity) LeonovKM-2015-CNSNS ; LeonovKM-2015-EPJST and others DancaKC-2016 ; DudkowskiJKKLP-2016 ; GarashchukSK-2018-HA ; GarashchukSK-2018-RCD .

We note that the Lorenz system (1) with parameters $r=28$, $\sigma=10$, $b=8/3$ is dissipative in the sense of Levinson, and for any initial data (except for equilibria) the trajectory tends to the attractor. Thus, system (1) generates a dynamical system \big{(}\{\varphi^{t}\}_{t\geq 0},(U\subseteq\mathbb{R}^{3},||\cdot||)\big{)}.

II Hidden transient chaotic sets in Lorenz system

In numerical computation of a trajectory over a finite-time interval it is difficult to distinguish a sustained chaos from a transient chaos (a transient chaotic set in the phase space, which can nevertheless persist for a long time) GrebogiOY-1983 ; LaiT-2011 , thus it is reasonable to give a similar classification for transient chaotic sets DancaK-2017-CSF ; ChenKLM-2017-IJBC : a transient chaotic set is called a hidden transient chaotic set if it does not involve and attract trajectories from a small neighborhood of equilibria; otherwise, it is called self-excited. In order to distinguish an attracting chaotic set (attractor) from a transient chaotic set in numerical experiments, one can consider a grid of points in a small neighborhood of the set and check the attraction of corresponding trajectories towards the set. There one can reveal a subset of points for which the trajectories leave the transient set.

For the Lorenz system (1), suppose that $\sigma=10$, $b=\tfrac{8}{3}$ are fixed and $r$ varies. If $r=28$, then all three equilibria $S_{0}$, $S_{\pm}$ are unstable and in the phase space there exists a self-excited attractor with respect to these equilibria. For $r<\sigma(\tfrac{\sigma+b+3}{\sigma-b-1})\approx 24.7368$, the equilibria $S_{\pm}$ become stable, and for $24.06<r<24.7368$, there exists a self-excited attractor with respect to equilibrium $S_{0}$. Near the point $r\approx 24.06$ it is possible to observe a long living transient chaotic set, which is hidden since it’s basin of attraction does not intersect with the small vicinities of equilibrium $S_{0}$. For example, for $r=24$ a hidden transient chaotic set can be visualized111 In this work, we use MATLAB’s standard procedure ode45 with default parameters.  YuanYW-2017-HA from the initial point $(2,2,2)$ (see Fig. 1). In MunmuangsaenS-2018-HA , hidden transient chaotic set was obtained in system (1) with $r=29$, $\sigma=4$, $b=2$.

The time of the transient process in this case depends strongly on the choice of the initial data, which complicates the task of distinguishing an attracting chaotic set (attractor) from a transient chaotic set in numerical experiments. E.g., for system (1) with parameters $r=24$, $\sigma=10$, $b=8/3$ and for initial point $u_{0}=(20,20,20)$ a transient chaotic behavior is observed222 The time of transient chaotic behavior is often estimated approximately by analyzing the sign of the largest Lyapunov exponent. For simplicity, here we approximate the transient behavior time by the time of entering a small ball with the center at the points $S_{\pm}$. on the time interval $[0,\,1.8\!\cdot\!10^{4}]$, for initial point $u_{0}=(-7,8,22)$ — on the time interval $[0,\,7.2\!\cdot\!10^{4}]$, for initial point $u_{0}=(2,2,2)$ — on the time interval $[0,\,2.26\!\cdot\!10^{5}]$, and for initial point $u_{0}=(0,-0.5,0.5)$ a transient chaotic behavior continues over a time interval of more than $[0,\,1\!\cdot\!10^{7}]$. In order to distinguish an attracting chaotic set from a transient chaotic set by computing trajectories on a reasonable time interval, one can consider a grid of points in a small neighborhood of the set and check the attraction of corresponding trajectories towards the set. There one can reveal a subset of points for which the trajectories leave the transient set.

Next, on the example of the Lorenz system (1) we will demonstrate difficulties in the reliable numerical computation of the finite-time Lyapunov exponents and finite-time Lyapunov dimension.

III Finite-time Lyapunov dimension of a transient chaotic set

Consider system (1) with parameters $r=24$, $\sigma=10$, $b=8/3$ and integrate numerically the trajectory with initial data $u_{0}=(20,\,20,\,20)$. We numerically approximate the finite-time Lyapunov exponents and finite-time Lyapunov dimension (see corresponding definitions, e.g. in  Kuznetsov-2016-PLA ; KuznetsovLMPS-2018 ). Integration with $t>T_{1}\approx 1.8\!\cdot\!10^{4}$ leads to the collapse of the ‘‘attractor’’, i.e. the ‘‘attractor’’ turns out to be a transient chaotic set. However, on the time interval $t\in[0,\leavevmode\nobreak\ T_{3}\approx 507883]$ we have $\operatorname{LE}_{1}(t,\,u_{\rm init})>0$ and, thus, may conclude that the behavior is chaotic, and for the time interval $t\in[0,\leavevmode\nobreak\ T_{2}\approx 262954]$ we have ${d}_{\rm L}^{\rm KY}(\{\operatorname{LE}_{i}(t,\,u_{\rm init})\}_{i=1}^{3})>2$. This effect is due to the fact that the finite-time Lyapunov exponents and finite-time Lyapunov dimension are the values averaged over the considered time interval. Since the lifetime of transient chaotic process can be extremely long and taking into account the limitations of reliable integration of chaotic ODEs, even long-time computation of the finite-time Lyapunov exponents and the finite-time Lyapunov dimension does not necessary lead to a relevant approximation of the Lyapunov exponents and the Lyapunov dimension.

On the one hand, computational errors (caused by a finite precision arithmetic and numerical integration of differential equations) and sensitivity to initial data allow one to get a more complete visualization of chaotic attractor (pseudo-attractor) by one pseudo-trajectory computed for a sufficiently large time interval. On the other hand, there arises a question of the reliability of calculating the trajectory itself and its various characteristics, such as finite-time Lyapunov exponents (FTLEs) and finite-time Lyapunov dimension (FTLD), over a long time interval. In KehletL-2013 for the Lorenz system (1) the time interval for reliable computation with $16$ significant digits and error $10^{-4}$ is estimated as $[0,36]$, with error $10^{-8}$ is estimated as $[0,26]$, and reliable computation for a longer time interval, e.g. $[0,10000]$ in LiaoW-2014 , is a challenging task.

For two different long-time pseudo-trajectories $\tilde{u}(t,u^{1}_{0})$ and $\tilde{u}(t,u^{2}_{0})$ visualizing the same attractor, the corresponding FTLEs can be, within the specified accuracy, similar due to averaging over time and similar sets of points $\{\tilde{u}(t,u^{1}_{0})\}_{t\geq 0}$ and $\{\tilde{u}(t,u^{2}_{0})\}_{t\geq 0}$. At the same time, one of the corresponding real trajectories $u(t,u^{1,2}_{0})$ may correspond to an unstable periodic orbit (UPO) which is embedded in the attractor and does not allow one to visualize it. The limitations of the possibilities of numerical integration procedures demonstrates KuznetsovM-2018-arXiv an example of the Rössler system Rossler-1976 :

[TABLE]

For system (2) it is possible to stabilize an unstable periodic orbit (UPO) $u^{\rm upo_{1}}$ with period $\tau\approx 5.8811$ embedded into attractor. The corresponding computations by the standard MATLAB numerical integration procedure with and without application of the Pyragas’ correction control Pyragas-1992 (see also KuznetsovLS-2015-IFAC ) of the largest FTLE, $\operatorname{LE}_{1}(t,u^{\rm upo_{1}}_{0})$, and FTLD $\dim_{\rm L}(t,u^{\rm upo_{1}}_{0})$ along a trajectory with initial data $u^{\rm upo_{1}}_{0}\in u^{\rm upo_{1}}$ over the time interval $[0,500]$ give us the following results. On the initial part of the time interval, one can indicate the coincidence of these values with a sufficiently high accuracy. For the UPO and for the unstabilized trajectory $\operatorname{LE}_{1}(t,u^{\rm upo_{1}}_{0})$ coincide up to the 5th decimal place inclusive on the interval $[0,30.4]$, up to the 4th decimal place inclusive on the interval $[0,53.8]$, up to the 3rd decimal place inclusive on the interval $[0,71.5]$. After $t>71.5$ the difference in values becomes significant and the corresponding graphics diverge in such a way that the part of the graph corresponding to the unstabilized trajectory is lower than the part of the graph corresponding to the UPO (see Fig. 2b). Thus, the application of the Pyragas’ procedure makes it possible to compensate round-off errors and to trace the UPO numerically333 There are well-known cases when the accumulation of errors in the computer representation of real data led to catastrophes (see, e.g. Skeel-1992 ) . Note that other UPOs could be revealed, e.g., by various evolutionary algorithms Zelinka-2015 .

IV Analytical localization of attractors via Lyapunov functions

In order to simplify the numerical search for attractors, one can apply an analytical localization approach related to the dissipativity in the sense of Levinson LeonovKM-2015-EPJST . Also this approach may help to obtain an exact formula or to estimate Lyapunov dimension in the entire phase space. As an example of the effectiveness of such approach, we consider the localization of attractors in the Vallis system describing the El Niño-Southern Oscillation (ENSO) phenomenon of irregular, anomalous, Christmas time warming of the coastal waters of Peru and Ecuador about every 3–6 years that affected weather on a global scale. A low-order model for the ENSO phenomenon was suggested by G.K. Vallis Vallis-1986 and has the following form:

[TABLE]

with parameters $B,C>0$, and real $p$. If $p=0$, then system (3) can be transformed to the Lorenz system (1). The number of equilibria in system (1) depends on the sign of the discriminant

[TABLE]

of the cubic equation

[TABLE]

If $B<C$, then for any real $p$ system (3) has only one equilibrium, Otherwise, if $B>C$, and also p\in(-p^{*},p^{*}),\quad p^{*}=\sqrt{\left(\tfrac{B^{2}}{8C^{2}}-\tfrac{5B}{2C}-1\right)+\tfrac{1}{8}\sqrt{\tfrac{B}{C}\big{(}\tfrac{B}{C}+8\big{)}^{3}}} we get $D<0$ and system (3) has three equilibria $O_{j}=(x_{j},y_{j},z_{j})$, $j=1,2,3$, where $x_{j}$ are the solutions of (4), and $y_{j}=\frac{C(x_{j}+p)}{B}$, $z_{j}=\frac{C(x_{j}+p)}{Bx_{j}}$. One can express $x_{j}$ as follows

[TABLE]

Using the direct Lyapunov method, we can prove the dissipativity of system (1) for $C\geq\frac{1}{2}$ and obtain the following ellipsoidal absorbing set:

[TABLE]

Let $u=(x,\,y,\,z)\in U=\mathbb{R}^{3}$, and the dynamical system $\{\varphi^{t}_{\rm V}\}_{t\geq 0}$, is generated by the Vallis system (3) with positive parameters $B$, $C$ and real parameter $p$, and $\mathcal{A}_{\rm V}\subset\mathbb{R}^{3}$ is a nonempty closed bounded set, which is invariant with respect to the dynamical system $\{\varphi_{\rm V}^{t}\}_{t\geq 0}$. i.e. $\varphi^{t}_{\rm V}(\mathcal{A}_{\rm V})=\mathcal{A}_{\rm V}$ for all $t\geq 0$. Using an effective analytical approach, proposed by Leonov Leonov-1991-Vest ; Kuznetsov-2016-PLA , which is based on a combination of the Douady-Oesterlé approach with the direct Lyapunov method we obtain the upper estimate of the Lyapunov dimension for the global attractor in system (3).

Theorem 1**.**

*For the Vallis dynamical system $\{\varphi^{t}_{\rm V}\}_{t\geq 0}$, generated by system (3) with $B,C>0$, and real $p$ we have the following estimate for the Lyapunov dimension of it global B-attractor *

[TABLE]

Let us compare the obtained estimate (7) with the values of local Lyapunov dimension at the equilibria $O_{1}$, $O_{2}$, $O_{3}$. E.g., for $B=102$, $C=3$, $p=0.83$ we get the following local Lyapunov dimensions

[TABLE]

and the corresponding estimate (7) is in accordance with the values of local Lyapunov dimensions

[TABLE]

Acknowledgement

The work is supported by the Russian Science Foundation project (14-21-00041). This article is based on the lecture at the VII International Conference ’’Problems of Mathematical Physics and Mathematical Modelling’’ (NRNU MEPhI, Moscow, 2018) and is accepted for publication in Journal of Physics: Conference Series.