1 Introduction

One of the main tasks of the investigation of dynamical systems is the study of established (limiting) behaviour of the system after transient processes, i.e. the problem of localization and analysis of attractors (limited sets of the system states, which are reached by the system from close initial data after transient processes) [1,2,3]. While trivial attractors (stable equilibrium points) can be easily found analytically or numerically, the search of periodic and chaotic attractors can turn out to be a challenging problem (see, e.g. famous 16th Hilbert problem [4] on the number of coexisting periodic attractors in two-dimensional polynomial systems, which was formulated in 1900 and is still unsolved, and its generalization for multidimensional systems with chaotic attractors [5]). For numerical localization of an attractor, one needs to choose an initial point in the basin of attraction and observe how the trajectory, starting from this initial point, after a transient process visualizes the attractor. Self-excited attractors, even coexisting in the case of multistability [6], can be revealed numerically by the integration of trajectories, started in small neighbourhoods of unstable equilibria, while hidden attractors have the basins of attraction, which are not connected with equilibria and are “hidden somewhere” in the phase space [7,8,9,10]. Remark that in numerical computation of 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 persist for a long time) [11]. The search and visualization of hidden attractors and transient sets in the phase space are challenging tasks [12].

In this paper, we study hidden attractors and transient chaotic sets in the Rabinovich system. It is shown that the methods of numerical continuation and perpetual point are helpful for localization and understanding of hidden attractor in the Rabinovich system.

For the study of chaotic dynamics and dimension of attractors, the concept of the Lyapunov dimension [13] was found useful and became widely spread [14,15,16,17,18]. Since only a finite time can be considered in the numerical analysis of dynamical system, in this paper we develop the concept of the finite-time Lyapunov dimension [19] and approaches to its reliable numerical computation. A new adaptive algorithm for the computation of finite-time Lyapunov dimension and exponents is used for studying the dynamics of the dimension. Various estimates of the finite-time Lyapunov dimension for the Rabinovich hidden attractor in the case of multistability are given.

2 The Rabinovich system: interaction between waves in plasma

Consider a system, studied in 1978 by Rabinovich [20] and Pikovski et al. [21],

$$\begin{aligned} \begin{aligned}&\dot{x} = h y - \nu _1 x - y z, \\&\dot{y} = h x - \nu _2 y + x z, \\&\dot{z} = - z + x y, \end{aligned} \end{aligned}$$
(1)

describing the interaction of three resonantly coupled waves, two of which are parametrically excited. Here, the parameter h is proportional to the pumping amplitude and the parameters \(\nu _{1,2}\) are normalized dumping decrements.

After the linear transformation (see, e.g. [22]):

$$\begin{aligned} \chi : (x,y,z) \rightarrow (\nu _1 \nu _2 h^{-1} y, \nu _1 x, \nu _1 \nu _2 h^{-1} z), \end{aligned}$$
(2)

and time rescaling:

$$\begin{aligned} t \rightarrow \nu _1^{-1} t, \end{aligned}$$
(3)

we obtain a generalized Lorenz system:

$$\begin{aligned} \begin{aligned}&\dot{x} = - \sigma (x - y) - a y z,\\&\dot{y} = r x - y - x z,\\&\dot{z} = -b z + x y, \end{aligned} \end{aligned}$$
(4)

where

$$\begin{aligned} \sigma = \nu _1^{-1} \nu _2,\, b = \nu _1^{-1},\, a = -\nu _2^2 h^{-2},\, r = \nu _1^{-1} \nu _2^{-1}h ^{2}. \end{aligned}$$
(5)

System (4) with \(a = 0\) coincides with the classical Lorenz system [23]. As it is discussed in [22], system (4) can also be used to describe the following physical processes: the convective fluid motion inside rotating ellipsoid, the rotation of rigid body in viscous fluid, the gyrostat dynamics, the convection of horizontal layer of fluid making harmonic oscillations and the model of Kolmogorov’s flow.

Note that since parameters \(\nu _1\), \(\nu _2\), h are positive, the parameters \(\sigma \), b, r are positive and parameter a is negative. From relation (5), we have:

$$\begin{aligned} \sigma = -ar. \end{aligned}$$
(6)

Further, we study system (4) under the assumption (6). If \(r < 1\), then system (4) has a unique equilibrium \({\mathrm {S_0}} = (0,0,0)\), which is globally asymptotically Lyapunov stable (global attractor) [18, 22]. If \(r > 1\), then system (4) has three equilibria: \({\mathrm {S_0}}= (0,0,0)\) and \( {\mathrm {S_{\pm }}} = (\pm x_1, \, \pm y_1, \, z_1), \) where

$$\begin{aligned} x_1 = \frac{\sigma b \sqrt{\xi }}{\sigma b + a \xi }, \quad y_1 = \sqrt{\xi }, \quad z_1 = \frac{\sigma \xi }{\sigma b + a \xi } \end{aligned}$$

and

$$\begin{aligned} \xi = \frac{\sigma b}{2 a^2} \left[ a (r-2) - \sigma + \sqrt{(\sigma - ar)^2 + 4a\sigma } \right] . \end{aligned}$$

The stability of equilibria \(S_{\pm }\) of system (4) depends on the parameters r, a and b. Using the Routh–Hurwitz criterion, we obtain the following

Lemma 1

The equilibria \(S_{\pm }\) of system (4) with parameters (5) are stable if and only if one of the following conditions holds:

  1. (i)

    \(0 \le \, ar + 1 \, {<} \frac{2 r}{r - \sqrt{r(r-1)}}\),

  2. (ii)

    \(ar + 1 < 0\),   \(b {>} b_{cr} = \frac{4 a (r - 1) (ar + 1) \sqrt{r(r-1)} + (ar - 1)^3}{(ar + 1)^2 - 4ar^2}.\)

Proof

The coefficients of the characteristic polynomial \(\chi (x,\,y,\,z) = \lambda ^3 + p_1(x,\,y,\,z) \lambda ^2 + p_2(x,\,y,\,z) \lambda + p_3(x,\,y,\,z)\) of the Jacobian matrix of system (4) at the point \((x,\,y,\,z)\) are the following:

$$\begin{aligned}&p_1(x,\,y,\,z) = b - ar + 1, \\&p_2(x,\,y,\,z) = x^2 + a y^2 - a z^2 - ar (b - r + 1)+ b, \\&p_3(x,\,y,\,z) = - a \big (2 x y z + r x^2 - y^2 + b z^2 - b r (r-1)\big ). \end{aligned}$$

One can check that inequalities \(p_1(x_1,\,y_1,\,z_1) > 0\) and \(p_3(x_1,\,y_1,\,z_1) > 0\) are always valid. If \(ar + 1 \ge 0\), then \(p_1(x_1,\,y_1,\,z_1) p_2(x_1,\,y_1,\,z_1) - p_3(x_1,\,y_1,\,z_1) > 0\), and if condition (i) also holds, then \(p_2(x_1,\,y_1,\,z_1) > 0\).

If \(ar + 1 < 0\), then \(p_2(x_1,\,y_1,\,z_1) > 0\), and if condition (ii) also holds, then \(p_1(x_1,\,y_1,\,z_1) p_2(x_1,\,y_1,\,z_1) - p_3(x_1,\,y_1,\,z_1) > 0\). \(\square \)

3 Attractors and transient chaos

Consider system (4) as an autonomous differential equation of general form:

$$\begin{aligned} \dot{u} = f({u}), \end{aligned}$$
(7)

where \(u=(x,y,z) \in \mathbb {R}^3\), and the continuously differentiable vector function \(f: \mathbb {R}^3 \rightarrow \mathbb {R}^3\) represents the right-hand side of system (4). Define by \({u}(t,{u}_0)\) a solution of (7) such that \({u}(0,{u}_0)={u}_0\). For system (7), a bounded closed invariant set K is

  1. (i)

    a (local) attractor if it is a locally attractive set (i.e. \(\lim _{t \rightarrow +\infty } \mathrm{dist} (K, {u}(t,{u}_0)) = 0\) \(\forall {u_0} \in K(\varepsilon )\), where \(K(\varepsilon )\) is a certain \(\varepsilon \)-neighbourhood of set K),

  2. (ii)

    a global attractor if it is a globally attractive set (i.e. \(\lim _{t \rightarrow +\infty } \mathrm{dist} (K, {u}(t,{u}_0)) = 0\) \(\forall {u_0} \in \mathbb {R}^3\)),

where \(\mathrm{dist}(K, {u}) = \inf _{{v} \in K} ||{v} - {u}||\) is the distance from the point \({u} \in \mathbb {R}^3\) to the set \(K \subset \mathbb {R}^3\) (see, e.g. [9]). Since the whole phase space is a global attractor and any finite union of attractors is again an attractor, it is reasonable to consider only minimal global and local attractors, i.e. the smallest bounded closed invariant set possessing the property (ii) or (i).

Note that system (4) [or (7)] is dissipative in the sense that it possesses a bounded convex absorbing set [9, 22]:

$$\begin{aligned} \mathcal {B}(r,a,b) = \left\{ u \in \mathbb {R}^3 ~|~ V(u) \le \frac{b (\sigma + \delta r)^2}{2 c (a + \delta )} \right\} , \end{aligned}$$
(8)

where \(V(u)=V(x,y,z) = x^2 + \delta y^2 + (a + \delta )\left( z - \frac{\sigma + \delta r}{a + \delta }\right) ^2\), \(\delta \) is an arbitrary positive number such that \(a + \delta > 0\) and \(c = \min (\sigma , 1, \frac{b}{2})\). Thus, the solutions of (4) exist for \(t \in [0,+\,\infty )\) and system (4) possesses a global attractor [9, 24], which contains the set of all equilibria.

Computational errors (caused by a finite precision arithmetic and numerical integration of differential equations) and sensitivity to initial data allow one to get a visualization of chaotic attractor by one pseudo-trajectory computed for a sufficiently large time interval. One needs to choose an initial point in the basin of attraction of the attractor and observe how the trajectory, starting from this initial point, after a transient process visualizes the attractor. Thus, from a computational point of view, it is natural to suggest the following classification of attractors, based on the simplicity of finding the basins of attraction in the phase space.

Definition 1

[7, 9, 10, 25] An attractor is called a self-excited attractor if its basin of attraction intersects with any open neighbourhood of an equilibrium; otherwise it is called a hidden attractor.

For a self-excited attractor, its basin of attraction is connected with an unstable equilibrium and, therefore, self-excited attractors can be found numerically by the standard computational procedure in which after a transient process, a trajectory, starting in a neighbourhood of an unstable equilibrium, is attracted to the state of oscillation and then traces it; then, the computations are being performed for a grid of points in vicinity of the state of oscillation to explore the basin of attraction and improve the visualization of the attractor. Thus, self-excited attractors can be easily visualized (e.g. the classical Lorenz and Hénon attractors are self-excited with respect to all existing equilibria and can be easily visualized by a trajectory from their vicinities).

For a hidden attractor, its basin of attraction is not connected with equilibria and, thus, the search and visualization of hidden attractors in the phase space may be a challenging task. Hidden attractors are attractors in the systems without equilibria (see, e.g. rotating electromechanical systems with Sommerfeld effect (1902) [26, 27]), and in the systems with only one stable equilibrium (see, e.g. counterexamples [7, 28] to Aizerman’s (1949) and Kalman’s (1957) conjectures on the monostability of nonlinear control systems [29, 30]). One of the first related problems is the second part of 16th Hilbert problem [4] on the number and mutual disposition of limit cycles in two-dimensional polynomial systems, where nested limit cycles (a special case of multistability and coexistence of periodic attractors) exhibit hidden periodic attractors (see, e.g. [7, 31, 32]). The classification of attractors as being hidden or self-excited was introduced by Leonov and Kuznetsov in connection with the discovery of the first hidden Chua attractor [25, 33,34,35,36,37,38] and has captured much attention of scientists from around the world (see, e.g. [39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68]).

Since in numerical computation of 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) [11, 69], a similar classification can be introduced for the transient chaotic sets.

Definition 2

[70, 71] A transient chaotic set is called a hidden transient chaotic set if it does not involve and attract trajectories from a small neighbourhood 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 neighbourhood of the set and check the attraction of corresponding trajectories towards the set.

For system (1) with parameters \(\nu _1 = 1\), \(\nu _2 = 4\), and increasing h, it is possible to observe [21] the classical scenario of transition to chaos (via homoclinic and subcritical Andronov-Hopf bifurcations) similar to the scenario in the Lorenz system. For \(4.84 \lessapprox h \lessapprox 13.4\) in system (1), there is a self-excited chaotic attractor (see e.g. Fig. 1), which coexists with two stable equilibria. The same scenario can be obtained for system (4) when parameters \(b > 0\) and \(a < 0\) are fixed and r is increasing. Besides self-excited chaotic attractors, a hidden attractor was found [71, 72] in the system. Note that in [9, 73], system (4) with \(a > 0\) was studied and a hidden attractor was also found numerically.

Fig. 1
figure 1

Multistability in the Rabinovich system (4) with the classical values of parameters \(\nu _1 = 1\), \(\nu _2 = 4\), \(h = 4.92\) from [21]: coexistence of three local attractors—two stable equilibria \(S_{\pm }\) and a chaotic self-excited attractor (self-excited with respect to the unstable zero equilibrium \(S_0\))

Further, we localize a hidden chaotic attractor in system (4) with \(a < 0\) by the numerical continuation method starting from a self-excited chaotic attractor. We change parameters, considered in [72], in such a way that the chaotic set is located not too close to the unstable zero equilibrium to avoid a situation, when numerically integrated trajectory oscillates for a long time and then falls on the unstable manifold of unstable zero equilibrium, leaves the chaotic set and tends to one of the stable equilibria.

3.1 Localization via numerical continuation method

One of the effective methods for numerical localization of hidden attractors in multidimensional dynamical systems is based on the homotopy and numerical continuation method (NCM). It is based on the assumption that the position of the attractor changes continuously with the parameters changing. The idea is to construct a sequence of similar systems such that for the first (starting) system, the initial point for numerical computation of oscillating solution (starting attractor) can be obtained analytically, for example, it is often possible to consider the starting system with a self-excited starting attractor; then, the transformation of this starting attractor in the phase space is tracked numerically while passing from one system to another; the last system corresponds to the system in which a hidden attractor is searched.

For the study of the scenario of transition to chaos, we consider system (7) with \(f(u) = f(u, \lambda )\), where \(\lambda \in \varLambda \subset \mathbb {R}^d\) is a vector of parameters, whose variation in the parameter space \(\varLambda \) determines the scenario. Let \(\lambda _\mathrm{end} \in \varLambda \) define a point corresponding to the system, where a hidden attractor is searched. Choose a point \(\lambda _\mathrm{begin} \in \varLambda \) such that we can analytically or numerically localize a certain nontrivial (oscillating) attractor \(\mathcal {A}^1\) in system (7) with \(\lambda = \lambda _\mathrm{begin}\) (e.g. one can consider a self-excited attractor, defined by a trajectory \({u}^1(t)\) numerically integrated on a sufficiently large time interval \(t \in [0, T]\) with the initial point \({u}^1(0)\) in the vicinity of an unstable equilibrium). Consider a pathFootnote 1 in the parameter space \(\varLambda \) , i.e. a continuous function \(\gamma ~:~ [0,\,1] \rightarrow \varLambda \) for which \(\gamma (0) = \lambda _\mathrm{begin}\) and \(\gamma (1) = \lambda _\mathrm{end}\), and a sequence of points \(\{\lambda ^j\}_{j=1}^k\) on the path, where \(\lambda ^1 = \lambda _\mathrm{begin}\), \(\lambda ^k = \lambda _\mathrm{end}\), such that the distance between \(\lambda ^j\) and \(\lambda ^{j+1}\) is sufficiently small. On each next step of the procedure, the initial point for a trajectory to be integrated is chosen as the last point of the trajectory integrated on the previous step: \({u}^{j+1}(0) = {u}^{j}(T)\). Following this procedure and sequentially increasing j, two alternatives are possible: the points of \(\mathcal {A}^j\) are in the basin of attraction of attractor \(\mathcal {A}^{j+1}\), or while passing from system (7) with \(\lambda = \lambda ^j\) to system (7) with \(\lambda = \lambda ^{j+1}\), a loss of stability bifurcation is observed and attractor \(\mathcal {A}^j\) vanishes. If under changing \(\lambda \) from \(\lambda _\mathrm{begin}\) to \(\lambda _\mathrm{end}\), there is no loss of stability bifurcation of the considered attractors, then a hidden attractor for \(\lambda ^k = \lambda _\mathrm{end}\) (at the end of the procedure) is localized.

3.2 Localization using perpetual points

The equilibrium points of a dynamical system are the ones at which the velocity and acceleration of the system simultaneously become zero. In this section, we show that there are points, termed as perpetual points [74], which may help to visualize hidden attractors.

For system (7), the equilibrium points \({u}_\mathrm{ep}\) are defined by the equation \(\dot{u} = f({u}_\mathrm{ep}) = 0\). Consider a derivative of system (7) with respect to time:

$$\begin{aligned} \ddot{u} = J({u}) \, f({u}) = g({u}), \end{aligned}$$
(9)

where \(J(u) = \left[ \frac{\partial f_i(u)}{\partial u_j} \right] _{i,j=1}^3\) is the Jacobian matrix. Here, g(u) may be termed as an acceleration vector. System (9) shows the variation of acceleration in the phase space.

Similar to the equilibrium points estimation, where we set the velocity vector to zero, we can also get a set of points, where \(\ddot{u} = g({u}_\mathrm{pp}) = 0\) in (9), i.e. the points corresponding to the zero acceleration. At these points, the velocity \(\dot{u}\) may be either zero or nonzero. This set includes the equilibrium points \({u}_\mathrm{ep}\) with zero velocity as well as a subset of points with nonzero velocity. These nonzero velocity points \({u}_\mathrm{pp}\) are termed as perpetual points [12, 74,75,76].

Lemma 2

Perpetual points \(S_\mathrm{pp} = (x_\mathrm{pp}, y_\mathrm{pp}, z_\mathrm{pp})\) of system (4) can be derived from the following system

$$\begin{aligned} \left\{ \begin{aligned}&(x_\mathrm{pp}^2 - a z_\mathrm{pp}^2 + a r^2 - 1 ) \left( -y_\mathrm{pp}^2 + z_\mathrm{pp}^2 + r^2 (a - 1) \right) \\&\qquad = r^2 (a r - 1)^2 - z_\mathrm{pp}^2 (a r - b - 1)^2, \\&2 x_\mathrm{pp} y_\mathrm{pp} z_\mathrm{pp} + r x_\mathrm{pp}^2 - y_\mathrm{pp}^2 + b z_\mathrm{pp}^2 - b r ( r - 1 ) = 0, \\&z_\mathrm{pp} = \frac{r x_\mathrm{pp}^2 + (a r - b - 1) x_\mathrm{pp} y_\mathrm{pp} - a r y_\mathrm{pp}^2}{x_\mathrm{pp}^2 + a y_\mathrm{pp}^2 - b^2}. \end{aligned} \right. \end{aligned}$$
(10)

The reason why perpetual points may lead to hidden states [perpetual point method (PPM)] is still discussed (see, e.g. [77, 78]).

4 Hidden attractor in the Rabinovich system

Next, we apply the NCM for localization of a hidden attractor in the Rabinovich system (4) and check whether the attractor can be also localized using PPM.

In this experiment, we fix parameter r and, using condition (ii) of Lemma 1, define parameters \(a = - \frac{1}{r} - \varepsilon _1\) and \(b = b_{cr} - \varepsilon _2\). For \(r = 100\), \(\varepsilon _1 = 10^{-3}\), \(\varepsilon _2 = 10^{-2}\), we obtain \(a = a_0 \equiv -\,1.1\,\times \, 10^{-2}\), \(b = b_0 \equiv 6.7454\,\times \, 10^{-2}\) and take \(P_0(a_0, \, b_0)\) as the initial point of line segment on the plane \((a, \,b)\). The eigenvalues of the Jacobian matrix at the equilibria \(S_0\), \(S_\pm \) of system (4) for these parameters are the following:

$$\begin{aligned}&S_0 \, : \quad 9.4382, \quad -0.0675, \quad -\,11.5382, \\&S_{\pm } \, : \quad 0.0037 \pm 3.6756i, \quad -\,2.1749. \end{aligned}$$

Consider on the plane \((a,\, b)\) a line segment, intersecting a boundary of stability domain of the equilibria \(S_{\pm }\) with the final point \(P_2(a_2, \, b_2)\), where \(a_2 = a_0 + 1.035\,\times \,10^{-3} = -\,9.965\,\times \,10^{-3}\), \(b_2 = b_0 + \varepsilon _2 = 7.7454\,\times \,10^{-2}\), i.e. the equilibrium \(S_0\) remains saddle and the equilibria \(S_\pm \) become stable focus-nodes

$$\begin{aligned} S_0 \,&: \quad 8.9842, \quad -\,0.0775, \quad -\,10.9807, \\ S_{\pm } \,&: \quad -\,0.0401 \pm 3.9152i, \quad -\,1.9937. \end{aligned}$$
Fig. 2
figure 2

Path \(P_0(a_0, \, b_0) \rightarrow P_1(a_1, \, b_1) \rightarrow P_2(a_2, \, b_2)\) in parameters plane \((a,\,b)\) for the localization of hidden attractor in system (4) with \(r = 100\). Here, \(a_0 = -\,1.1\,\times \,10^{-2}\), \(b_0 = 6.7454 \,\times \,10^{-2}\), \(a_1 = -\,1.049\,\times \,10^{-2}\), \(b_1 = 7.2454\,\times \,10^{-2}\), \(a_2 = -\,9.965\,\times \,10^{-3}\), \(b_2 = 7.7454\,\times \,10^{-3}\); \((\bullet )\) \(P_0(a_0, \, b_0)\): self-excited attractor with respect to \(S_{0}\), \(S_{\pm }\); \((\bullet )\) \(P_1(a_1, \, b_1)\): self-excited attractor with respect to \(S_{0}\); \((\bullet )\) \(P_2(a_2, \, b_2)\): hidden attractor. Stability domain is defined using Lemma 1

The initial point \(P_0(a_0, \, b_0)\) corresponds to the parameters for which in system (4), there exists a self-excited attractor. Then for the considered line segment, a sufficiently small partition step is chosen, and at each iteration step of the procedure, an attractor in the phase space of system (4) is computed. The last computed point at each step is used as the initial point for the computation at the next step. In this experiment, we use NCM with 3 steps on the path \(P_0(a_0, \, b_0) \rightarrow P_1 (a_1, \, b_1)\rightarrow P_2 (a_2, \, b_2)\), with \(a_1 = \frac{1}{2}(a_0 + a_2)\), \(b_1 = \frac{1}{2}(b_0 + b_2)\) (see Fig. 2). At the first step, we have self-excited attractor with respect to unstable equilibria \(S_0\) and \(S_\pm \); at the second step, the equilibria \(S_\pm \) become stable but the attractor remains self-excited with respect to equilibrium \(S_0\); at the third step, it is possible to visualize a hidden attractor of system (4) (see Fig. 3).

Fig. 3
figure 3

Localization, by NCM, of a hidden attractor in system (4) with \(r = 100\), \(a = -\,9.965\,\times \, 10^{-3}\), \(b = 7.7454\,\times \,10^{-2}\). Trajectories \(u^i(t) = (x^i(t), y^i(t), z^i(t)\) (blue) are defined on the time interval \([0, T= 10^3]\), and initial point (yellow) on \((i+1)\)-th iteration is defined as \(u_0^{i+1} := u_T^{i}\) (violet arrows), where \(u_T^{i} = u^{i}(T)\) is a final point (yellow). (Color figure online)

Using Lemma 2 for parameters \(r = 100\), \(a = -\,9.965\,\times \,10^{-3}\), \(b = 7.7454\,\times \,10^{-2}\), we obtain one perpetual point \(S_\mathrm{pp} = (-\,0.2385, 49.1403, -\,101.4613)\), which allows one to localize a hidden attractor (see Fig. 4). Hence, here both NCM and PPM allow one to find this hidden attractor.

Around equilibrium \(S_0\), we choose a small spherical vicinity of radius \(\delta \) (in our experiments, we check \(\delta \in [0.1,0.5]\)) and take N random initial points on it (in our experiment, \(N = 4000\)). Using MATLAB, we integrate system (4) with these initial points in order to explore the obtained trajectories. We repeat this procedure several times in order to get different initial points for trajectories on the sphere. We get the following results: all the obtained trajectories attract to either the stable equilibrium \(S_+\) or the equilibrium \(S_-\) and do not tend to the attractor. This gives us a reason to classify the chaotic attractor, obtained in system (4), as the hidden one.

Remark that there exist hidden chaotic sets in the Rabinovich system, which cannot be localized by PPM. For example, for parameters \(r = 6.8\), \(a = -0.5\), \(b \in [0.99,\, 1]\) [72], the hidden attractor obtained by NCM is not localizable via PPM (see Fig. 5).

5 Computation of the Lyapunov dimension

5.1 Finite-time and limit values of the Lyapunov dimension and Lyapunov exponents

The Lyapunov exponents [2] and Lyapunov dimension [13] are widely used for the study of attractors (see, e.g. [14,15,16,17,18]). Nowadays, various approaches to the definition of Lyapunov dimension are used. Since in numerical experiments we can consider only finite time, in this paper we develop the concept of finite-time Lyapunov dimension [19] and an approach to its reliable numerical computation. This approach is inspirited by the works of Douady and Oesterlé [79], Hunt [80], and Rabinovich et al. [81].

Fig. 4
figure 4

Localization of hidden attractor in system (4) with \(r = 100\), \(a = -9.965\,\times \,10^{-3}\), \(b = 7.7454\,\times \,10^{-2}\) from the perpetual point \(S_\mathrm{pp} = (-\,0.2385, 49.1403, -\,101.4613)\)

Fig. 5
figure 5

Hidden attractor in system (4) with \(r = 6.8\), \(a = -\,0.5\), \(b = 0.99\), which can be localized via NCM (initial point \({u}_0 = (-\,0.1629,\, -\,0.2154, \, 1.9553)\)), but cannot be localized from the perpetual point \(S_\mathrm{pp} = (0.7431, -\,12.6109, -\,7.4424)\) (yellow). (Color figure online)

For a fixed \(t\ge 0\) let us consider the map \(\varphi ^t: \mathbb {R}^3 \rightarrow \mathbb {R}^3\) defined by the shift operator along the solutions of system (7): \(\varphi ^t(u_0) = u(t,u_0)\), \(u_0 \in \mathbb {R}^3\). Since system (7) possesses an absorbing set [see (8)], the existence and uniqueness of solutions of system (7) for \(t \in [0,+\infty )\) take place and, therefore, the system generates a dynamical system \(\{\varphi ^t\}_{t\ge 0}\). Let a nonempty closed bounded set \(K \subset \mathbb {R}^3\) be invariant with respect to the dynamical system \(\{\varphi ^t\}_{t\ge 0}\), i.e. \(\varphi ^t(K) = K\) for all \(t \ge 0\) (e.g. K is an attractor). Further, we use compact notations for the finite-time local Lyapunov dimension: \(\dim _\mathrm{L}(t,u) = \dim _\mathrm{L}(\varphi ^t, u)\), the finite-time Lyapunov dimension: \(\dim _\mathrm{L}(t,K) = \dim _\mathrm{L}(\varphi ^t, K)\) and for the Lyapunov dimension (or the Lyapunov dimension of dynamical system \(\{\varphi ^t\}_{t\ge 0}\) on invariant set K): \(\dim _\mathrm{L}K = \dim _\mathrm{L}(\{\varphi ^t\}_{t \ge 0}, K)\).

Consider linearization of system (7) along the solution \(\varphi ^t(u)\):

$$\begin{aligned} \begin{aligned}&\dot{v} = J(\varphi ^t(u))v, \quad J(u) = Df(u), \end{aligned} \end{aligned}$$
(11)

where J(u) is the \(3\times 3\) Jacobian matrix, the elements of which are continuous functions of u, and suppose that \(\det J(u) \ne 0 \quad \forall u \in \mathbb {R}^3\). Consider a fundamental matrix of solutions of linear system (11), \(D\varphi ^t(u)\), such that \(D\varphi ^0(u) = I\), where I is a unit \(3\times 3\) matrix. Denote by \(\sigma _i(t,u) = \sigma _i(D\varphi ^t(u))\), \(i = 1,2,3\), the singular values of \(D\varphi ^t(u)\) (i.e. \(\sigma _i (t, u) > 0\) and \({\sigma _i (t, u)}^2\) are the eigenvalues of the symmetric matrix \(D\varphi ^t(u)^*D\varphi ^t(u)\) with respect to their algebraic multiplicity),Footnote 2 ordered so that \(\sigma _1(t,u) \ge \sigma _2(t,u) \ge \sigma _3(t,u) > 0\) for any \(u \in \mathbb {R}^3\) and \(t > 0\).

A singular value function of order \(d \in [0,3]\) is defined as

$$\begin{aligned} \begin{aligned}&\omega _d(D\varphi ^t(u)) = \sigma _1(t, u)\cdots \sigma _{\lfloor d \rfloor }(t, u) \sigma _{\lfloor d \rfloor +1}(t, u)^{d-\lfloor d \rfloor }, \\&\omega _0(D\varphi ^t(u)) = 1, \omega _3(D\varphi ^t(u)) = \sigma _1(t, u)\sigma _2(t, u)\sigma _3(t, u), \end{aligned} \end{aligned}$$

where \({\lfloor d \rfloor }\) is the largest integer less or equal to d. For a certain moment of time \(t \ge 0\), the finite-time local Lyapunov dimension at the point u is defined as [19]

$$\begin{aligned} \dim _\mathrm{L}(t,u) = \max \{d \in [0,3]: \omega _{d}(D\varphi ^t(u)) \ge 1 \} \end{aligned}$$
(12)

and the finite-time Lyapunov dimension of K is defined as

$$\begin{aligned} \dim _\mathrm{L}(t, K) = \sup \limits _{u \in K} \dim _\mathrm{L}(t,u). \end{aligned}$$
(13)

The Douady–Oesterlé theorem [79] implies that for any fixed \(t \ge 0\), the Lyapunov dimension of the map \(\varphi ^t\) with respect to a closed bounded invariant set K, defined by (13), is an upper estimate of the Hausdorff dimension of the set K: \(\dim _\mathrm{H}K \le \dim _\mathrm{L}(t, K)\).

For the estimation of the Hausdorff dimension of invariant closed bounded set K, one can use the map \(\varphi ^t\) with any time t (e.g. \(t=0\) leads to the trivial estimate \(\dim _\mathrm{H}K \le 3\)) and, thus, the best estimation is \( \dim _\mathrm{H}{K} \le \inf _{t\ge 0}\dim _\mathrm{L} (t, K). \) The following property

$$\begin{aligned} \inf _{t\ge 0}\sup \limits _{u \in K} \dim _\mathrm{L}(t,u) = \liminf _{t \rightarrow +\infty }\sup \limits _{u \in K} \dim _\mathrm{L}(t,u) \end{aligned}$$
(14)

allows one to introduce the Lyapunov dimension of K as [19]

$$\begin{aligned} \dim _\mathrm{L} K = \liminf _{t \rightarrow +\infty }\sup \limits _{u \in K} \dim _\mathrm{L}(t,u) \end{aligned}$$
(15)

and get an upper estimation of the Hausdorff dimension:

$$\begin{aligned} \dim _\mathrm{H}{K} \le \dim _\mathrm{L} K. \end{aligned}$$

Recall that a set K with noninteger Hausdorff dimension is referred to as a fractal set [15] and, when such set K is an attractor, it is called a strange attractor [82, 83].

Consider a set of finite-time Lyapunov exponentsFootnote 3 at the point u:

$$\begin{aligned} {{\mathrm{LE}}}_i(t,u) = \frac{1}{t}\ln \sigma _i(t,u), \ t > 0, \quad i=1,2,3. \end{aligned}$$
(16)

Here, the set \(\{{{\mathrm{LE}}}_i(t,u)\}_{i=1}^3\) is ordered by decreasing (i.e. \({{\mathrm{LE}}}_1(t,u) \ge {{\mathrm{LE}}}_2(t,u) \ge {{\mathrm{LE}}}_3(t,u)\) for all \(t>0\)). Then for \(j(t, u) = \lfloor \dim _\mathrm{L} (t, u) \rfloor <3\) and \(s(t, u) = \dim _\mathrm{L} (t, u) - \lfloor \dim _\mathrm{L} (t, u) \rfloor \), we have \(0 = \frac{1}{t} \ln (\omega _{j(t,u)+s(t,u)}(D\varphi ^t(u))) = \sum _{i=1}^{j(t,u)}{{\mathrm{LE}}}_i(t,u) + s(t, u){{\mathrm{LE}}}_{j(t,u)+1}(t,u) \) and \(j(t,u) = \max \{m: \sum _{i=1}^{m}{{\mathrm{LE}}}_i(t,u) \ge 0\}\). Thus, we get an analog of the Kaplan–Yorke formula [13] with respect to the set of finite-time Lyapunov exponents \(\{{{\mathrm{LE}}}_i(t,u)\}_{i=1}^3\) [19]:

$$\begin{aligned} d_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LE}}}_i(t,u)\}_{i=1}^3)= j(t,u)+\tfrac{{{\mathrm{LE}}}_1(t,u) + \cdots + {{\mathrm{LE}}}_{j(t,u)}(t,u) }{|{{\mathrm{LE}}}_{j(t,u)+1}(t,u)|},\nonumber \\ \end{aligned}$$
(17)

which gives the finite-time local Lyapunov dimension:

$$\begin{aligned} \dim _\mathrm{L}(t, u) = d_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LE}}}_i(t,u)\}_{i=1}^3). \end{aligned}$$

Thus, in the above approach, the use of Kaplan–Yorke formula (17) with the finite-time Lyapunov exponents \(\{{{\mathrm{LE}}}_i(t,u)\}_{i=1}^3\) is rigorously justified by the Douady–Oesterlé theorem.

Note that the Lyapunov dimension is invariant under Lipschitz diffeomorphisms [19, 86], i.e. if the dynamical system \(\{\varphi ^t\}_{t\ge 0}\) and closed bounded invariant set K under a smooth change of coordinates \(w = \chi (u)\) are transformed into the dynamical system \(\{\varphi _\chi ^t\}_{t\ge 0}\) and closed bounded invariant set \(\chi (K)\), respectively, then \(\dim _\mathrm{L}(\{\varphi ^t\}_{t\ge 0},K) = \dim _\mathrm{L}(\{\varphi _\chi ^t\}_{t\ge 0},\chi (K)). \) Also the Lyapunov dimension is invariant under positive time rescaling \(t \rightarrow at, a>0\).

5.2 Adaptive algorithm for the computation of the finite-time Lyapunov dimension and exponents

Applying the statistical physics approach and assuming the ergodicity (see, e.g. [13, 87,88,89]), the Lyapunov dimension of attractor \(\dim _\mathrm{L} K\) is often estimated by the local Lyapunov dimension \(\dim _\mathrm{L} (t, u_0)\), corresponding to a “typical” trajectory, which belongs to the attractor: \(\{u(t,u_0), t \ge 0 \},\ u_0 \in K\), and its limit value \(\lim _{t\rightarrow +\infty }\dim _\mathrm{L} (t, u_0)\). However, from a practical point of view, the rigorous verification of ergodicity is a challenging task [84, 90] and hardly it can effectively be done in a general case (see, e.g. discussions in [91, 92, p. 118], [93, 94, p. 9], [95, p. 19], and the works [96, 97] on the Perron effects of the largest Lyapunov exponent sign reversals). Examples of the rigorous use of the ergodic theory for the computation of the Lyapunov exponents and dimension can be found, for example, in [87, 98,99,100]. In the work by Yorke et al. [88, p. 190], the exact limit values of finite-time Lyapunov exponents \(\{\lim \nolimits _{t\rightarrow +\infty }{{\mathrm{LE}}}_i(t,u_0)\}_{i}^3 = \{{{\mathrm{LE}}}_i(u_0)\}\), if they exist and are the same for all \(u_0 \in \mathbb {R}^3\), i.e.

$$\begin{aligned} \begin{aligned}&\lim \limits _{t\rightarrow +\infty }{{\mathrm{LE}}}_i(t,u_0) = {{\mathrm{LE}}}_i(u_0) \equiv {{\mathrm{LE}}}_i \quad \forall u_0 \in \mathbb {R}^3, \ i=1,2,3, \end{aligned} \end{aligned}$$

are called the absolute ones, and it is noted that the absolute Lyapunov exponents rarely existFootnote 4 (in this case, one also has \( \dim _\mathrm{L} K = d_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LE}}}_i\}_{1}^3) = j + \tfrac{{{\mathrm{LE}}}_1 + \cdots + {{\mathrm{LE}}}_{j}}{|{{\mathrm{LE}}}_{j+1}|} \)).

Note also that even if a numerical approximation (visualization) \(\widetilde{K}\) of the attractor K is obtained, it is not straightforward how to get a point on the attractor itself: \(u_0 \in K\). Thus, an easy way to get reliable estimation of the Lyapunov dimension of attractor K is to localize the attractor \(K \subset K^{\varepsilon }\), to consider a grid of points \(K^{\varepsilon }_\mathrm{grid}\) on \(K^{\varepsilon }\), and to find the maximum of the corresponding finite-time local Lyapunov dimensions for a certain time \(t=T\): \(\max \limits _{u \in K^{\varepsilon }_\mathrm{grid}} \dim _\mathrm{L}(\varphi ^T,u) =\max \limits _{u \in K^{\varepsilon }_\mathrm{grid}} j(T,u) + \frac{{{\mathrm{LE}}}_1(T,u) + \cdots + {{\mathrm{LE}}}_{j(T,u)}(T,u)}{|{{\mathrm{LE}}}_{j(T,u)+1}(T,u)|}\).

Concerning the time T, remark that while the time series obtained from a physical experiment are assumed to be reliable on the whole considered time interval, the time series produced by the integration of mathematical dynamical model (7) can be reliable on a limited time interval onlyFootnote 5 due to computational errors (caused by finite precision arithmetic and numerical integration of ODE). Thus, in general, the closeness of the real trajectory \(u(t,u_0)\) and the corresponding pseudo-trajectory \(\tilde{u}(t,u_0)\) calculated numerically can be guaranteed on a limited short time interval only. For the numerical visualization of a chaotic attractor, the computation of a pseudo-trajectory on a longer time interval often allows one to obtain a more complete visualization of the attractor (pseudo-attractor) due to computational errors and sensitivity to initial data. However, 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 finite-time LEs can be, within the considered error, similar due to averaging over time [see (16)] and similar sets of points \(\{\tilde{u}(t,u^1_0)\}_{t \ge 0}\) and \(\{\tilde{u}(t,u^2_0)\}_{t \ge 0}\). At the same time, the corresponding real trajectories \(u(t,u^{1,2}_0)\) may have different LEs (e.g. \(u_0\) may correspond to an unstable periodic trajectory \(u(t,u_0)\), which is embedded in the attractor and does not allow one to visualize it). Here, one may recall the question [105, p. 98] (known as Eden conjecture) whether the supremum of the local Lyapunov dimensions is achieved on a stationary point or an unstable periodic orbit embedded in the strange attractor. Thus, the numerical computation of trajectory for a longer time may not lead to a more precise approximation of the Lyapunov exponents and dimension (see, also example (24) below). Also any long-time computation may be insufficient to reveal the limit values of LEs if the trajectory belongs to a transient chaotic set, which can be (almost) indistinguishable numerically from sustained chaos (see Figs. 9a, 10a). Note that there is no rigorous justification of the choice of t and it is known that unexpected jumps of \(\dim _\mathrm{L}(t,K)\) can occur (see, e.g. Fig. 7). Thus, for a finite time interval \(t \in [0,T)\), it is reasonable to compute \(\inf _{t \in [0,T)}\dim _\mathrm{L}(t,K)\) instead of \(\dim _\mathrm{L}(T,K)\), but, at the same time, for any T the value \(\dim _\mathrm{L}(T,K)\) gives an upper estimate of \(\dim _\mathrm{H}K\).

Finally, in the numerical experiments, based on the finite-time Lyapunov dimension definition (15) [19] and the Douady–Oesterlé theorem [79], we have

$$\begin{aligned} \begin{aligned}&\dim _\mathrm{H}K \le \dim _\mathrm{L}K \approx \inf _{t \in [0,T]} \max _{u \in K^{\varepsilon }_\mathrm{grid}} \dim _\mathrm{L}(t, u) \\&\quad \quad \quad \quad =\inf _{t \in [0,T]}\max \limits _{u \in K^{\varepsilon }_\mathrm{grid}} \left( j(t,u) + \tfrac{{{\mathrm{LE}}}_1(t,u) + \cdots + {{\mathrm{LE}}}_{j(t,u)}(t,u)}{|{{\mathrm{LE}}}_{j(t,u)+1}(t,u)|} \right) \\&\quad \quad \quad \quad \le \max _{u \in K^{\varepsilon }_\mathrm{grid}} \dim _\mathrm{L}(T, u) \approx \dim _\mathrm{L}(T,K). \end{aligned} \end{aligned}$$
(18)

Additionally, we can consider a set of random points in \(K^{\varepsilon }\). If the maximums of the finite-time local Lyapunov dimensions for random points and grid points differ more than by \(\delta \), then we decrease the distance between grid points. This may help to improve the reliability of the result and at the same time to ensure its repeatability.

Nowadays, there are several widely used approaches to numerical computation of the Lyapunov exponents; thus, it is important to state clearly how the LEs are being computed [92, p. 121]. Next, we demonstrate some differences in the approaches. For a certain \(u_0 \in \mathbb {R}^3\), to compute the finite-time Lyapunov exponents according to (16), one has to find the fundamental matrix \(\varPhi (T,u_0) = D\varphi ^T(u_0)\) of (11) from the following variational equation:

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{u}(t,u_0)=f(u(t,u_0)), \ \ u(0,u_0) = u_0, \\ \dot{\varPhi }(t,u_0)=J(u(t,u_0)) \, \varPhi (t,u_0), \ \varPhi (0,u_0)=I, \end{array}\right. } t \in [0,T], \end{aligned}$$
(19)

and its singular value decomposition (SVD)Footnote 6:

$$\begin{aligned} \varPhi (T,u_0) \overset{\text {SVD}}{=}U(T,u_0){\Sigma }(T,u_0)\mathrm{V}(T,u_0)^{*}, \end{aligned}$$

where \(U(T,u_0)^*U(T,u_0) \equiv I \equiv \mathrm{V}(T,u_0)^*\mathrm{V}(T,u_0)\), \({\Sigma }(T,u_0)=\text {diag}\{\sigma _1(T,u_0),\sigma _2(T,u_0),\sigma _3(T,u_0)\}\) is a diagonal matrix composed by the singular values of \(\varPhi (T,u_0)\), and compute the finite-time Lyapunov exponents \(\{{{\mathrm{LE}}}_i(T,u_0)\}_1^3\) from \({\Sigma }(T,\,u_0)\) as in (16). Further, we also need the QR decompositionFootnote 7

$$\begin{aligned} \varPhi (T,u_0) \overset{\text {QR}}{=} Q(T,u_0)R(T,u_0), \end{aligned}$$

where \(R(T,u_0)\) is an upper triangular matrix with nonnegative diagonal elements \(\{R[i,i]=R[i,i](T,u_0)\}_1^3\) and \(Q(T,u_0)^*Q(T,u_0) \equiv I\).

To avoid the exponential growth of values in the computation, the time interval has to be represented as a union of sufficiently small intervals, for example, \((0, T] = (0,\tau ]\cup (\tau ,2\tau ]\cdots \cup ((k-1)\tau ,k\tau =T]\). Then, using the cocycle property, the fundamental matrix can be represented as

$$\begin{aligned} \varPhi (k \tau ,u_0) = \varPhi (\tau ,u_{k-1}) \, \dots \, \varPhi (\tau ,u_{1}) \, \varPhi (\tau ,u_{0}). \end{aligned}$$
(20)

Here, if \(\varPhi (m \tau ,u_0)\) and \(u_m=u(m \tau ,u_0)\) are known, then \(\varPhi ((m+1) \tau ,u_0) = \varPhi (\tau ,u_{m}) \varPhi (m \tau ,u_0)\), where \(\varPhi (\tau ,u_{m})\) is the solution of initial value problem (19) with \(u(0) = u_m\) on the time interval \([0, \tau ]\).

By sequential QR decomposition of the product of matrices in (20), we get

$$\begin{aligned} \begin{aligned}&\varPhi (k \tau ,u_0) = \varPhi (\tau ,u_{k-1}) \ldots \varPhi (\tau ,u_{1}) \, \boxed {\varPhi (\tau ,u_{0})} \\&\quad = \varPhi (\tau ,u_{k-1}) \ldots \boxed {\varPhi (\tau ,u_{1}) \, Q^0_1} \, R^0_1\\&\quad = \cdots \overset{\text {QR}}{=} \, \overbrace{Q^0_k}^{Q} \,\overbrace{R^0_k \ldots R^0_1}^{R}. \end{aligned} \end{aligned}$$

Then, the matrix with singular values

$$\begin{aligned} \varSigma (k\tau , u_0) = U^*(k\tau ,u_0) \,\varPhi (k \tau ,u_0) \, V(k\tau ,u_0) \end{aligned}$$

in the SVD can be approximated by sequential QR decomposition of the product of matrices:

$$\begin{aligned} \begin{aligned}&\varSigma ^0 = \varPhi (k \tau ,u_0)^* \, Q^0_k = (R^0_1)^* \ldots (R^0_k)^* \overset{\text {QR}}{=} \, Q^1_k \, R^1_k \ldots R^1_1, \\&\varSigma ^1 = (Q^0_k)^* \, \varPhi (k \tau ,u_0) \, Q^1_k = (R^1_1)^* \ldots (R^1_k)^* \overset{\text {QR}}{=} Q^2_k R^2_k \ldots R^2_1, \\&\cdots \end{aligned} \end{aligned}$$

where

$$\begin{aligned} \begin{aligned}&\varSigma ^{p} = (R^{p}_1)^* .. (R^{p}_k)^* =\;\left( \begin{matrix} \sigma _1^{p} &{}\quad 0 &{}\quad 0 \\ \cdot &{}\quad \sigma _2^{p} &{}\quad 0 \\ \cdot &{}\quad \cdot &{}\quad \sigma _3^{p} \end{matrix} \right) \end{aligned} \end{aligned}$$

is a lower triangular matrix and [107, 108]

$$\begin{aligned} \sigma _i^{p} = R^{p}_1[i,i] \ldots R^{p}_k[i,i] \underset{p \rightarrow \infty }{\longrightarrow } \sigma _i(k\tau ,u_0), \quad i=1,2,3. \end{aligned}$$

Thus, the finite-time Lyapunov exponents can be approximated as

$$\begin{aligned} {{\mathrm{LE}}}_i(T, u_0)\approx {{\mathrm{LE}}}_i^{p}(k\tau , u_0)=\tfrac{1}{t}\ln \sigma _i^{p}= \tfrac{1}{k\tau }\sum _{l = 1}^{k} \ln R_{l}^p[i,i]. \end{aligned}$$
(21)

The MATLAB implementation of the above method for the computation of finite-time Lyapunov exponents with the fixed number of iterations p can be found, for example, in [9]. For large k, the convergence can be very rapid: for example, for the Lorenz system with the classical parameters (\(r = 28\), \(\sigma = 10\), \(b = 8/3\), \(a = 0\)), \(k = 1000\) and \(\tau = 1\), the number of approximations \(p=1\) is taken in [108, p. 44]. For a more precise approximation of the finite-time Lyapunov exponents, we can adaptively choose \(p=p(l)\), \(l=1,\ldots ,k\) so as to obtain a uniform estimate:

$$\begin{aligned} \begin{aligned} \max _{i=1,2,3}\left| {{\mathrm{LE}}}_i^{p(l)-1}(l\tau , u_0)-{{\mathrm{LE}}}_i^{p(l)}(l\tau , u_0)\right| < \delta .\quad l=1,\ldots ,k. \end{aligned} \end{aligned}$$
(22)

Remark that there is another widely used definition of the “Lyapunov exponents” via the exponential growth rates of norms of the fundamental matrix columns \(\big (v_1(t,u_0),v_2(t,u_0),v_3(t,u_0)\big )=\varPhi (t,u_0)\): the finite-time Lyapunov characteristic exponents \(\{{{\mathrm{LCE}}}_i(t,u_0)\}_{1}^3\) are the set \(\{\frac{1}{t}\ln ||v^i(t,u_0)||\}_{1}^3\) ordered by decreasing.Footnote 8 Relying on ergodicity [84], Benettin et al. [109] (Benettin’s algorithm; see also Wolf et al. [110]) approximate the LCEs by (21) with \(p=0\):

$$\begin{aligned} {{\mathrm{LCE}}}_i(k\tau , u_0)\approx {{\mathrm{LE}}}_i^{0}(k\tau , u_0) = \frac{1}{k\tau }\sum _{l = 1}^{k} \ln R_{l}^0[i,i]. \end{aligned}$$
(23)

The LCEs may differ from LEs, thus, the corresponding Kaplan–Yorke formulas with respect to LEs and LCEs: \(\dim _\mathrm{L}(t, u_0) =d_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LE}}}_i(t,u_0)\}_1^3)\) and \(d_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LCE}}}_i(t,u_0)\}_1^3)\),Footnote 9 may not coincide. The following artificial analytical example demonstrates the possible difference between LEs and LCEs. The matrix [19, 86]

$$\begin{aligned} R(t)=\left( \begin{array}{cc} 1 &{} g(t)-g^{-1}(t)\\ 0 &{} 1 \\ \end{array} \right) , \ g(t)=\exp (\tfrac{t}{10}) \end{aligned}$$

has the following exact limit values

$$\begin{aligned} \begin{aligned}&{{\mathrm{LCE}}}_1 =\lim \limits _{t \rightarrow +\infty } t^{-1} \ln g(t) = 0.1, \quad {{\mathrm{LCE}}}_2 = 0, \\&{{\mathrm{LE}}}_{1,2} = \lim \limits _{t \rightarrow +\infty }t^{-1} \ln g^{\pm 1}(t) = \pm 0.1, \end{aligned} \end{aligned}$$

where \({{\mathrm{LCE}}}_2 \ne {{\mathrm{LE}}}_2\). For the finite-time values we have

$$\begin{aligned} \begin{aligned}&{{\mathrm{LCE}}}_1(t)=\tfrac{1}{t} \ln \big ((g(t)-\frac{1}{g(t)})^2+1\big )^{\tfrac{1}{2}} \in (0, 0.1],\\&\quad {{\mathrm{LCE}}}_2(t)\equiv 0, \\&{{\mathrm{LE}}}_{1,2}(t) \equiv {{\mathrm{LE}}}_{1,2} = \pm 0.1. \end{aligned} \end{aligned}$$

Approximations by the above algorithm with \(k=1\) are given in Table 1.

Table 1 Approximation of the finite-time Lyapunov exponents

Remark that here the approximation of LCEs and LEs by Benettin’s algorithm, i.e. by (23), becomes worse with increasing time:

$$\begin{aligned} \begin{aligned}&{{\mathrm{LCE}}}_1(t) \underset{t \rightarrow +0}{\longrightarrow } 0 \equiv {{\mathrm{LE}}}_1^{0}(t) = \frac{1}{t} \ln 1 \equiv 0, \\&{{\mathrm{LCE}}}_1(t) \underset{t \rightarrow +\infty }{\longrightarrow } 0.1 \equiv {{\mathrm{LE}}}_1(t) \ne {{\mathrm{LE}}}_1^{0}(t) = \frac{1}{t}\ln 1\equiv 0. \end{aligned} \end{aligned}$$
(24)

Thus, although the notions of LCEs and LEs often do not differ (see, e.g. Eckmann and Ruelle [15, p. 620, p. 650], Wolf et al. [110, p. 286, pp. 290–291] and Abarbanel et al. [17, p. 1363, p. 1364]) relying on the Oseledec ergodic theorem [84], in general case, the computations of LCEs by (21) and LEs by (23) may give unreliable results. See also numerical examples below and counterexamples in [111, p. 289], [97, p. 1083].

5.3 Estimation of the Lyapunov dimension without integration of the system and the exact Lyapunov dimension

While analytical computation of the Lyapunov exponents and Lyapunov dimension is impossible in a general case, they can be estimated by the eigenvalues of the symmetrized Jacobian matrix [79, 112]. Let \(\{\lambda _i(u)\}_{i=1}^{3}\) be the eigenvalues of the symmetrized Jacobian matrix \(\frac{1}{2} \left( J(u) + J(u)^{*}\right) \), ordered so that \(\lambda _1(u) \ge \lambda _2(u) \ge \lambda _3(u)\). The Kaplan–Yorke formula with respect to the ordered set of eigenvalues of the symmetrized Jacobian matrix [19] gives an upper estimation of the Lyapunov dimension:

$$\begin{aligned} \dim _\mathrm{L}K \le \sup _{u \in K}d_\mathrm{L}^\mathrm{KY}\big (\{\lambda _{i}(u)\}_{i=1}^3\big ). \end{aligned}$$

In a general case, one cannot get the same values of \(\{\lambda _{i}(u)\}_{i=1}^3\) at different points u; thus, the maximum of \(d_\mathrm{L}^\mathrm{KY}\big (\{\lambda _{i}(u)\}_{i}^3\big )\) on K has to be computed. To avoid numerical localization of the set K, we can consider an analytical localization, e.g. by the absorbing set \(\mathcal {B} \supset K\). Then for the corresponding grid of points \(\mathcal {B}_\mathrm{grid}\), we expect in numerical experiments the following:

$$\begin{aligned}&\dim _\mathrm{H}K \le \dim _\mathrm{L}K \le \sup _{u \in K} d_\mathrm{L}^\mathrm{KY}\big (\{\lambda _{i}(u)\}_{1}^3\big ) \nonumber \\&\quad \le \sup _{u \in \mathcal {B}} d_\mathrm{L}^\mathrm{KY}\big (\{\lambda _{i}(u)\}_{1}^3\big ) \approx \max _{u \in \mathcal {B}_\mathrm{grid}} \Big (j(u) \nonumber \\&\quad \quad + \tfrac{\lambda _{1}(u) + \cdot \cdot + \lambda _{j(u)}(u)}{|\lambda _{j(u)+1}(u)|}\Big ). \end{aligned}$$
(25)

If the Jacobian matrix \(J(u_\mathrm{eq})\) at one of the equilibria has simple real eigenvalues: \(\{\lambda _i(u_\mathrm{eq})\}_{i=1}^3\), \(\lambda _{i}(u_\mathrm{eq}) \ge \lambda _{i+1}(u_\mathrm{eq})\), then [19, 86] the invariance of the Lyapunov dimension with respect to linear change of variables implies

$$\begin{aligned} \dim _\mathrm{L}u_\mathrm{eq} = d_\mathrm{L}^\mathrm{KY}(\{\lambda _i(u_\mathrm{eq})\}_{i=1}^3). \end{aligned}$$
(26)

If the maximum of local Lyapunov dimensions on the global attractors, which involves all equilibria, is achieved at an equilibrium point: \(\dim _\mathrm{L} u^{cr}_\mathrm{eq} = \max _{u \in K} \dim _\mathrm{L} u\), then this allows one to get analytical formula for the exact Lyapunov dimension.Footnote 10 In general, a conjecture on the Lyapunov dimension of self-excited attractor [19, 114] is that for a typical system, the Lyapunov dimension of a self-excited attractor does not exceed the Lyapunov dimension of one of the unstable equilibria, the unstable manifold of which intersects with the basin of attraction and visualizes the attractor.

To avoid numerical computation of the eigenvalues, one can use an effective analytical approach [19, 22, 115], which is based on a combination of the Douady-Oesterlé approach with the direct Lyapunov method: for example, in [22] for system (4) with \(b = 1\), it analytically obtained the following estimate

$$\begin{aligned} \dim _\mathrm{L} K \le 3 - \frac{2 (\sigma + 2)}{\sigma + 1 + \sqrt{(\sigma - 1)^2 + \frac{16\,r}{3\,\sigma }}}. \end{aligned}$$

The proof of the above conjecture and analytical derivation of the exact Lyapunov dimension formula for system (4) with all possible parameters is an open problem.

In [116, 117], it is demonstrated how the above technique can be effectively used to derive constructive upper bounds of the topological entropy of dynamical systems.

Fig. 6
figure 6

Localization of the hidden attractor of system (4) with \(r = 100\), \(a = -9.965\,\times \,10^{-3}\), \(b = 7.7454\,\times \,10^{-2}\) by the absorbing set \(\mathcal {B}^{h}\) with \(\delta = -a + 0.1\) and cuboid \(C^{h} =[-\,11, 11]\times [-17, 19]\times [80, 117]\)

Table 2 Numerical estimation of finite-time Lyapunov dimension in the case of hidden attractor (see Fig. 3)

6 The finite-time Lyapunov dimension in the case of hidden attractor and multistability

Consider the dynamical system \(\{\varphi ^t\}_{t\ge 0}\), generated by system (4) with parameters (5), and its attractor K. Here, \(\varphi ^t\big ((x_0,y_0,z_0)\big )\) is a solution of (4) with the initial data \((x_0,y_0,z_0)\). Since the dynamical system \(\{\varphi ^t_\mathrm{R}\}_{t\ge 0}\), generated by the Rabinovich system (1), can be obtained from \(\{\varphi ^t\}_{t\ge 0}\) by the smooth transformation \(\chi ^{-1}\), inverse to (2), and inverse rescaling time (3) \(t \rightarrow \nu _1 t\), we have \(\dim _\mathrm{L}(\{\varphi ^t\}_{t\ge 0},K)=\dim _\mathrm{L}(\{\varphi ^t_\mathrm{R}\}_{t\ge 0},\chi ^{-1}(K))\). In our experiments, we consider system (4) with parameters \(r = 100\), \(a = -9.965\,\times \,10^{-3}\), \(b = 7.7454\,\times \,10^{-2}\) corresponding to the hidden chaotic attractor.

In Fig. 6 is shown the grid of points \(C^{h}_\mathrm{grid}\) filling the hidden attractor: the grid of points fills cuboid \(C^{h} = [-11,11]\times [-17,19]\times [80,117]\) with the distance between points being equal to 0.5. The time interval is \([0,\, T=100]\), \(k=1000\), \(\tau =0.1\), and the integration method is MATLAB ode45 with predefined parameters. The infimum on the time interval is computed at the points \(\{t_i\}_{1}^{k}\) with time step \(\tau =t_{i+1}-t_i=0.1\). Note that if for a certain time \(t=t_i\) the computed trajectory is out of the cuboid, the corresponding value of finite-time local Lyapunov dimension is not taken into account in the computation of maximum of the finite-time local Lyapunov dimension (there are trajectories with initial data in cuboid, which are attracted to the zero equilibria and belong to its stable manifold, e.g. system (4) with \(x=y=0\) is \(\dot{z} = -bz\)). For the finite-time Lyapunov exponents (FTLEs) computation, we use MATLAB realization from [9] based on (21) with \(p=2\). For computation of the finite-time Lyapunov characteristic exponents (FTLCEs), we use MATLAB realization from [118] based on (23).

For the considered set of parameters, we compute:

  1. (i)

    Finite-time local Lyapunov dimensions \(\dim _\mathrm{L}(100,\cdot )\) at the point \(O_1 = (0, \, 1, \, 98)\), which belongs to the grid \(C^{h}_\mathrm{grid}\), and at the point \(O_2 = (0.0099, \, 0.0995, \, 0)\) on the unstable manifold of zero equilibrium \(S_0\);

  2. (ii)

    Maximum of the finite-time local Lyapunov dimensions at the points of grid, \(\max _{u \in C^{h}_\mathrm{grid}} \dim _\mathrm{L}(t,u)\), for the time points \(t=t_k=0.1\,k\) \((k=1,\ldots ,1000)\);

  3. (iii)

    The corresponding values, given by the Kaplan–Yorke formula with respect to finite-time Lyapunov characteristic exponents.

The results are given in Table 2. The dynamics of finite-time local Lyapunov dimensions for different points and their maximums on the grid of points are shown in Fig. 7.

For the absorbing set \(\mathcal {B}^{h}\) and the corresponding grid of points \(\mathcal {B}^{h}_\mathrm{grid}\) (the distance between grid points is 5), by estimation (25), we get the following estimate:

$$\begin{aligned}&\dim _\mathrm{H}K \le \dim _\mathrm{L} K \le \sup _{u \in \mathcal {B}^{h}}d_\mathrm{L}^\mathrm{KY}\big (\{\lambda _{j}(u)\}_{i=1}^3\big ) \nonumber \\&\quad \approx \sup _{u \in \mathcal {B}^{h}_\mathrm{grid}}d_\mathrm{L}^\mathrm{KY}\big (\{\lambda _{j}(u)\}_{i=1}^3\big ) = 2.97001... \, . \end{aligned}$$
(27)

Assuming \(\sigma + 1 \ge b\), the eigenvalues of the unstable zero equilibrium \(S_0\):

$$\begin{aligned}&\lambda _{1,3}(S_0) = -\tfrac{1}{2} \left[ (\sigma + 1)\mp \sqrt{(\sigma - 1)^2 + 4 \sigma r}\right] ,\nonumber \\&\quad \lambda _{2}(S_0) = - b, \end{aligned}$$

are ordered by decreasing, i.e. \(\lambda _{1}(S_0) > \lambda _{2}(S_0) \ge \lambda _{3}(S_0)\). Therefore, for the considered values of parameters by (26), we get

$$\begin{aligned} \dim _\mathrm{L} S_0= & {} d_\mathrm{L}^\mathrm{KY}(\{\lambda _i(S_0)\}_{i=1}^3) = 2 + \tfrac{\lambda _{1}(S_0) + \lambda _{2}(S_0)}{|\lambda _{3}(S_0)|} \nonumber \\= & {} 3 - \tfrac{2 (\sigma + b + 1) }{(\sigma + 1) + \sqrt{(\sigma - 1)^2 + 4 \sigma r}} =2.8111...\ .\nonumber \\ \end{aligned}$$
(28)
Fig. 7
figure 7

Dynamics of the finite-time local Lyapunov dimension estimations on the time interval \(t \in [0, 100 ]\): the maximum on the grid of points (dark red), at the point \(O_1 = (0, \, 1, \, 98) \in C^{h}_\mathrm{grid}\) (light red), at the point \(O_2 = (0.0099, \, 0.0995, \, 0)\) from the one-dimensional unstable manifold of \(S_0\) (blue). (Color figure online)

The above numerical experiments lead us to the following important remarks. While the Lyapunov dimension, unlike the Hausdorff dimension, is not a dimension in the rigorous sense [119] (e.g. the Lyapunov dimension of the saddle point \(S_0\) in (28) is noninteger), it gives an upper estimate of the Hausdorff dimension. If the attractor K or the corresponding absorbing set \(\mathcal {B} \supset K\) is known [see, e.g. (8)] and our purpose is to demonstrate that \(\dim _\mathrm{H} K \le \dim _\mathrm{L} K < 3\), then it can be achieved without integration of the considered dynamical system [see, e.g. (27)]. If our purpose is to get a precise estimation of the Hausdorff dimension, then we can use (18) and compute the finite-time Lyapunov dimension, i.e. to find the maximum of the finite-time local Lyapunov dimensions on a grid of points for a certain time. To be able to repeat a computation of finite-time Lyapunov dimension, one need to know the initial points of considered trajectories \(\{u_i\}_{i=1}^{k}=K_\mathrm{grid}\) on the set K, time interval \((0,T] = \bigcup _{i=0}^{k-1} (t_i,t_{i+1}]\) and how the finite-time Lyapunov exponents were computed (e.g. by (23), (21), or (22) with the parameter \(\delta \)).

7 Computation of the Lyapunov dimension and transient chaos

Consider an example, which demonstrates difficulties in the reliable numerical computation of the Lyapunov dimension (i.e. numerical approximation of the limit values of the finite-time Lyapunov dimensions).

Consider system (4) with parameters \(r = 6.485\), \(a = -0.5\), \(b = 0.85\) for which equilibrium \(S_0\) is a saddle point and equilibria \(S_\pm \) are stable focus-nodes and integrate numericallyFootnote 11 the trajectory with initial data \(u_{0} = (0.5,\,0.5,\,0.5)\) in the vicinity of the \(S_0\). We discard the part of the trajectory, corresponding to the initial transition process (for \([0,t_\mathrm{tp} = 25{,}000]\)), and get the point \(u_\mathrm{init} = (\) \(-2.089862710574761\), \(-2.500837780529156\), 2.776106323157132). Further, we numerically approximate the finite-time Lyapunov exponents and dimension for the time interval [0,  T] by (23) (see MATLAB code in [120]).

Fig. 8
figure 8

The trajectory forms a chaotic set, which looks like an “attractor” (navy blue) and then tends to \(S_+\) (cyan). a Trajectory \(u(t,u_\mathrm{init})\) for \(t \in [0,~T_1]\),   \(T_1 \approx 15{,}295\). b Trajectory \(u(t,u_\mathrm{init})\) for \(t \in [0,~T]\),   \(T = 500{,}000\). (Color figure online)

Fig. 9
figure 9

Numerical computation of \({{\mathrm{LCE}}}_1(t, \, u_\mathrm{init})\) and \({d}_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LCE}}}_i(t, \, u_\mathrm{init})\}_{i=1}^3)\) for the time interval \([0, T_1\approx 15{,}295]\). a \({{\mathrm{LCE}}}_1(t, \, u_\mathrm{init})\), \(t \in [0,~T_1]\),   \(T_1 \approx 15{,}295\). b \({d}_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LCE}}}_i(t, \, u_\mathrm{init})\}_{i=1}^3)\), \(t \in [0,~T_1]\),   \(T_1 \approx 15{,}295\)

Fig. 10
figure 10

Numerical computation of \({{\mathrm{LCE}}}_1(t, \, u_\mathrm{init})\) and \({d}_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LCE}}}_i(t, \, u_\mathrm{init})\}_{i=1}^3)\) for the time interval \([0,\, 5 \cdot 10^5]\). a \({{\mathrm{LCE}}}_1(t, \, u_\mathrm{init})\), \(t \in [0,~T]\),   \(T = 5 \cdot 10^5\). b \({d}_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LCE}}}_i(t, \, u_\mathrm{init})\}_{i=1}^3)\), \(t \in [0,~T]\),   \(T = 5\,\cdot \,10^5\)

The trajectory computed on the time interval \([0, T_1\approx 15295]\) traces a chaotic set in the phase space, which looks like an “attractor” (see Fig. 8a, Fig. 9a). Further integration with \(t > T_1\) leads to the collapse of the “attractor”, i.e. the “attractor” (see Fig. 8a) turns out to be a transient chaotic set (see Fig. 8b). However on the time interval \(t \in [0,~T_3\approx 431{,}560]\), we have \({{\mathrm{LCE}}}_1(t, \, u_\mathrm{init}) > 0\) (see Fig. 10a) and, thus, may conclude that the behaviour is chaotic, and for the time interval \(t \in [0,~T_2\approx 223{,}447]\), we have \({d}_\mathrm{L}^\mathrm{KY}(\{{{\mathrm{LCE}}}_i(t, \, u_\mathrm{init})\}_{i=1}^3) > 2\) (see Fig. 10b). 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 numerical computation of the finite-time Lyapunov exponents and the finite-time Lyapunov dimension does not necessarily lead to a relevant approximation of the Lyapunov exponents and the Lyapunov dimension [see also effects in (24)].

8 Conclusion

In this work, the Rabinovich system, describing the process of interaction between waves in plasma, is considered. It is shown that the methods of numerical continuation and perpetual point are helpful for the localization and understanding of hidden attractor in the Rabinovich system. For the study of dimension of the hidden attractor, the concept of the finite-time Lyapunov dimension is developed. An approach to reliable numerical estimation of the finite-time Lyapunov exponents [see (21) and (22)] and finite-time Lyapunov dimension [see (18)] is suggested. Various numerical estimates of the finite-time Lyapunov dimension for the hidden attractor in the case of multistability are given.