On a High-Precision Method for Studying Attractors of Dynamical Systems and Systems of Explosive Type

The author of this article considers a numerical method that uses high-precision calculations to construct approximations to attractors of dynamical systems of chaotic type with a quadratic right-hand side, as well as to find the vertical asymptotes of solutions of systems of explosive type. A special case of such systems is the population explosion model. A theorem on the existence of asymptotes is proved. The extension of the numerical method for piecewise smooth systems is described using the Chua system as an example, as well as systems with hysteresis.


Introduction
Let us consider a dynamical system with quadratic nonlinearities of the forṁ where X(t) = [x 1 (t) . . . x n (t)] T is a vector function of time t with values in space R n , B 0 ∈ R n is a given column vector, ϕ(X) = [ϕ 1 (X) . . . ϕ n (X)] T , ϕ p (X) = Q p X, X , B 1 and Q p (p = 1, . . . , n) are the matrices (n × n) of real numbers. There has been a strong interest in the systems of the form (1) in scientific literature starting somewhere in the middle 20th century. The most popular object of study from the class of systems (1) is the Lorenz system [1]: (2) For this system, the matrices have the form: Such systems are of interest to researchers because they are systems with energy dissipation. There is compression of the volume of the phase space in them, and all trajectories are uniformly bounded. Thus, there is an attracting set called the attractor of the dynamical system to which all solutions tend at t → ∞. Consequently, the attractors of dissipative systems determine their behavior over long time intervals. However, for many systems of the form (1), the limit solutions are unstable. Therefore, correct numerical integration requires the high-precision methods that make it possible to flexibly control the accumulation of calculation errors over a long time interval. For example, with an increase of the order of the Runge-Kutta methods, the complexity of the formulas for calculating the approximate values of the phase coordinates (they have a tree structure) [2,3] grows.
Recently, to simulate of dynamic chaos in systems of the type (1), some researchers have begun to use effective modifications of the power series method. Gibbons in his article [4] proposed the quick formulas for calculating the coefficients of power series for the main types of right-hand sides of differential equations in the scalar form. Today this idea was generalized in a recursive procedure (called as automatic differentiation) to compute the values of the derivatives for power series [5]. An advantage over the general Taylor series method is that the calculations can be constructed by fast formulas in comparison to the direct symbolic differentiation of right-hand sides of nonlinear ODEs which requires a lot of computer memory for high-precision calculations. The method of power series in [6][7][8] is applied as the Adomian decomposition method (ADM). The Clean Numerical Simulation (CNS) [9] is based on the method of power series at arbitrary-order and used the multiple-precision data, plus a check of solution by means of an additional computation using even smaller numerical noises.
There is the FGBFI numerical method (the firmly grounded backward-forward integration method) described in [10][11][12][13]. In this articles, the recurrence relations for calculating of the coefficients of expansion of local solutions in a power series are received for any dynamic system with quadratic nonlinearities in the vector form (it distinguishes the FGBFI-method from automatic differentiation and CNS). This form of coefficients of the power series is simpler and faster to compute than in ADM (because it does not contain factorials). The author of this article derived the simple formulas of calculating the length of the integration step in the general form (it distinguishes the FGBFI-method from CNS). The criteria for checking the accuracy of the approximate chaotic solution are obtained.
This article will consider this method which uses high-precision calculations to construct approximations to unstable solutions of the system (1), as well as its extension to piecewise smooth systems composed of quadratic nonlinearities (for example, the Chua's piecewise linear system).
It is also noteworthy that for certain types of matrices B 1 and Q p and the vector B 0 , the system (1) can have non-extendable solutions of the explosive type. For example, for scalar systems of the forṁ with a particular solution x 1 (t) = tan t one of the asymptotes is The existing classical numerical methods available in the standard implementation of such mathematical packages including Mathcad, MATLAB, Maple, etc., cannot be used to build approximations to such solutions, since the numerical scheme will jump over the asymptote, "not feeling" it. But as will be shown in this article, the FGBFI method allows to get arbitrarily close to the asymptote without skipping the time t as , i.e., localizing it.

The FGBFI Method
The instability of the limiting solutions of the system (1), in fact, makes incorrect to use classical numerical methods over long time intervals. The limit sets of dynamical systems are constructed on such segments. We can use high-precision calculations but then we will have rather small integration steps. The solution to the problem is to apply a modification of the power series method with recurrent calculation of expansion coefficients.
Let us represent the solution of the system (1) as a series where the vector Υ 0 is determined by the given initial condition. The remaining expansion coefficient vectors (5) are calculated by the following recursive formulas [11]: . . , n, The convergence region [−τ(Υ 0 ), τ(Υ 0 )] of the series (5) is limited for many nonlinear systems of the form (1), depends on the initial condition, and can be estimated by the following formulas [11]: where δ is any positive number (can be chosen arbitrarily small). In these formulas, one of the compatible matrix norms · 1 or · ∞ with a vector norm is meant. The recursive formulas (6) for calculating the coefficients of the power series expansion were obtained from the substitution of power series into the system (1), as well from the products in Cauchy form. A theorem on the convergence of a power series is proved in the article [11]. The proof is based on the method of mathematical induction. The specified region from formulas (7) is obtained as a corollary from this proof.
Next, we describe an algorithm for constructing a trajectory arc on a given time interval [0, T].
In the papers [12,13], the numerical method based on the representations (5)- (7) is called the FGBFI method. Next, we describe the algorithm underlying it.
Let T is the given length of the integration interval, S a is the ball containing the attractor of the system (1). Let us set such a representation of a real number that the machine epsilon where ε pw is an accuracy of evaluation of the common term of the series (5). Thus, the summation by the Formula (5) stops at such a value i = i * when where ∆t is an integration step. Note that for the convergence of the series, the value ∆t must be chosen as Since the convergence region of the series (5) is limited, we cannot take ∆t = T (as mentioned above, the value of T is large). Therefore, the trajectory arc on the time interval [0, T] will consist of the connected arcs for which the series (5) converges.
Since in the computer simulation we can only operate with a discrete set of phase coordinate vectors (instead of a continuous function X(t)), then, by analogy with the Runge-Kutta methods, we will obtain approximate values of phase coordinates through an integration step, only calculations will be carried out at each step using Formulas (5)- (8) for Note that here the error accumulation at each step is much smaller than for the fixed-order Runge-Kutta methods, since we can quickly enough calculate the missing terms in the sum using the Formulas (5) and (6) by reducing the value of ε pw . It can be seen from the formulas (7) that with such a numerical integration, the step will be variable.
Let the value way is the direction of integration: way = 1 for forward time, way = −1 for reverse time.
Next, we show an algorithm for obtaining approximate values of phase coordinates.

Algorithm 1
The FGBFI Method 1: Set the vector X 0 ∈ S a of initial condition and the value way 2: ct := 0 The current time 3: ended := false Are we finishing the algorithm? 4: l := 1 Number of closed interval of series convergence (5) 5: Calculate the value of integration step ∆t = τ(X l−1 ) using the formulas (7) 6: if ∆t > T − ct then 7: ∆t := T − ct 8: ct := T 9: else 10: ct := ct + ∆t 11: end if 12: Υ 0 := X l−1 13: X l := X l−1 14: p := 1 The product of powers of t 15: i := 0 Number of the current member of the series 16: do 17: i := i + 1 18: p := p · way · ∆t Calculate the current degree of t = ∆t taking into account the direction of integration; p is negative for odd powers when way = −1

19:
Calculate the vector Υ i using the formulas (6) 20: Add to X l the current member of the series 21: The right side of inequality (8) 22: while L > ε pw 23: Print the value of time way · ct and the vector X l of the received phase coordinates 24: if X l / ∈ S a then 25: We moved beyond the ball S a The criteria for checking the accuracy of the resulting approximate solution are described in [10,11,13]. Note that the backward time pass is used in some of them.

The Software Implementation of the FGBFI Method and Numerical Simulations
For the software implementation of Algorithm 1 with and for the study of the limit points of dynamic systems of the form (1) for Poisson stability and the calculation of the Lyapunov exponents by the modified Benettin's algorithm [12,13], a software in the C++ language for Linux was developed. The source codes are available on GitHub [14].
To perform matrix operations in C++, the algorithms and template types of the uBLAS library of the linear algebra collection of the Boost class libraries were used. Note that the uBLAS library gives a large gain in time.
Usually, in calculations, many researchers work with single or double precision subnormal real numbers represented in the IEEE 754 format. The main disadvantage here is the fixed accuracy of the representation of real numbers which may not allow us to numerically construct approximations to unstable solutions of systems of differential equations over long time intervals. Therefore, for high-precision calculations, the MPFR C++ library [15] is used.
Algorithm 1 was applied to check the accuracy of the found unstable cycle [16] in the system (2): The period value is obtained equal to T c ≈ 1.558652210.
The initial values (9) were checked on the period in the computer program [14] that implements the FGBFI method with ε pw = 10 −25 , 100 bits for mantissa of real number and ε m = 1.57772 · 10 −30 . With such parameters of the method, the approximate values x 1 (T c ), x 2 (T c ) and x 3 (T c ) were also verified by the same numerical method, but in backward time. The values at the backward time in the end of the trajectory arc coincide with (9) up to the 9th character inclusive after the point. The resulting values of x 1 (T c ), x 2 (T c ) and x 3 (T c ) coincide with (9) up to the 8th character inclusive.
The cycle corresponding to (9) is shown in Figure 1. D. Viswanath [17] found in the system (2) an unstable cycle with a long period and with 99 decimal places use the Lindstedt-Poincaré technique: The author of this article decided to check the found values with high accuracy in the forward time by the FGBFI method. For this, it was taken 390 bits for mantissa of real number. Then ε m = 7.93107 · 10 −118 . The value of accuracy for the power series is ε pw = 10 −110 . The maximum degree of the approximating polynomial at variable integration steps is 78. The computation time was approximately 6.2 min.
After the computational experiment, it was found that at the end of the period the values x 1 (T c ), x 2 (T c ) and x 3 (T c ) coincide up to 38 decimal places with (10). At the same time, the increase in accuracy (decreasing values ε m and ε pw ) does not affect the obtained results where we can conclude that some decimal places of the values (10) are incorrect in [17], or need more correct decimal places.
Note that the classical numerical methods (such as the Euler method, Runge-Kutta methods, Adams methods, etc.) do not make it possible to carry out such a high-precision check of the values (10) due to the instability of the found cycle by D. Viswanath.
We also note that the use of high-precision calculations is important in data encryption. For example, the results obtained in [18] show that when using the high-precision arithmetic, the generated sequences provide good randomness and security with more iterations of the implemented discrete-time chaotic systems compared to the received sequences in the IEEE 754 format.

The Systems of Explosive Type and Application of the FGBFI Method
The dynamic system of the explosive type is a system for which the state X(t) → ∞ is reached in a finite time. The examples of such systems are the system (3), and the population explosion model where the coefficient k > 0. In other words, in this section of the article we will consider the autonomous systems of differential equations whose solutions have the vertical asymptotes. Such systems have been little-studied in the literature. We note the papers [19,20] where the scalar case is described: the asymptote existence theorem is proved and a modification of the Euler method is given which allows one to localize the asymptotes of solutions without jumping them.
For the multidimensional case, no studies have been found in the known literature. Most likely, this is due to the fact that for these systems, at X(t) → ∞, the energy also tends to infinity (which cannot be in real systems according to the energy conservation law); however, the systems of explosive type can be models of real systems in some approximation. Therefore, it is important for the researcher to have an idea at what time t as a sharp increase of phase coordinates will occur when it is not possible to find solutions of the system in quadratures.
Before proving the theorem on the existence of vertical asymptotes, we prove the following statement: Lemma 1. Let the number sequence r (1) , r (2) , ... converge to some number r, the sequence w (1) , w (2) , ... converge to number w, and there is the number for all m ∈ N, then r > w.
Proof. Let us represent the sequences as r (m) = r + ζ (m) , where ζ (m) and ϑ (m) are the infinitesimal sequences. From the inequality (11) it follows that Note that the sequence is also infinitesimal. Therefore, there exists a number m = m * such that for all m > m * moreover, we choose the small positive number ε such that because the number d is fixed. Then from (12) whence it follows that r > w.
Note that the proof of this lemma can be carried over to the case when the sequences r (m) (t) and w (m) (t) are the sequences of scalar functions that converge uniformly over some time interval.
Let us now prove a theorem that allows us to establish the existence of the vertical asymptotes in the system (1) at t > 0.

Theorem 1.
If some equation with the number q of the system (1) can be represented as (k is some number of the phase coordinate x k , k = q)ẋ where the symmetric matrixQ q is positive definite (D + Q ) or negative definite (D − Q ) and is obtained by deleting a row and a column with number k from matrix Q q , the initial conditions x k (0) = 0, x q (0) > a q,q for D + Q or x q (0) < a q,q for D − Q , then the q-th component of the solution of the system (1) has a vertical asymptote at t > 0.
Proof. For definiteness, we will consider the situation D + Q . A similar proof can be carried out for D − Q . Then the inequality ( [21], pp. 34, 35) is valid where λ q > 0 is the smallest eigenvalue of the matrixQ q . Since Let us introduce the equationẏ = c q + λ q (y − a q,q ) 2 (15) with the initial condition y(0) for which By the existence and uniqueness theorem, the solution of the Equations (13) and (15) exists and is unique on some time interval [0, t * ].
Let us construct the integral equations equivalent to the Equations (13) and (15): Since the function x k (t) is continuous on the interval [0, t * ], x k (0) = 0 and υ > 0, then is a maximum value of the function R(t) on the segment [0, t * ]. Let us write the function x q (t) as Let us now consider the schemes of successive approximations obtained from the Formula (17): where the vector functionX (m) (η) is obtained fromX(η) by replacing the function Since the existence and uniqueness theorem is proved by the method of successive approximations, then for the constructed approximations lim m→∞ x (m) By virtue of the inequality (14), at the first iteration we have From the positiveness of the function R(t), the inequalities (16) and we have i.e., x Also from the the inequalities (16) and (18), we have . In a similar way, one can show that At that max Then, by virtue of the proved lemma at t ∈ [0, t * ].
Note that the solution to the Cauchy problem (15)- (16), which is easy to find in quadratures, has a vertical asymptote. Then it follows from the inequality (19) that the function x q (t) also has a vertical asymptote.
Let us consider an example of a system that satisfies the conditions of this theorem. In the article [22], the following system was introduced: in which the equilibrium positions lie on the circle x 2 1 + x 2 2 = r 2 in the plane x 3 = 0, a is the parameter of the dynamical system. Usually the functions f 1 and f 2 are chosen as where b, β and γ are also system parameters.
Let r = 0, β = 0 and where α > 0 is a parameter. We get The system (21) belongs to the class of systems (1) and has the unique equilibrium (0, 0, 0). For the last equation of the system the matrixQ 3 is symmetric and positive-definite. Then the component x 3 (t) of the solution of the system (21) has a vertical asymptote for x 3 (0) > 0 and x 1 (0) = 0. Note that the authors of the article [22] do not provide a qualitative analysis of the system (20). As shown above, it may happen that it has unbounded solutions. Therefore, before looking numerically for attractors in dynamical systems, the researcher needs to show the existence of bounded solutions, for example, using the Lyapunov functions.
Note that by the proof of the above theorem, it is not necessary that all equations of the system (1) have a quadratic right-hand side. It suffices to have at least one equation that satisfies the conditions of this theorem.
It follows from the Formula (7) that the application of the FGBFI method allows one to approach the vertical asymptote arbitrarily close without jumping over it, since at Υ 0 → ∞, i.e., the Formula (7) give a guaranteed region of convergence of the power series, and for larger Υ 0 the integration step is smaller. Thus, from step to step applying Algorithm 1, we will come as arbitrarily close to time t as localizing the moment of explosion. Note also that the application of the FGBFI method for the numerical integration of the equation (3) gives an approximate value of the number π localizing the asymptote (4).

The Piecewise Smooth Systems with Quadratic-Type Nonlinearities. The Systems with Hysteresis
In practice, the problem often arises of studying piecewise smooth systems, for example, the Chua system [23]: where ρ, θ, χ, m 0 and m 1 are system parameters, the function The authors of the article [23] use the method of a describing function to localize hidden attractors in the system (22), and the values of the initial conditions are also given. Therefore, the task of developing the high-precision numerical methods for constructing approximate solutions to systems of the form (22) is important. Note that the described idea can also be transferred to models with hysteresis which were considered in [24][25][26].
We can extend the FGBFI method to piecewise smooth systems with quadratic nonlinearities. Let us describe the idea of the method using the system (22) as an example.
From the form of the function φ(x 1 ), the phase space is divided into three parts, where the right side of the (22) system is smooth. Similarly, the phase space can be divided into regions of smoothness for systems with hysteresis. For example, for systems with a non-smooth potential [25]. The joint point coordinates are By the initial condition for the coordinate x 1 , we determine in which part of the phase space the initial point is located. The modules in the Formula (23) expand accordingly. Next, we take a forward step in time ∆t according by Algorithm 1. In this case, it is necessary to remember the obtained polynomialsx 1 (t),x 2 (t) andx 3 (t) approximating the corresponding phase coordinates x 1 (t), x 2 (t) and x 3 (t) on the time interval [0, ∆t]. Ifx 1 (∆t) does not belong to the current part of the phase space, then it is necessary to find the time t pl with a high accuracy when the trajectory intersects one of the planes (24). For this, for example, one of the equationsx 1 (t pl ) = ±1 is numerically solved by the secant method.
Further, the vector of initial conditions is taken equal to and for the new part of the phase space, the modules are again expanded in the Formula (23). The advantage of the described algorithm is that we find the moments of intersection of the trajectory of the system (22) with the planes (24) with a high accuracy, slightly accumulating calculation errors.

Conclusions
In this article, the high-precision FGBFI method for constructing approximations to unstable solutions of dynamic systems of the form (1) was described, and a link to the developed program software was also provided. The arbitraryprecision numbers were represented by mpreal data type of the MPFR C++ library with overloaded arithmetic operations, as well as friendly mathematical functions.
For systems of explosive type, a theorem is proved that allows one to establish the existence of asymptotes in the system (1), and the advantage of the numerical FGBFI method in this case is also shown.
The FGBFI method has also been extended to piecewise smooth systems with quadratic-type nonlinearities using the Chua system as an example, as well as systems with hysteresis.