A new way to compute the Lyapunov characteristic exponents for non-smooth and discontinues dynamical systems

One of the most important problems of nonlinear dynamics is related to the development of methods concerning the identification of the dynamical modes of the corresponding systems. The classical method is related to the calculation of the Lyapunov characteristic exponents (LCEs). Usually, to implement the classical algorithms for the LCEs calculation, the smoothness of the right-hand sides of the corresponding equations is required. In this work, we propose a new algorithm for the LCEs computation in systems with strong nonlinearities (these nonlinearities cannot be linearized) including the hysteresis. This algorithm uses the values of the Jacobi matrix in the vicinity of singularities of the right-hand sides of the corresponding equations. The proposed modification of the algorithm is also can be used for systems containing such design hysteresis nonlinearity as the Preisach operator. (Thus, this modification can be used for all members of the hysteresis family.) Moreover, the proposed algorithm can be successfully applied to the well-known chaotic systems with smooth nonlinearities. Examples of dynamical systems with hysteresis nonlinearities, such as the Lorentz system with hysteresis friction and the van der Pol oscillator with hysteresis block, are considered. These examples illustrate the efficiency of the proposed algorithm.


Chaos and methods to identify it
It is well-known that dynamical systems can demonstrate not only regular evolution, but also a stochastic or chaotic mode. The main feature is the fact that chaotic dynamics take place in the absence of random factors and are completely determined by the internal properties of the system and initial conditions. One of the first papers determined to chaos in the continuoustime system was published in 1963 by Lorenz [37]. This paper was a start point for such a wide field of modern science as dynamical chaos. To these days, the list of publications devoted to the dynamics of various chaotic systems include dozens of thousands (see, for instance, review papers [1,13,32], etc.).
Diagnostics of the chaotic behavior of systems is one of the most important stages in the study of dynamics. At this point, the variety of methods for identifying dynamical chaos is quite wide: spectral analysis, method of LCEs, 0-1 test, and method of the fast Lyapunov indicator (FLI). In this list, the LCEs method seems much powerful compared to other ones. As is commonly known, the classical LCEs method requires the smoothness of the right-hand sides of the corresponding equations. At the same time, models of important from a practical point of view technical systems contain non-smooth and discontinuous terms including terms of a hysteresis nature [9,15,46]. Thus, diagnostics of the dynamical regimes of such systems require a generalization of standard methods based on LCEs to include non-smooth and operator nonlinearities. Various approximations (such as functional approximations) of the right-hand sides with non-smoothness leads to unpredictable errors in the calculation of LCEs. Things are even worse in systems with a design model of hysteresis. The possibility of an adequate approximation employing functional nonlinearities in such systems is omitted. Therefore, it seems important to generalize the methods of diagnostics of the dynamical modes (initially based on the calculation of the spectrum of LCEs) to the study of systems with non-smooth right-hand sides, as well as with hysteresis nonlinearities.
The list of publications on systems with hysteresis nonlinearities encountered in scientific and technical problems numbers thousands of works. We mention only one review devoted to the Bouc-Wen model [24], as well as recent works [21,30,35,36,39,44,49]. Work [29] proposed a new method for obtaining random sequences using hysteresis converters under thermal fluctuations.
The Fourier analysis method is considered to be one of the simplest methods for detecting chaos in dynamical systems [22,40]. The power spectrum of signal s(t) is defined as the square of the modulus of the Fourier spectrum F[s(t)], which is equivalent to calculating the spectrum of the autocorrelation function s(t)s(t +t ) . To carry out numerical calculations, the continuous signal is presented as a discrete sequence (a transition to a time series with elements s j = s(t j ), j ∈ N), then numerical implementations of the algorithm are used (usually the fast Fourier transform (FFT)). This method makes it easy to distinguish quasiperiodic dynamics, which corresponds to a discrete spectrum from chaotic, whose spectrum is continuous. However, this method does not allow to reveal the nature of chaotic behavior. This behavior may be also caused by the presence of purely random components in the studied signal. In other words, the Fourier analysis method can only be used for a qualitative analysis of the behavior of dynamical systems, because of insufficient specificity in the sense of false-positive results from noisy time series.
In recent times, a method for detecting a chaotic regime using the 0-1 test has been increasingly encountered in the literature. The absence of requirements for information about the system parameters is its distinctive feature. To implement this method, we only need a set of values for the time series of the analyzed process [12,20]. Let we have sequence ϕ( j) for the discrete parameter j = 1, 2, . . . , N . Let us construct sequences p(n) and q(n) following rules: where c ∈ (0, 2π) does not depend on parameter n. Next, the time mean square deviation of sequence ϕ( j) is determined: For n ∈ N, the exponential growth rate is: The value of K is can be used to estimate the presence of chaotic behavior: K = 1 for chaotic dynamics, and K = 0 for regular evolution of ϕ( j). But again, as for spectral analysis, the 0-1 test does not allow to distinguish between deterministic chaos and stochastic sequences. As was mentioned above, the method based on the calculation of LCEs is the most powerful tool for diagnosing deterministic chaos. It is known that trajectories of deterministic systems with chaotic behavior diverge locally with exponential speed. LCEs allow quantifying this speed [6,19,23,26,43]. Indeed, if we consider two initial points in the phase space located at distance ε 0 , for relatively small values of time t, the distance between trajectories grows following the exponential law: where λ is a positive constant called the largest LCE.
A dynamical system with n degrees of freedom has exactly n such exponents that correspond to changes in the distance between trajectories along each of the coordinate axis. Here, just a brief description of LCEs and more detailed consideration are given in Sect. 1.2. Together with the considered criteria for chaos in dynamical systems, there exists a relatively simple criterion which is called as the fast Lyapunov indicator (LPI) [17,18,27,34]. Consider the system of ordinary differential equations: where x = (x 1 , x 2 , . . . , x n ) ∈ R n is a vector of variables, and F is some vector field. The corresponding initial conditions will be denoted as x(0) = x 0 . Introduce an auxiliary vector v defined by the equation: By definition, the LPI indicator is In some cases, for instance, when the LPI is oscillatory, it is convenient to use the global maximum of the LPI in the time interval [0, T ] instead of the instant value of ln |v(t)| [4]: In [34], it was shown that dependence of the LPI on time is logarithmic for regular trajectories, while for chaotic trajectories, the dependence is linear. Thus, this fact makes it possible to "recognize" chaos in the system.

Spectra of LCEs
The LCEs are defined as the eigenvalues of the matrix of infinitesimal perturbation of the trajectory in the phase space. They are usually considered as a measure of the system's sensitivity to small perturbations of the initial conditions. The largest LCE is often used to distinguish between chaotic and regular modes [16,51]. Let the dynamical system be described by the following vector differential equation: where x ∈ R n , F : R n → R n . Consider a trajectory in the configuration space starting at point x 0 : The Jacobi matrix is the derivative of the particular solution ϕ(x 0 ) with respect to components x 0 : It shows the change in the trajectory as a result of an infinitesimal perturbation of the initial conditions δx 0 : The norm of vector δx(t) in R n is written through initial conditions as follows: We denote the eigenvalues of matrix Φ T (x 0 )Φ(x 0 ) as m i , where i = 1, . . . , n, and its eigenvectors as v i , i = 1, . . . , n. It can be shown [51] that scalar product (v i , δx(t)), which is defined in the standard way in R n , is proportional to (v i , δx 0 ): The spectrum of LCEs is defined as an ordered set (λ 1 , λ 2 , . . . , λ N ) where components are: Definition (16) follows the asymptotic estimate: where notation Θ(g(n)) = { f (n) : ∃c ∈ (0, ∞), n 0 ∈ N is used such that for all n n 0 , c 1 g(n) f (n) c 2 g(n)} is executed.

Algorithms to compute LCEs
There are several algorithms for the numerical estimation of the LCEs. We briefly consider some realizations of such algorithms.

Benettin method
Consider some point of the phase trajectory x 0 of a three-dimensional dynamical system at time instant t 0 . Choose ε > 0 and an arbitrary pointx 0 so that the equality |x 0 −x 0 | = ε is satisfied. Let us trace evolution in time of points x 0 andx 0 on a short time interval [t 0 , t 0 +T ], while the basic and variational systems will be transformed into new points x 1 andx 1 , respectively. (The variational system is obtained by linearizing the basic in a vicinity of the trajectory of the basic system.) To estimate the largest LCE use the following expression: Time interval T is chosen so that the amplitude of the perturbation is less than the linear dimension of the inhomogeneity in the phase space. We normalize the displacement vector Δx 1 = ε(Δx 1 − Δx 0 )/|Δx 1 | and calculate new point x 1 = x 1 + Δx 1 . Next, we repeat the above procedure having previously replaced x 0 and x 0 by x 1 andx 1 , respectively.
After repeating this procedure M times, λ is determined as the expectation of estimatesλ l at each step: The numerical implementation of this procedure is performed in [10,11].
Note that time interval T should not be taken too large otherwise all perturbed trajectories will be deflected in the direction of the vector corresponding to the maximum exponent.

Wolf method
This method is suitable for calculating nonnegative elements of the Lyapunov spectrum. The main idea of the Wolf method is that the largest LCE is calculated using only the time series of the system [53].
Consider time series x(t), t ∈ {1, 2, . . . , N }, which is a discretization of chaotic process obtained at equal time intervals τ . The mutual information method allows to determine the value of the time interval τ and the false neighbors' method gives dimension n of the embedded space. As a result, we obtain a set of points in space R n : where i = ((n − 1)τ + 1), . . . , N . Take some point from series (20) and denote it by x 0 . In the series, we can find pointx 0 such that |x 0 − x 0 | = ε 1 < ε. Let us trace the evolution of these points in time until the distance between points reaches value ε max . The new points will be denoted by x 1 and x 1 , the distance is equal to ε 0 , and the corresponding time interval will be denoted by T 1 . Then, consider the sequence again and find pointx 1 : By repeating this procedure M times, the largest LCE is estimated as: This method can be applied also in the case when exact equations of the system are unknown, and we have only the time series. This is the specific feature of the Wolf method.

Rosenstein method
This method calculates function: where x j is the initial point of the phase trajectory, There is a relationship between d j and the LCE: The largest LCE is calculated by estimating the slope of the linear part of this function.

Kantz method
The algorithm proposed by Kantz [25] calculates the largest LCE following expression: where x t is an element of the time series, U t is a neighborhood of x t , point x i belongs to this neighborhood; τ is the relative time multiplied by the sample rate; T is the sample size, S(τ ) is a stretch coefficient in the linear growth area. The angular slope of the straight line is equal to the largest LCE, i.e., e λτ ∼ e S(τ ) . Note that the list of known numerical methods for estimating LCEs presented in this section is not limited. There are also known methods based on the use of neural networks, wavelet analysis, etc. Note that the above-mentioned methods are applicable only for systems with smooth right-hand sides. Thus, when analyzing systems with non-smooth (as well as discontinuous) nonlinearities in the overwhelming majority of publications, these nonlinearities were approximated by smooth functions. Moreover, the question on the correctness of this operation remains open. Special attention should be paid to the work by Balcerzak et al. [7], where a method for calculating LCEs for systems with nondifferentiable right-hand sides is proposed. Also, we note [42] where the author introduced the transition conditions at the instant discontinuity. However, for systems with operator nonlinearities (including hysteresis nonlinearities formalized within the design approach such as the Preisach and Ishlinskii operators), correct methods for calculating LCEs have not yet been reported. Note that in works [7,42], discontinuities such as the Coulomb friction were considered and did not report about results for other strong nonlinearities including the operator ones. Thus, new instrumental methods for computing LCEs in such systems seem important and promising.

Non-ideal relay
The classical approach to the description of hysteresis nonlinearities developed by Krasnoselsky and Pokrovsky [31] considers hysteresis operators as converters defined on the space of continuous functions, the dynamics of which is described by the relations: "input-state" and "state-output." Denote by R[α, β, x 0 , t 0 ] a hysteresis converter corresponding to a non-ideal relay with threshold numbers α and β where x 0 ∈ {−1, 1} is an initial state of the converter, t 0 is an initial time instant. The state space of the non-ideal relay is the two-element set {−1, 1}. The input of the system is function u(t), which is continuous for t t 0 , and the output is the step function x(t), defined by the operator relation: In this case, initial state x 0 of the converter satisfies the condition: As a result, the output has a value on [t 1 , . We will say that the relay is on if the output is equal to one and that the relay is off otherwise. A detailed description of the non-ideal relay converter is given in [31], as well as its properties.
The relationship between the input and output of the R[α, β, x 0 , t 0 ] operator is schematically shown in Fig. 1.
The Preisach operator is a continuous analog of the family of non-ideal relays connected in parallel. The state space of this converter consists of the pairs {u(t), z(α, β, t)} where u(t) is the input value at time instant t and z(α, β, t) is the function of a subset on the half-plane α < β taking values ±1 [31,47]. The inputoutput relations of the Preisach operator are determined where z 0 = z(α, β, t 0 ), and the output-state (see Fig. 2) As a rule, in applied problems, the carrier of the measure of the Preisach operator is bounded. Thus, the state space consists of functions on bounded sets of the corresponding plane.

The derivative of the Preisach operator with respect to input at time instant t
By the weak differential, or the Gâteaux differential of the mapping, F(u, ϕ): where ϕ(t) ∈ C 0 (−L , L). In this case, convergence is meant in the sense of convergence by norm [5,28]. The weak derivative (the Gâteaux derivative) of an operator is the linear with respect to ϕ(t) part of the weak differential [14]: Note if equality (30) is satisfied, the derivative of operator F(u, ϕ) is called as the "directional derivative with respect to ϕ".
In the general case, the Preisach operator does not even have a weak Gâteaux derivative. However, there are inputs and corresponding directions (in the functional space) for which an analog of the Gâteaux derivative takes place. The corresponding constructions are described below.
In the simplest situation when the carrier of the measure of the Preisach operator is placed in triangle {(α, β) : α ∈ [−L , L], β ∈ [α, L]}, the derivatives with respect to the parameter introduced above admit analytical representation.
The state of operator (28) and τ ∈ [0, T ] that the output of the operator with initial state z 0 (α, β) ≡ 0 is defined as: In other words, there is such a piece-wise monotone input that transfers the operator with initially offed relays to state {u(τ ), z(α, β, τ )}.
if u(t) increases monotonically in the left neighborhood of t = τ ; if u(t) decreases monotonically in the left neighborhood of t = τ and u (2k L], and the following identities take place: Proof We fix an arbitrary time instant τ ∈ (t 0 , T ) and consider the 2k-th monotonicity segment where u(t) increases monotonically and for which the inequalities are true u (2k−1) < u(τ ) < u (2k) (see the left panel in Fig. 3).
From the geometric meaning of the double integral, implies that as μ increases, integral I has a small increment corresponding to the area of the trapezoid highlighted by a solid fill and equal to For μ → +0, the value of the derivative is The right panel in Fig. 3 shows the case when input u(τ ) + μϕ(t) reaches a value equal to the local extremum u (2k−1) , the sides trapezoids as functions of μϕ(t) have a discontinuity, and the area becomes equal to Herewith, the following operator relations take place and so on. Thus, the derivative with respect to parameter (the input at time instant t) is obtained for increasing inputs. Further, after reaching the value of u (2k) on the 2k + 1-th section of monotonicity where function u(t) monotonically decreases, similar reasoning leads to relations of the form: if u(t) decreases monotonically in the left neighborhood of t = τ and u (2k) < u(t) < u (2k+1) . The theorem on the differentiability of the Preisach operator with respect to parameter u(t) is proved.
Note It is easy to see that the following relations take place: u (2k+1) > u (2k) and −L u (k) L for all 1 k M.
Example 1 For the Preisach operator in a state characterized by the distribution function of the switched on elementary relays z(α, β) ≡ 0 and for monotonic function u(t) ≡ (t − t 0 ) − L, belonging to set C 0 (−L , L) at t 0 t t 0 + 2L, we have the following relations: Thereby, the differentiation by parameter leads to Further, the derivative of the Preisach operator with respect to the input D u(t) P[z(α, β, t), u(t), t 0 ]u(t) means the right-side derivative at a monotonically increasing input and the left-sided derivative at a decreasing one.  As it can be easily seen, derivative D u(t) P[z(α, β, τ ), u(τ ), t 0 ]u(τ ) is a piece-wise linear function that has discontinuities of the first kind at points of the local extremum of u(t), as well as at those points where there is an intermittent in the output of the converter due to the geometry of the on and off elementary hysterons. Particularly, when the input increases from u(129) = 0.4 to u(132) = 0.5, a set of three elementary hysterons is switched on (in the fourth row from the top, see Fig. 5), which leads to a discontinuity of the derivative at point t = 132.

Phenomenological description of the hysteresis
One of the most commonly used phenomenological hysteresis models is the Bouc-Wen model [24].
Consider a one-dimensional mechanical system in which the restoring force Φ BW (x, t) is represented as the sum of two terms: where x(t) is the displacement of the body along the coordinate axis, αkx(t) is an elastic restoring force, k is the stiffness coefficient, (1−α)Dkz(t) is the hysteresis term, α is the relative stiffness, i.e., the ratio of deformation outside the elastic threshold to the deformation under the elastic limit, D is the plastic displacement deformation. The auxiliary dimensionless function z(t) satisfies the first-order differential equation: where the right-hand side contains the time derivativė constructed in this way is called the Bouc-Wen hysteresis model [52]. Parameters A, β and γ are dimensionless and determine the shape of the hysteresis loop, parameter n determines the smoothness of the transition from elastic to plastic deformation. (It is usually supposed to be integer.) Based on the physical meaning of Eq. (44), the following conditions for real parameters D, k, α are usually take place D > 0, k > 0 and 0 < α < 1. As a result, these seven parameters describe the hysteresis properties of the system within this model. Such a multiparameter model allow to describe a wide class of phenomena taking into account a hysteresis relationship between their parameters. First, we calculate the next value x ( j) . Then, we compute the Jacobi matrix Φ Δt j (x ( j−1) ) and use it for calculation of new perturbations y Because of significantly different displacements from the unperturbed trajectory in different spatial directions, the Gram-Schmidt orthogonalization procedure is used to preserve the accuracy of calculations: Performing the above steps for a sequence of split points t 1 , t 2 , …,t j , …,t J , where j = 1, 2, . . . , J , we have an estimation of the Lyapunov spectrum: The considered algorithm is also correct for systems whose right-hand side contains nondifferentiable functions. In this case, the derivatives in the calculation of the Jacobi matrix are approximated by finite differences for some sufficiently small (in the sense of constructing) real Δ. Calculation of matrix Φ t (x) on values x perturbed by the displacement of Δ is performed in the direction of each of the coordinate axes, i.e., along the unit vectors e 1 , …, e n . In [7,8], it was shown that if the trajectory φ Δt j (x) is differentiable at points x ( j−1) even in the presence of nondifferentiable singularities (regardless of their type), the algorithm based on (52) is correct.
To apply the considered algorithm to systems with phenomenological hysteresis (such as the Bouc-Wen, Dahl, and Duhem models) or discontinues nonlinearities, we propose the following specific modification.
In contrast to the classical Benettin algorithm where instants t 1 , t 2 , …, t J are fixed during all time interval of the LCEs calculation, in the proposed algorithm, the set of these instants is expending by instants τ d , d = 1, . . . , D, where τ d corresponds to the discontinuity or the smoothness loss points in the right side of (6). To be specific, assume that right-hand side F(x) = (F 1 (x), . . . , F n (x)) in (10) contains singularity in the first component of set F 1 (x). Suppose that the smoothness loss arise on surface S ϕ which is determined by ϕ(x) = 0. Assuming that trajectories of (6) are transversal (cross the surface under nonzero angle) with respect to S ϕ (i.e., there are no sliding modes), there is only a finite set of instants τ d that correspond to the cross of the trajectories with the non-smoothness (or discontinues) surface during the finite time interval. These instants should be added to set t 1 , t 2 , …, t J . Herewith, the Jacobi matrix calculation starts at time instants t d + 0 and stops at t d+1 − 0. Particularly, in the case when the right side of (6) has the form of the Preisach operator the corresponding smoothness loss arising at instants where the monotonicity of input x 1 (t) changes. Then, ϕ(x) = 0 is the same asẋ 1 (t) = F 1 (x) = 0.

Lorenz system with hysteresis friction
To illustrate the performance of the proposed algorithm we consider one of the most famous examples of dynamical systems demonstrating chaotic behavior. This is the Lorenz system [37]: where σ , r , b are nonnegative.
These equations are encountered in the study of dynamical systems of various physical nature. Among other factors, Eq. (53) follows from a system of partial differential equations first used to describe convection flows in the lower atmosphere represented as a model of a continuous medium that is uniformly heated from below and cooled from above. Parameters σ , r , and b have a clear physical meaning. The first of them σ is the ratio of the kinematic viscosity (momentum diffusivity) coefficient to the thermal diffusivity (Prandtl number). The normalized Rayleigh number is denoted by r , and b is the dimensionless scale of the convective cell. Parameter b is determined only by the geometry of the system, while r and σ depend not only on the geometry, but also on the hydrodynamic properties of the medium.
Physically interesting behavior of the Lorentz system manifests when r changes, which corresponds to a change in the external excitation. Note that variables x, y and z are not Cartesian coordinates, but relate to the temperature distribution (y and z are horizontal and vertical parts, respectively) and velocity (x) of points of the medium.
Optical systems (the amplitude of wave oscillations in a single-mode laser) and even mechanical systems such as a tube closed in a ring filled with a liquid, and also a water wheel lead to a system of equations similar to (53). In the problem of a water wheel, the wheel mounted on a horizontal axis. Small containers are attached to the rim of the wheel in which water pouring from above is retained. If containers allow water loss due to its flowing down, then with a sufficiently intense flow, the dynamics of the wheel rotation shows chaotic behavior. In the problem of a rotating wheel, friction can be naturally introduced in the form of an additive term in the first equation of (53): where F(x) is the friction in the wheel axle. As shown in [48], dry friction (Coulomb friction) leads to stabilization of the system at any nonzero coefficient of friction. Nevertheless, hysteresis friction models (such as the Bouc-Wen and Dahl models) permit controlled chaos regimes.
Due to the fact that the compression ratio of the phase volume is the Lorentz system is not Hamiltonian [3]. In configuration space (x, y, z), the phase trajectory of the system forms the so-called strange attractor (the Lorentz attractor).
Consider dynamical system defined as follows: where P[z(α, β, t), x(t), t 0 ] is the Preisach operator. As is well-known, for the values of parameters σ = 10, r = 28, and b = 8/3, the classical Lorentz system (parameter a is equal to zero) exhibits chaotic behavior. In this case, the smoothness loss surface is determined by −σ x +σ y +a P[z(α, β, t), x(t), t 0 ]x = 0. In general case, it is impossible to obtain the corresponding solution and determine the smoothness loss surface because an output of the Preisach operator depends on the whole prehistory. (It is determined by the dynamics of the system at all previous time instants). Thus, to compute the LCE, it seems more convenient to track the sign ofẋ and determine the corresponding instants τ d as the instants where the derivative changes sign. Figure 6 shows the Lyapunov spectrum and phase portraits for (56).
The output of the Preisach operator has the form of a characteristic hysteresis loop at relatively small values of amplitude a (see Fig. 7), in other words, the hysteresis block "operates" during the entire simulation time. However, at a 50 due to the presence of a stable focus in the phase space, a monotonic dependence x(t) is observed, which means that the input of the hysteresis operator is stabilized. Note that numerical results shown in Fig. 6 (panels d and f ) in the case of constant external excitation (F(x) ≡ const) also show stabilization.
Note also that as a increases, a transition from chaos to a deterministic regime is observed, as shown in Fig. 8 (the dependence of the largest LCE on a).
Another widely used phenomenological model of hysteresis is the Dahl model. Add to the right-hand side of the first equation of the Lorenz, an inhomogeneity corresponding to friction within the Dahl model is: where a is the amplitude of the hysteresis force, , F c are parameters of the Dahl model, sign is the standard Fig. 9 Lyapunov spectra of the Lorentz system with the Dahl term(left panels) and phase portraits of the Lorentz system with the Dahl term (right panels) at following parameters: = 3.0, F c = 2.0, a = 0, 1, 3, 5 signum function. In this case, surface S ϕ determining a set where the right side loses the smoothness has the form: σ y −σ x −a F = 0. Figure 9 shows the Lyapunov spectrum and phase portraits for (57). As it can be seen, the chaotic regime is replaced by a deterministic one when a increases. At the same time, the phase volume collapses when a increases.

System with phenomenological hysteresis
Another example of nonlinear systems is an oscillator with nonlinear damping. For certain values of the parameters, stable relaxation oscillations exist in such a system (it is the well-known van der Pol system [41,45]). This system plays special role in a wide range of problems in applied mathematics and simulations [33,38,50].
To the right side of the differential equation describing the damped van der Pol oscillator subjected to external harmonic excitation (59) Figure 10 shows the Lyapunov spectra and phase portraits of the system at ω = 1.0, μ = 5.0, λ = 1.0 for various values of parameter a . The frequency of the external harmonic excitation is chosen equal to b = 2.467, and the amplitude is equal to B = 5 (as in [2]).
The bifurcation diagram for the considered system depending on a is shown in Fig. 11 together with the corresponding Lyapunov spectra. An unexpected behavior of the system is observed: when a increases from a = 0 to a ≈ 3.95, a deterministic regime takes place, while it changes at a 3.95 to chaos. The same behavior can be seen from the Lyapunov spectra.

Conclusions
In this paper, we proposed a new algorithm for identifying dynamical modes of systems with strong nonlinearities (i.e., nonlinearities that do not allow the possibility of linearization) including nonlinearities of the hysteresis type. This algorithm is based on calculating the spectrum of LCEs. For systems with hysteresis blocks that are formalized within phenomenological models (Dahl, Bouc-Wen), the proposed algorithm is based on a modification of the Benettin algorithm. In contrast to the classical algorithm focused on calculating the LCEs of systems with smooth right-hand sides, in the modified algorithm, the Jacobi matrix is calculated on the "smooth sections" of the right-hand sides. Herewith, the calculation of the accumulated sums stops at points corresponding to singularities of the right-hand side; at the next cycle, the Jacobi matrix is calculated on the next section of the smoothness of the right-hand side. In other words, the modified algorithm makes it possible to "bypass" singularities (discontinuity surface, as well as the nondifferentiability surface) of the righthand side while maintaining the accuracy of the LCEs estimation. The paper also proposes a new algorithm for identifying LCEs in systems with the design model of the hysteresis nonlinearity, such as the Preisach operator. This algorithm is also based on a modification of the Benettin algorithm and use the monotonicity property of the Preisach operator, as well as the possibility of smooth approximation of its output at monotonic input.