Features of the Computational Implementation of the Algorithm for Estimating the Lyapunov Exponents of Systems with Delay

The computational implementation of the algorithm for estimating the spectrum of Lyapunov exponents for systems of differential equations with delayed arguments is considered. Taking into account that for such systems, as well as for boundary value problems, it is not possible to prove the well-known Oseledets theorem, which makes it possible to efficiently calculate the required values, we only have to talk about estimates of characteristic exponents, in a sense close to Lyapunov ones. In this paper, we propose two methods for processing solutions of systems linearized on an attractor, one of which is based on a set of impulse functions; the other one is based on a set of trigonometric functions. The flexibility of application of the indicated algorithms is demonstrated in the case of quasi-stable structures, when several Lyapunov exponents are close to zero. The developed methods are tested on the logistic equation with delay. The results illustrate the “closeness” of the estimated characteristics and Lyapunov exponents.


INTRODUCTION
An extension of the standard algorithm [1] for calculating the first few Lyapunov exponents for systems of differential equations with delayed arguments is considered.
Note that Lyapunov exponents for systems with a delay may not quite be correctly estimated numerically. The fact is that for finite-dimensional systems there is a well-known Oseledets theorem [2], which states that a system linearized on a stable solution is always Lyapunov correct. This makes it possible to replace the upper limit with the usual one in determining the Lyapunov exponents [3] and numerically estimate these values. In the case of equations with delayed arguments and boundary value problems, such a theorem cannot be proven. Therefore, when developing algorithms for calculating Lyapunov exponents, it is important to have a model equation with delay for which the spectrum can be calculated in some other way. The presence of such a task allows you to test the developed algorithm and make sure that it works. In the articles [4][5][6], the spectra of Lyapunov exponents are calculated, but the authors do not provide a justification for the proposed algorithm, as well as a testing example, in contrast to work [7], which served as the basis for the methods presented in this article. We will mention separately [8], one of the first works on this topic. Note that all the results described in this article are experimental.

DESCRIPTION OF THE ALGORITHM
Let us describe the process of obtaining estimates for the first Lyapunov exponents in the case of systems of differential equations with delayed arguments of the following form: (1) where for , is the dimension of the system, ( ), and .
As a phase space, we take the space of functions in continuous on the segment , namely .
For the numerical solution of system (1) with initial conditions (2) we will use the fifth-order Dormand-Prince method (DOPRI54) with a variable step length [9].
Thus, we will solve system (1) with the corresponding initial condition (2) by the chosen method until the moment of time , sufficient for the solution trajectory to approach the attractor under study. In this case, on the interval , we obtain the function that will become a new initial condition of system (1).
Let us supplement the system of equations (1) with the following identical systems: where , is the solution of the system of equations (1) with the initial condition at . They represent systems of equations (1) linearized on the solution . In what follows, we will denote for which we introduce the norm (4) For each system of equations (3), we use the initial conditions in the form of orthonormal impulse functions, for example: (5) or in the form of orthonormal trigonometric functions: (6) Jointly solving system (1) with the initial condition and systems of equations (3) with initial conditions (5) or (6) on the interval , where , we obtain the solution , for each of the linearized systems.
Given that the solutions behave exponentially, it is necessary to renormalize them at regular time intervals. Note that both the unlimited growth of solutions and their tendency to zero are a problem. Thus, over the interval , we average the calculated solutions of the linearized systems within each of equal time intervals of length , as a result of which we obtain piecewise continuous functions , respectively, which we use in one of the methods described below. The method of impulse functions.
• We apply the Gram-Schmidt method [10] to . • In this case, after the procedure of orthogonalization of each function, but before its renormalization, we calculate the values , where is the norm defined in (4), is the orthogonalized system of the functions . • We continue to solve systems (1) and (3), while using the resulting orthonormal system of functions as initial conditions for linearized systems.
To apply this method, it is necessary to choose a system of orthonormal impulse functions (5) as the initial conditions for linearized systems.
The method of trigonometric functions.
• First, we transform the functions into a system of vectors according to the following rule: , , , . • To the resulting system of vectors, we apply the discrete Fourier transform [11], as a result of which we obtain complex-valued vectors .
• We divide the vectors into pairs of real-valued vectors consisting of real and imaginary parts .
• We apply the Gram-Schmidt method to the vector system . • In this case, after the procedure of orthogonalization of each vector, but before its renormalization, we calculate the values where is the orthogonalized system of the vectors , . • The obtained orthonormal real-valued vectors are converted back to complex-valued ones, to which we apply the inverse discrete Fourier transform, thus obtaining the vectors .
• We build a system of functions according to the rule: = at , , , .
• We continue to solve systems (1), (3), while using the system of the functions as the initial conditions for linearized systems.
To apply this method, it is necessary to choose a system of orthonormal trigonometric functions (6) as the initial conditions for linearized systems.
We repeat the described process at time intervals , , as a result of which the corresponding solutions are processed by the algorithm. The estimation of Lyapunov exponents in this case is calculated by the formula (7) Note that the renormalization time can be selected in two different ways: at regular intervals or dynamically [12]. In the latter case, at each step of the algorithm, you will have to store not only , but also . Then the formula for estimating Lyapunov exponents will have the form The described structure of the algorithm makes it possible to start calculations even before the solution reaches the attractor under study, especially when the process of approaching the solution to the attractor itself can be associated with great computational difficulties, as shown in [13] for the case of several Lyapunov exponents close to zero for quasi-stable structures. In this case, the first steps corresponding to the time intervals in which the solution has not yet approached the attractor by a sufficient distance should be discarded and not taken into account in any way in the sum from formula (7). It is also recommended to discard the first few steps in the general case, since the process of forming a new orthonormal basis of linearized systems can significantly affect the spectrum of exponents with not too many steps taken into account (see a similar technique in [14]). Now, let us move on to the results of testing the described methods.
(1) ( )  [15], which has the following form: (8) It is shown in [16] that nonzero solutions of Eq. (8) are asymptotically stable at , moreover, the solution tends to one at monotonically, and at oscillating. In addition, the   unit solution has global stability at r ≤ 37/24 [16][17][18], and in [19,20] an algorithm is given that allows to improve this estimate. In the case of the unit equilibrium state at r < π/2, the Lyapunov exponents for the Hutchinson equation coincide with the real parts of the roots of the characteristic quasi-polynomial , . To calculate them, a system of algebraic equations is used:   The rounded components of the solution τ i of this system at different values of the parameter r are presented in the second column of Tables 1-3. We will call them reference values.
Equation (8)    . For convenience, the calculated estimates of the Lyapunov exponents, the absolute and relative differences, have been rounded to the fourth decimal place.
As can be seen from Tables 1-3, the accuracy of calculating the exponents depends on the size of the selected partition. When increasing by 10 times from 100 to 1000, the relative error decreased by an order of magnitude. In the extreme case, when the number of partition points is equal to the number of calculated Lyapunov exponents, the achieved accuracy is low. Note that in some cases, for example, for the modes of the model from [13], the method of trigonometric functions may be inapplicable due to the loss of information at the stage of applying the discrete Fourier transform.

CONCLUSIONS
Thus, the conducted numerical experiments show that, when choosing a sufficient number of partition points , the estimated characteristics can be qualitatively close to the Lyapunov exponents. In the case of a small dimension of system (1), the choice of the method affects the speed of calculations, depending on the nature of the solution under study. In particular, if the solution is relatively smoothed, then the method of trigonometric functions works faster than the method of impulse functions. If the solution contains a large number of sections with sufficiently sharp peaks, then the method of impulse functions is preferable in terms of execution speed. In addition, with an increase in the number of calculated Lyapunov exponents and the number of partition points , as well as the dimension of the system under study, the use of multiprocessor parallel systems becomes effective.