{'fraction_non_alphanumeric': 0.043835829001474126, 'fraction_numerical': 0.013965396850027154, 'mean_word_length': 4.409148132605959, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimation of Lyapunov exponents for quasi-stable attractors of dynamical systems with time delay

Lyapunov exponents play an important role among the invariant characteristics of dynamical systems. Analysis of the spectrum of Lyapunov exponents is widely used to study the complex dynamics in systems of ordinary differential equations and models that can be reduced to maps. In this paper we consider the question of numerical evaluation of Lyapunov exponents for delay differential equations. We used a new algorithm with FFT, tested it on Hutchinson equation with known Lyapunov exponents spectrum and compared with results of the old algorithm. The so-called quasi-stable behavior is typical for a number of gene networks models and neuronal associations that have been studied recently. The phenomenon of cycle (k-dimensional torus) quasistability in dynamical system is characterized by the fact that some of its multipliers are asymptotically close to the unit circle, and the remaining multipliers are modulo less than one (with the exception of a simple root 1 (k such roots)). In some cases, it is possible to prove the existence and give an asymptotic estimate of the multipliers of considered system by using the large parameter methods. However, if these methods are not applicable, it is necessary to obtain a tool for the numerical estimation of multipliers. Such a tool is provided by algorithms for Lyapunov exponents estimation. Given that the equations with delay are often applied in the models of neural and gene networks, the algorithm for Lyapunov exponents estimation for such systems will be in demand.


Introduction
Lyapunov exponents and Lyapunov dimension play an important role among the invariant characteristics of dynamical systems attractors. The definition of Lyapunov exponents can be found in the book [1]. For a linear system of n equations written in vector forṁ where x ∈ R n and A(t) is n × n matrix, Lyapunov exponents are determined by formula λ(x) = lim t→∞ t −1 ln |x(t)|. Analysis of Lyapunov exponents spectrum is widely used to study complex dynamics in systems of ordinary differential equations and models that can be reduced to maps. Cases where they can be found analytically are quite rare. The Benettin method is commonly used to calculate the largest exponent [2]. Further development of this method was obtained in [3]. Authors added renormalization of the initial conditions by the Gram-Schmidt algorithm [4], which allowed to calculate the spectrum of Lyapunov exponents. As follows from the Oseledets theorem [5] in the finite-dimensional case, the linearized on attractor system of  (1) is always Lyapunov proper and thus the upper limit may be replaced by a normal limit, that allows us to effectively compute Lyapunov exponents [6]. This theorem cannot be proved in the cases of equations with delay and boundary value problems. Therefore, during the development of Lyapunov exponents estimation algorithms it is important to have a model equation with delay, which spectrum of Lyapunov exponents can be computed in some other way. It allows us to test the algorithm and make sure that it is correct. Spectrum of Lyapunov exponents has been estimated in [7,8]. However, authors do not justify their algorithms and testing examples.
In this article we describe a standard algorithm extension for computing the Lyapunov exponents spectrum for systems of differential equations with delay. We show the testing results of the algorithm applied to Hutchinson equation as an example [9], and illustrate the algorithm application to Lyapunov exponents spectrum computing for some systems described in [10,11,12,13].
In the strict sense, the estimated invariant characteristics of the maps are not Lyapunov exponents. But for the case of Hutchinson equation equilibrium state, as will be shown below, we obtained results that can illustrate the "proximity" of the estimated characteristics and Lyapunov exponents. It makes it possible to suppose that our assumptions hold in other cases.

Description of Lyapunov exponents evaluation algorithm for delay differential equations
Let us describe an algorithm for obtaining the first K Lyapunov exponents for systems of differential equations with time delay of the following form: where ∀t x(t) ∈ R n , n is a system dimension, h i ∈ R n (i = 1, . . . , s), h 1 > h 2 > . . . > h s > 0. We take the space of continuous functions on the interval [−h 1 , 0] in R n , namely C([−h 1 , 0]; R n ), as the phase space. As a numerical method for solving the system (2) with initial conditions we choose the fifth-order Dormand-Prince method (DOPRI54) [14] with variable time step. We solve the system (2) with the corresponding initial condition (3) by the chosen method until the time Θ when the solution is sufficiently close to attractor. Next we obtain the function x 0 (t) ∈ R n on the interval t ∈ [Θ − h 1 , Θ], which will subsequently become the new initial condition. Using it we get the solution x * (t) ∈ R n , on which we will estimate Lyapunov exponents.
We add to the system (2) with the initial condition x 0 (t) the following K identical systemṡ where j = 1, . . . , K. They represent the systems (2) linearized on solution x * (t).
For each system (4) we use the initial conditions in the form of orthonormal impulse functions, for example: Here M is a number of solution partition points on the interval with length equals to h 1 , K ≤ nM . Also we can use a system of orthonormal trigonometric functions: Solving together the system (2) with the initial condition x(ϕ) = x 0 (ϕ) and the system (4) with the initial conditions (5) or (6) on the interval t ∈ [Θ, Θ + T ], T ≥ h 1 , we obtain for linearized systems u 1 j (t) ∈ R n (j = 1, . . . , K). Considering that u 1 j usually behave exponentially, it is necessary to renormalize them from time to time. Note that the problem is the solution's unlimited growth or vanishing. Time of renormalization T can be chosen constant or variable [15].
Further, on the interval t ∈ [Θ + T − h 1 , Θ + T ] we discretize the obtained solutions into M equal parts and use one of the algorithms described below.
Impulse functions method. We do orthonormalization of u 1 j (t), j = 1, . . . , K by Gram-Schmidt method [4]. At the same time, after the orthogonalization procedure and before the normalization, we calculate and keep the values Then we re-solve the systems (2), (4) with obtained orthonormal solutions used as initial conditions. Trigonometric functions method. Before the orthogonalization procedure we apply the Fast Fourier Transform to u 1 j , obtaining the complex coefficients c j = c j,i , i = 1, . . . , M/2 + 1. Then we split coefficients into real and imaginary parts (at cosine and sine) d j = d j,i , i = 1, . . . , M +2. We apply the Gram-Schmidt method to d j . At the same time, after the orthogonalization procedure and before the normalization, we compute and keep the values Finally, we use the inverse Fast Fourier Transform to obtain the initial conditions on the time delay interval. Thus, we obtain an approximation of Lyapunov exponents: Let us proceed to the results of the testing algorithms.

Algorithms test results
Computational experiments were carried out for the Hutchinson equation [9]: Non-zero solutions of the equation (10) are asymptotically stable for r ∈ [0, π/2] (monotonically stable for r ∈ [0, e −1 ]) and tend to 1 in an oscillatory manner for r ∈ [e −1 , π/2]. In this case, the Lyapunov exponents of (10) are equal to real parts of the roots of characteristic quasi-polynomial P (λ) ≡ λ + r exp(−λ), λ = τ + iω. To calculate them, a system of algebraic equations is used: τ + re −τ cos ω = 0, ω − re −τ sin ω = 0.  Table 1. The first 10 computed Lyapunov exponents (λ i -by impulse functions method, and λ i -by trigonometric functions method) for the equation (10) at r = 1.0 and relative difference σ i between them and reference values. The rounded components τ i of this system's solution at r = 1.0 are presented in the second column of Table 1. We will call them reference values.
Linearized equations for (10) have the form: For all experiments the following parameters were used: number of computed Lyapunov exponents K = 10, time to approach the attractor Θ = 150, renorming time T = 4, number of Lyapunov exponents recomputations L = 5000, as initial conditions we used the system (5) in the case of impulse functions method and the system (6) in the case of trigonometric functions method.
As shown in Table 1, the accuracy of exponents computing depends on the chosen partition points count. With an increase of M by 20 times from 100 to 2000, the relative difference decreased by an order. When the partition points number is equal to the number of computed Lyapunov exponents, the achieved accuracy is low.

Lyapunov exponents computing in case of dynamical systems with a quasi-stable attractor
The study of dynamic systems from neurobiological and physical applications leads to singularly perturbed systems, which are characterized by so-called relaxation oscillations. Recently, among the problems with relaxation dynamics, we can distinguish a subclass of equations, which stable modes have multipliers that are modulo less than 1, but still close (sometimes asymptotically) to the unit circle on the complex plane. Attractors having this property can be called quasistable [10]. This fact make numerical analysis of such problems difficult, since the approach to the stable mode is very slow and, based on numerical analysis, it is not possible to find out the form of stable modes. At the same time, even a long computation can lead to incorrect conclusions about the stability of some part of the transition process.
In such situation, computing the Lyapunov exponents along the studying problem solutions can help. The presence of any number of values close to zero in the Lyapunov exponents spectrum indicates the problem's quasi-stability. In many cases the model's quasi-stability allows to make conclusions about adequacy or inadequacy of this problem to the model.
As an example of this situation, we consider a mathematical model of synaptically connected  Table 2. The first 4 calculated Lyapunov exponents of different modes of the system (13) at a = 1.5, b = 7, c = −4.5, h = 7, λ = 100 after 4000 steps of renorming rounded to the 4th digit.
This model was proposed in [12] (see also [13]). There u 1 (t), u 2 (t) > 0 are normalized membrane potentials of neurons. Each neuron is modeled by singularly perturbed differential-difference equation with delayu = λf (u(t − 1))u, the parameter λ 1 characterizes the rate of electrical processes, the term bg u j−1 (t−h) ln(u * /u j )u j models the synaptic interaction. Communication is assumed to be a threshold with a time lag h > 1.
Examples of functions that satisfy these properties are f (u) = (a − au)(a + u) −1 , g(u) = u(1 + u) −1 . Solutions of the system (12) at large λ are intense peaked oscillations, which are alternated by area where the solution is close to zero (see [10,11,12,13]). For simplicity of the problem's numerical analysis, we perform an exponential substitution u j = exp(λx j ), j = 1, 2, and the system (12) takes the forṁ For the system (13) at a = 1.5, b = 7, c = −4.5, h = 7, λ = 100 we have found at least 6 (without symmetry) different coexisting stable cycles. However, the transition process leading to the first five of them takes a lot of time. The use of the methods described in the previous paragraph allows us to understand the reasons for this phenomenon. As it turned out, in addition to one zero exponent, which is determined by the autonomy of the system (13), the first three self-symmetric cycles have two additional exponents close to zero (Fig. 1 a), b), c)) and nonsymmetric cycles have one additional exponent close to zero (Fig. 1 d), e)). The case presented in Fig. 1 f) stands alone, as it corresponds to only one zero exponent and the transition to this cycle from its attraction area is incomparably faster than in previous cases.

Conclusions
Quasi-stable structures are detected in mathematical models with relaxation behavior of various nature. Often enough, the models in such problems are equations or systems of equations with time delay. In this regard, the tools of numerical analysis of such systems are very important. Calculating the spectrum of Lyapunov exponents makes it possible to understand that the studied attractor is quasi-stable and we have to do sufficiently long computing and use special initial conditions to get reliable results. Note that in the case when dynamical system attractors are cycles or tori, there is no need in such statements as Oseledets theorem. The significant problems with the described algorithms applicability can arise only in the case of chaotic behaviour of corresponding system solutions.