On the Poisson Stability to Study a Fourth-Order Dynamical System with Quadratic Nonlinearities

This article discusses the search procedure for the Poincar\'e recurrences to classify solutions on an attractor of a fourth-order nonlinear dynamical system using a previously developed high-precision numerical method. For the resulting limiting solution, the Lyapunov exponents are calculated using the modified Benettin's algorithm to study the stability of the found regime and confirm the type of attractor.


Introduction
As it is known [1], the calculation of the Lyapunov exponents is used for the classification of attractors of dynamical systems. The combinations of signs of such values determine the attractor type: an equilibrium position, limit cycle, multidimensional torus or strange attractor. The dynamics of the corresponding solutions are called static, periodic, quasiperiodic, or chaotic. The computational procedure for determining the Lyapunov exponents is mainly based on the Benettin's algorithm [2]. However, in [3,4] the points of limiting solutions were investigated for the Poisson stability, which made it possible to understand whether we have a quasiperiodic or chaotic regime. In [4], it was done for the Jafari-Sprott system [5].
It is noteworthy [6] that a point y of the phase space is called positively Poisson stable (P + ), if for any neighborhood U of the point y and for any T P > 0 there is a time value t ≥ T P for which the trajectory of the dynamical system enters into the neighborhood U. Similarly, if there is t ≤ −T P such that the trajectory enters into the neighborhood U, then the point y is negatively Poisson stable (P − ). A point is called Poisson stable, if it is P + and P − -stable.
The boundedness of the limiting solutions of dissipative systems implies [6][7][8][9][10] that any steady-state oscillation mode is described by Poisson-stable trajectories. This also applies to dynamical chaos. A trajectory different from the equilibrium position is said to be Poisson stable if it verifies the property of returning in an arbitrarily small ε-neighborhood of each of its points an infinite number of times. Such returns are called the Poincaré recurrences. In [10, p. 146], the author notes that "the study of the statistics of the Poincaré recurrences is a powerful tool for the analysis and classification of dynamic modes. Apparently, the potentialities of this approach have not yet been fully exhausted in modern nonlinear dynamics". For example, the returns follow one another regularly for quasiperiodic regimes. Then [10, p. 145] "the dynamic chaos is a situation when the Poincaré recurrences to the ε-neighborhood of the initial point do not show regularity, the time interval between two successive returns turns out to be different each time, and some statistical distribution of the times of return is arised" [11,12]. An example of the Poincaré recurrences analysis based on the Kac's theorem [13][14][15][16][17][18], for a discrete dynamical system with a chaotic non-hyperbolic attractor is given in [16].
As it is known, for the most of dissipative systems, the formulas for a general solution of the system in a class of any known functions with well-studied properties have not yet been found. Therefore, numerical methods are used. The use of classical numerical methods (such as Euler method, Runge-Kutta 4th order method, Adams methods, etc.) for constructing approximate solutions in attractors of dynamical systems leads to significant errors over large sections of time due to the instability of the studied chaotic regimes. In general, the problem of numerical modeling with control of the accuracy of obtained solution and the choice of a platform for computer implementation is relevant today [19], since small errors introduced at each integration step cause exponential divergence of close trajectories. More recently, the numerical FGBFI method (Firmly Grounded Backward-Forward Integration) has been proposed in [3,4,20], based on the power series method for dynamical systems with quadratic nonlinearities. The main advantages of this method are as follows: 1.
The recurrence relations are obtained for calculating the coefficients of the expansion of solutions in a power series for any dynamical system with quadratic nonlinearities in a general form; 2.
The convergence of the power series is studied. A simple formula for calculations is derived (in comparison with that [21] obtained in the known literature) for calculating the length of the integration step in a general form; 3.
The criteria for checking the accuracy of the approximate solution are obtained. The control of the accuracy and configuration of the approximate solution of a dynamical system uses forward and backward time, which makes the numerical method reliable (degrees of piecewise polynomials, the value of the maximum integration step, etc.); 4.
The FGBFI method allows to construct high-precision approximations to non-extendable solutions of a system of autonomous differential equations with a quadratic right-hand side. Like for instance the system of the forṁ In this case, the numerical solution computed with FGBFI will never cross the asymptote and will approach it arbitrarily close.
Thus, it is a high-precision method for constructing the trajectory arc of the system for any time interval. This makes it possible to track the Poincaré recurrences to any neighborhood of the trajectory points, since the resulting computational error can be smaller than the radius of the monitored neighborhoods.
This article considers a fourth-order dissipative system [22], in which, in the opinion of its authors, there is a hyperchaos. The goals of the research are: (1) to find the Poincaré recurrences in the attractor of this system for the approximate solutions obtained by the FGBFI method and to analyze the statistics of return times; (2) to apply the modified Benettin's algorithm [20] to refine the computation of the Lyapunov exponents.

Bohr's almost periodic functions
It is noteworthy [23, pp. 368, 418, 419] that the function (trajectory of a dynamical system) F(t) ∈ R n is called Bohr's almost periodic in the sense that if for any ε > 0 there exists a relatively dense set of almost periods τ F = τ F (ε) of the function F(t) with ε-accuracy, i.e. there is the positive number L = L(ε) such that any interval [α; α + L] contains at least one number τ F for which the inequality |F(t + τ F ) − F(t)| < ε holds at t ∈ R. Note that for an almost periodic function (and, in general, for a Poisson stable trajectory), different from a periodic function, with decreasing number ε the number L(ε) must increase indefinitely, otherwise the function would be periodic.
Usually (see, for example, [10]) in nonlinear dynamics, almost periodic functions are called quasiperiodic. Note an important theorem in the theory of dynamical systems [10, p. 147]: if a non-periodic trajectory is Poisson and Lyapunov stable, then it is almost periodic.
Since we will be investigating the Poincaré recurrences, it is noteworthy that for an almost periodic trajectory (as follows from the above definition), such returns must form a relatively dense set. An example of a sequence that is not a relatively dense set is: 0, 1 2 , 2 2 , 3 2 , . . . , i.e. in the given example, the distance between neighboring elements grows with the growth of q. A similar situation was observed in a numerical experiment for the Poincaré recurrences in the case of a chaotic solution of the Chen system [3].

Dissipative system of the fourth order
Let us consider a fourth-order system [22]       ẋ where a, b, c, d, e, f and k are positive parameters. The divergence of the vector field G defined by the right-hand side of the system (1) is equal to div G = −(a + 1 + c + d) < 0.
Then the system (1) is dissipative. Hence, there is the ball B a ⊂ R 4 , into which each trajectory of the system (1) is immersed forever, after a while. Therefore, there is a limit set (attractor) such that all trajectories of the dynamical system are attracted when t → ∞ [6]. Let us investigate the dynamics of the system (1) for the values of the parameters a = 7, b = 50, c = 3, d = 10, e = 5, f = 5 and k = 1.5. In this case, it is possible to determine the position and radius of the ball B a from the data obtained in [22]. Also, in this article the Lyapunov exponents were determined: LE 1 = 2.1040, LE 2 = 0.3563, LE 3 = 0 and LE 4 = −23.4508, indicating a hyperchaos in the system (1). Therefore, its solutions are highly unstable, which requires the use of special numerical procedures for large time intervals.
Let us consider the application of the FGBFI method for constructing approximations to the solutions of the system (1) using high-precision calculations.

The FGBFI method
Following [4], we will rewrite the system (1) asẊ where Let us expand the solution where X(0) = Υ 0 is an initial condition for the system (2), the vector Υ 0 is given, The recurrence relations for calculating the coefficients of the power series have the form [4] for j ∈ N.
Since the criteria for checking the accuracy of the obtained numerical solution require repeated calculations in forward and backward time, we need to have a guaranteed estimate of the convergence interval of the power series (3) for a given vector Υ 0 . In [4], a theorem is proved on the estimate of this interval for systems with a quadratic right-hand side.
The norms are calculated: Next, two auxiliary numbers are calculated as functions of Υ 0 : δ is any positive number (can be taken very small). It can be seen from the formula (5) that the length of the convergence interval depends on the choice of the initial condition for the system (2). This fact is confirmed by computational experiments, for example, for the Lorenz system [24].
It is worth noting that a similar approach can be applied to systems with cubic nonlinearities, but when the derivative on the left side of the system is of the second (or higher) order. For example, in [25] this was done for the Duffing equation.
Usually, in calculations, many researchers work with subnormal real numbers of single or double precision, presented in the IEEE 754 format [26]. The main drawback here is the fixed accuracy of the representation of real numbers, which may not allow us numerically constructing approximations to unstable solutions of systems of differential equations on the large time intervals. Therefore, for high-precision calculations, the GNU MPFR library [27,28] is used.
Next, consider the algorithm [20] for the numerical construction of an approximate solution of a dynamical system in forward and backward time:

1.
Set the number b m of bits under the mantissa of a real number and the precision ε pw for estimate of the common term of power series. Note that b m defines the machine epsilon ε m . Let us choose b m so that the precision of representation of the real number is with a margin, i.e. ε m ε pw ; Set X(0) ∈ B a for the system (2) and value the number way that determines the direction in time: way = 1 is gone forward in time, way = −1 is gone backward in time; 4.
Set T: the length of the time interval on which the numerical integration will be performed; 5.
As can be seen, the algorithm is universal in the direction of numerical integration in time. Next, we consider the criteria for verifying the accuracy of the resulting solution. For this we assume where l is an index of the interval [t l−1 , t l ], on which the series (3) is converged, N {way} is a number of such intervals for the direction way on time. Since on each interval [t l−1 , t l ] during calculations we truncate the series (3) to some polynomial X l (t) approximating the solution of the system (2) on it (i.e. locally), then we denote by n {way} l the degree of the polynomial X l (t). Then we introduce the notation: In [4], the following criteria for checking the accuracy of the obtained solution were proposed: 1.
The accuracy ε a of approximation at way = 1 is a frequently used criterion in applications of numerical methods for solving differential equations. When the inequality (6) is true, it is necessary to increase the degrees of all polynomials n {1} l , obtaining the next approximation, and compare the distance δ a between the obtained approximate solutions on the interval [0, T] with the value ε a . If δ a > ε a , then we increase the powers of n l , otherwise we have to use the obtained solution;

2.
The radius ε R of the neighborhood of the initial point, to which the approximate solution should return in backward time is another criterion. In other words, we need to select the accuracy of ε pw so that the following inequality holds where N = N {−1} . Note here that we do not know the exact solution to the system (2), but we do know its initial point. Then we can set how many digits are in each component of the vector X after the point must match the digits of the corresponding components X(0). The problem is that a large amount of computation is required (due to repetitions of the forward and backward traverse in time). The solutions of the system (2), as a rule, are highly unstable in the backward time: they immediately leave from the attractor, because in our calculations, we are near it, no matter how accurate ε pw is taken. Therefore, in the above algorithm, we control the finding of the trajectory within the boundaries of the set B a ; 3.
The comparison of configurations of approximate solutions in forward and backward time is necessary to determine the numbers (7), describing approximate solutions. Next, check Unlike criteria 1 and 2, we control here the arguments of the approximate solutions.

Analysis of the Poincaré recurrences for a fourth-order system
In [22, p. 3221], the authors indicated that there is a hyperchaos for the considered parameters of the system (1). For example, the authors found it for the initial condition Note that in the numerical experiment by to the first criterion for checking the accuracy, decreasing the value of ε pw , the number of indicated correct signs in coordinates is preserved.
Next, we take the point (9) as the initial point and track the Poincaré recurrences for it. In order to reduce the amount of calculations without moving backward in time (so as not to take too small a value of ε pw ), we look at the times when there is a return to the neighborhood of the initial point for different values ε pw . This is necessary to control the accuracy of the obtained points, choosing such the value ε pw that these moments of time do not differ significantly.
To determine the Poincaré recurrences, we choose the small time step ∆t P > 0. Let us divide the time interval [0, T P ], where T P is the end time at which the search for the Poincaré recurrences is terminated, into intervals of length ∆t P (let N P be the number of such intervals). Let us denote X k the coordinates of the k-th point of the considered trajectory arc, corresponding to the time moment t k = k∆t P , k = 1; N P .
We introduce the distance d k = |X k − X(0)|.
To track the Poincaré recurrences, we fix such values k = k * , moving in the trajectory arc, when there are the local rapprochements with the point X(0), i.e.
For the point (9) for T P = 10 and ∆t P = 10 −4 , the returns are shown in Tabl. 1. As can be seen from this table, the rapprochement with the initial point occurs at approximately the same time (i.e., we have regularity), which suggests the idea of a limit cycle in the system (1). Decrease in the value of ∆t P does not cause significant changes in Tabl. 1.
In this case to show the absence of chaos in the system (1), we take the time value T = 40 (instead of T = 15, which was used to obtain the coordinates of the point (9)) to reach the attractor. We get the following point at the end time: Similarly, for the point (10), the Poincaré recurrences shown in Tabl. 2 for T P = 10 and ∆t P = 10 −4 . As we can be seen from this table, returns are getting closer. Therefore, we will decrease the value ∆t P , setting it equal to ∆t P = 10 −7 . We get the returns shown in Tabl. 3. From this iterative procedure, we see that the limiting trajectory in this case is a cycle with a period T c ≈ 0.3655. Its projections are shown in Figs. 1 and 2.
To select the value ∆t P , this step must be reduced until the rapprochement distance d k * changes significantly. This approach was also applied for other initial conditions. As a result, we get the same cycle. Next, we investigate the stability of the resulting cycle by calculating the Lyapunov exponents.

Calculation of the Lyapunov exponents
In [20], a modification of the algorithm for calculating the Lyapunov exponents for a model describing the growth of a cancerous tumor was considered. Next, we will consider a generalization of this algorithm for dynamical systems with a quadratic right-hand side of the n-th order. Note that an other algorithm for calculating the Lyapunov exponents based on calculating the fundamental matrix of the dynamical system and its QR decomposition is proposed in [29]. The procedure described below can also be applied to it.
First, we define the form of the linearized system of equations. Let us denote the perturbations by We define the form of the following expression using the rule for differentiating the sums and products: where p = 1; n, j = 1; n. Let us find the product Next, we extend the system (1) by adding to it a linearized system, which also has a quadratic right-hand side. To do this, we introduce the extended matrixÃ: The vector functionX(t) is equal toX Based on the formula (11), we have a quadratic form on the right-hand side of the linearized system for each equation. Then, to reduce the system to general form, we introduce the following extended matrix: LE p := 0, p = 1; n (the initial values of the sums at calculating the Lyapunov exponents);
Note that for all values of p the vector Y (k) is the same. End of cycle 11. For p from 1 to n with step 1: Beginning of cycle End of cycle 12. If k < M, then Goto step 9; 13. LE p := LE p M · τ M , p = 1; n; 14. Print LE p , p = 1; n.
For the point (10), the Lyapunov exponents were found using the described algorithm for T = 100, M = 20000 and τ M = 0.005. Next, we give them the approximate values: As can be seen, their signature corresponds to a stable limit cycle in the system (1), which confirms the presence of the cycle, found through the study of the Poincaré recurrences. Note that in [30] an alternative numerical-analytical scheme for constructing periodic solutions of systems with a quadratic right-hand side was proposed using an example of the Lorenz system. Now this scheme is actively developing in [31][32][33] to find periodic solutions of nonlinear systems of differential equations.
The graphs of changes in the Lyapunov exponents over time are shown in Figs. 3 -6.
The source codes of all computing programs are located at [34].

Conclusion
In this article, the author investigates the limiting regimes corresponding to different initial points of the phase space. As a result, we get the same cycle shown in Fig. 1 and 2. Consequently, there is no hyperchaos in system (1).
The analysis of the Poincaré recurrences, described in the goals of research, showed that the rapprochement with the initial point occurs approximately at the same time (i.e., we have a regularity). This indicates a limit cycle in the system (1). Figs. 3 -6 indicates the correctness of finding their limiting values over a relatively short time interval [0, 100].

Stabilization of the behavior of the values of the Lyapunov exponents in
Thus, the limiting trajectory investigated using the analysis of the Poincaré recurrences and the high-precision calculation of the Lyapunov exponents (i.e., different methods) is the limit cycle which contradicts the results in [22] for the values of the parameters a = 7, b = 50, c = 3, d = 10, e = 5, f = 5 and k = 1.5 of the system (1). The obtained data is available at [34]. Most likely, the authors of the article [22] made a typo indicating such values of the parameters at which they found a hyperchaos.
A significant novelty of the article is the high-precision method for searching for the Poincaré recurrences on attractors of dynamical systems with a quadratic right-hand side, also a modification of Benettin's algorithm for finding the Lyapunov exponents in order to check the found limiting regimes.
The methods described in the article can be used to study not only periodic regimes, but also quasiperiodic and chaotic regimes for dynamical systems with a quadratic right-hand side. According to the statistics of Poincaré recurrences, we can determine the dimension of the attractor and compare it with the dimension obtained when calculating the Lyapunov exponents (for example, the articles [17,18]).
The advantage of the considered algorithm of calculating the Lyapunov exponents is combination of the linearized system of equations and the original dynamical system in general form (12) for finding the state and perturbation vectors together.

Acknowledgments
The reported study was funded by RFBR for the research project No. 20-01-00347.