Avoiding Dynamical Degradation in Computer Simulation of Chaotic Systems Using Semi-Explicit Integration: Rössler Oscillator Case

: Dynamical degradation is a known problem in the computer simulation of chaotic systems. Data type limitations, sampling, and rounding errors give rise to the periodic behavior. In applications of chaotic systems in secure communication and cryptography systems, such effects can reduce data storage security and operation. In this study, we considered a possible solution to this problem by using semi-explicit integration. The key idea is to perturb the chaotic trajectory by switching between two integrators, which possess close but still different numerical solutions. Compared with the traditional approach based on the perturbation of the bifurcation parameter, this technique does not signiﬁcantly change the nonlinear properties of the system. We verify the efﬁciency of the proposed perturbation method through several numerical experiments using the well-known Rössler oscillator.


Introduction
In recent decades, chaotic systems have been widely used in many engineering and scientific tasks, including secure communication [1][2][3][4], neural networks [5,6], watermarking [7], cryptography [8][9][10], robotics [11,12], sensors [13], signal processing [14], etc. Usually, such practically oriented studies and applications are performed using computer simulation. Converting mathematical models to executable ones inevitably causes an imprecision. In addition to the identification error of the mathematical model, the discrepancy between the numerical results and the prototype behavior is enhanced by the data type precision, discretization methods error, and round-off results of arithmetic operations [15]. In chaotic systems, the combination of these factors can lead to the degradation of chaotic dynamics [16]. In [17], K.J. Persohn and R.J. Povinelli explicitly demonstrated that for four-bit binary fractions representation, the trajectory of the known logistic map starts repeating after two iterations with a period of three for any initial value.
Several methods have been proposed to avoid cycles in the chaotic trajectories obtained in computer simulations [16,[18][19][20][21][22][23][24][25]. One of the most straightforward and easily implementable approaches is perturbing the trajectory of a discrete chaotic system or the bifurcation parameter [16,[18][19][20]. Usually, the cycle breaking is performed using an auxiliary discrete chaotic map of a small dimension. A significant period prolongation with this technique was demonstrated for delayed chaotic maps [21]. To determine the disturbance moments, a linear feedback shift register (LFSR) or a chaotic analog system can be used [22][23][24]. In [25], it was proposed to introduce perturbation only when the trajectory becomes a loop. This approach requires additional memory.
All of the above approaches effectively cope with the period increasing of the trajectory of a chaotic system in long-term simulation. These techniques are especially useful in cryp-

•
the new technique for extending the period of chaotic sequences is given; • the proposed method does not introduce additional perturbations to the chaotic oscillations in comparison with traditional methods based on switching nonlinearity parameters; • the considered approach can be implemented in embedded systems without using additional hardware resources.
The rest of the paper is organized as follows. Section 2 describes the proposed technique and shows how it can be applied to the sample chaotic system. Section 3 presents the algorithm for detecting cycles in numerical sequences. We experimentally investigated the proposed method in comparison with a common technique based on parameter switching. Finally, some conclusions and discussions are given in Section 4.

Materials and Methods
A common practice for avoiding cycles in the trajectories of a chaotic system in computer simulations is switching the values of the bifurcation parameters. In this approach, two main factors prevent looping. First, the steady-state oscillation mode can be violated, and an undesirable situation can arise when the trajectory leaves the attractor or the chaotic behavior changes to a periodic mode. The second aspect concerns floating-point calculations, which involve normalizing numbers in binary operations [32]. Two bifurcation parameter values can be normalized in different ways, which leads to various errors at each integration step. We employed a similar trick of switching between two semi-explicit finite-difference schemes that simulate the same chaotic system in the proposed approach.
Let us consider the idea of semi-explicit integration first.

Semi-Explicit Integration
Semi-implicit integration implies the reuse of already approximated values of state variables at the (n + 1) step to calculate the rest of the state variables within one integration step [32,33]. Being compared with explicit and implicit numerical methods, semi-explicit methods can include additional terms of the Taylor series. The set and number of such expressions depend on the calculation order of the right-hand sides of the equations.
The simplest semi-explicit integration technique is the first-order Euler-Cromer method, also known as the semi-explicit Euler or symplectic Euler method [33,34]. It can be introduced for the two-dimensional Hamiltonian system: where f and g are given right-hand side functions. The Euler-Cromer method produces an approximate discrete solution for System (1) by iterating: where h is the integration step. The key difference the Euler-Cromer method (2) makes following a conventional explicit Euler method is that the equation for v n+1 uses an already calculated value of x n+1 . Such calculation inevitably leads to the appearance of an additional term in the Taylor expansion of variable v. If we reverse the order of calculations, we get the modification of the Euler-Cromer method of the same order of accuracy: with an extra term in the Taylor series for variable x. Thus, by changing the sequence of the approximation of the state variables, we can vary the number of additional terms of the Taylor expressions when solving the differential equations of arbitrary order. One can notice that for a system consisting of N equations, in the general case, there are N! variants of finite-difference schemes that can be obtained with the Euler-Cromer (3) integration. Let us consider how one can use this feature of semi-explicit methods to avoid chaotic degradation in discrete simulations of continuous chaotic systems.

Perturbation Technique
The proposed technique is based on switching between two finite-difference schemes obtained through the Euler-Cromer integration with different calculations. The Taylor series reconstructed by such methods coincide up to the first term but are different in the remaining terms. Therefore, switching between two various schemes obtained via the Euler-Cromer integration can significantly vary the truncation error accumulated on each integration step. In addition to the methodical error, the different order of calculations also introduces some extra perturbation since floating-point operations are not completely associative [32].
Thus, we propose applying the perturbation scheme shown in Figure 1 to avoid the cycles in the computer simulation of chaotic systems. At each integration step, the switching law generates 0 or 1, which corresponds to the choice of the corresponding Euler-Cromer integration method with a certain order of calculations. For example, one can use the forward and backward approximation, as is shown in Figure 1.
To determine the moment of switching, any of the known approaches can be applied, including LFSR or a pseudo-random bit generator (PRBG) [24]. According to the LFSR rule, each finite-difference scheme is matched to the switching law output. Figure 2 represents one of the possible schemes with LFSR given by a primitive polynomial x 8 + x 4 + x 3 + x + 1 with initial state (0, 0, 0, 0, 0, 0, 0, 1).  It should be noted that the developed perturbation technique can be applied using high-order methods as well, including the semi-implicit composition and extrapolation methods [35,36]. However, in practical applications-e.g., communication systems based on chaotic synchronization phenomenon-it is required to use finite-difference schemes with a minimum number of arithmetic operations [37]. The use of high-order methods in solving such problems leads to a decrease in the data transfer rate. In addition, in embedded systems the dimension of the implemented finite-difference scheme is limited by hardware resources, since multiple chaos generators and receivers are implemented. Therefore, the present study focused on studying methods of a small order of algebraic accuracy that satisfy this criterion.
Let us consider the results of applying the proposed approach to simulate a chaotic system.

Experimental Results
As the sample system, we use the well-studied Rössler oscillator [38] described by the following system of ordinary differential equations (ODE): where a, b, and c are bifurcation parameters.
Let us consider the finite-difference scheme for System (4) obtained by the Euler-Cromer integration with a straightforward order of state variables calculation: To determine the moment of switching, any of the known approaches can be applied, including LFSR or a pseudo-random bit generator (PRBG) [24]. According to the LFSR rule, each finite-difference scheme is matched to the switching law output. Figure 2 represents one of the possible schemes with LFSR given by a primitive polynomial x 8 + x 4 + x 3 + x + 1 with initial state (0, 0, 0, 0, 0, 0, 0, 1). To determine the moment of switching, any of the known approaches can be applied, including LFSR or a pseudo-random bit generator (PRBG) [24]. According to the LFSR rule, each finite-difference scheme is matched to the switching law output. Figure 2 represents one of the possible schemes with LFSR given by a primitive polynomial x 8 + x 4 + x 3 + x + 1 with initial state (0, 0, 0, 0, 0, 0, 0, 1).  It should be noted that the developed perturbation technique can be applied using high-order methods as well, including the semi-implicit composition and extrapolation methods [35,36]. However, in practical applications-e.g., communication systems based on chaotic synchronization phenomenon-it is required to use finite-difference schemes with a minimum number of arithmetic operations [37]. The use of high-order methods in solving such problems leads to a decrease in the data transfer rate. In addition, in embedded systems the dimension of the implemented finite-difference scheme is limited by hardware resources, since multiple chaos generators and receivers are implemented. Therefore, the present study focused on studying methods of a small order of algebraic accuracy that satisfy this criterion.
Let us consider the results of applying the proposed approach to simulate a chaotic system.

Experimental Results
As the sample system, we use the well-studied Rössler oscillator [38] described by the following system of ordinary differential equations (ODE): where a, b, and c are bifurcation parameters.
Let us consider the finite-difference scheme for System (4) obtained by the Euler-Cromer integration with a straightforward order of state variables calculation: It should be noted that the developed perturbation technique can be applied using high-order methods as well, including the semi-implicit composition and extrapolation methods [35,36]. However, in practical applications-e.g., communication systems based on chaotic synchronization phenomenon-it is required to use finite-difference schemes with a minimum number of arithmetic operations [37]. The use of high-order methods in solving such problems leads to a decrease in the data transfer rate. In addition, in embedded systems the dimension of the implemented finite-difference scheme is limited by hardware resources, since multiple chaos generators and receivers are implemented. Therefore, the present study focused on studying methods of a small order of algebraic accuracy that satisfy this criterion.
Let us consider the results of applying the proposed approach to simulate a chaotic system.

Experimental Results
As the sample system, we use the well-studied Rössler oscillator [38] described by the following system of ordinary differential equations (ODE): where a, b, and c are bifurcation parameters. Let us consider the finite-difference scheme for System (4) obtained by the Euler-Cromer integration with a straightforward order of state variables calculation: Substituting x n+1 into the other equations and making a cancellation, one can obtain: It can be seen from Equation (6) that the solution of System (4) is accurately approximated to the first term of the Taylor series, while the resulting expansions partially contain terms from the rest of the series. Applying the Euler-Cromer method with reverse calculation order, we obtain another finite-difference scheme: which can be rewritten as: As one can see from Formulas (6) and (8), the reconstructed series coincide up to the first term of the Taylor series but differ in the remaining terms.

Comparison with Known Perturbation Techniques
In the first experimental part of our study, we compared the proposed approach with known techniques of avoiding cycles in chaotic sequences. First, we considered the widely used method of switching the bifurcation parameter. The second investigated method was changing the order of calculations in the right-hand side function of the differential equations. For example, for System (4), the equation of the z variable can be rewritten as: Nepomuceno et al. found that such a trick in chaotic systems leads to entirely different chaotic trajectories [39]. It is expected that this approach will demonstrate poorer performance than the conventional parameter perturbation and the proposed method. However, this approach can assess the degree of influence of rounding errors of different numbers or the order of operations on the extension of the sequence period.
For the approaches under comparison, the switching moments are determined by the LFSR shown in Figure 2. We estimated the period length of the trajectories of the Rössler oscillator simulated with a 16-bit fixed-point data type with an 8-bit set for the word length. The simulation was performed on time interval T = 300 s and the integration step equaled 0.0117188 s. We considered three sets of initial conditions. Each sample set was generated as follows. We equated one of the initial conditions to 0.96875. The rest of the two initial values were in the range [−0.0315; 0.96875] with an increment of 0.00390625, which coincided with the minimum representable value in the chosen 16-bit fixed-point data type.
In the parameter-perturbation method, switching was performed for parameter b between values 0.1875 and 0.210938, for which chaotic behavior was observed. The approach that is based on a different order of operations uses the original Rössler system (4) and its representation (9). In both cases, simulations were performed by the Euler-Cromer method with the forward order of computation. In the proposed technique, we applied finite-difference schemes (5) and (7), respectively.
To detect cycles' appearance in chaotic sequences, we used the algorithm described in [40]. The experimental results obtained with the set of initial conditions where y 0 is fixed are given in Figure 3. The number of periodic sequences for the sample set with y 0 = 0.96875 for different perturbation techniques is shown in Table 1. The estimates obtained with the other two sets are presented in Appendix A (Figures A1 and A2). The abscissa and ordinate axes in Figures 3, A1 and A2 correspond to the values of the initial conditions from which the simulation began. Pairs of initial values for which equal numbers were found in the generated sequence are marked black. It can be noted that the number of pairs of initial conditions leading to periodic solutions was reduced for all the considered techniques compared to the original model (case (a) in Figures 3, A1 and A2). The efficiency of the proposed approach was close to the traditional method based on switching the bifurcation parameter (cases (c) and (d) in Figures 3, A1 and A2). In addition, it can be seen that by changing the order of calculations only, one cannot obtain the required period extension (case (b) in Figures 3, A1 and A2). It can be seen that none of the considered techniques disturbed the trajectory enough to change this mode to chaotic behavior. that is based on a different order of operations uses the original Rössler system (4) and its representation (9). In both cases, simulations were performed by the Euler-Cromer method with the forward order of computation. In the proposed technique, we applied finite-difference schemes (5) and (7), respectively. To detect cycles' appearance in chaotic sequences, we used the algorithm described in [40]. The experimental results obtained with the set of initial conditions where y0 is fixed are given in Figure 3. The number of periodic sequences for the sample set with y0 = 0.96875 for different perturbation techniques is shown in Table 1. The estimates obtained with the other two sets are presented in Appendix A (Figures A1 and A2). The abscissa and ordinate axes in Figures 3, A1 and A2 correspond to the values of the initial conditions from which the simulation began. Pairs of initial values for which equal numbers were found in the generated sequence are marked black. It can be noted that the number of pairs of initial conditions leading to periodic solutions was reduced for all the considered techniques compared to the original model (case (a) in Figures 3, A1 and A2). The efficiency of the proposed approach was close to the traditional method based on switching the bifurcation parameter (cases (c) and (d) in Figures 3, A1 and A2). In addition, it can be seen that by changing the order of calculations only, one cannot obtain the required period extension (case (b) in Figures 3, A1 and A2). It can be seen that none of the considered techniques disturbed the trajectory enough to change this mode to chaotic behavior.

Perturbation Technique Number of Periodic Sequences
Original 16-bit model 10,038 Switching between two forms of the right-hand side functions 428 Perturbation of the bifurcation parameter 13 Proposed technique 21

Comparison with Other Finite-Difference Schemes
In the second part of the experimental study, we checked whether the perturbation was caused by different additional terms of the Taylor series in switchable semi-explicit schemes. We used the same flowchart shown in Figure 1 but applied other known integration methods, namely the traditional Euler integration, the explicit midpoint method, and the Heun (Runge-Kutta 2) algorithm.
For the sample System (4), the conventional Euler integration gives the following finite-deference scheme: Formula (10) approximates the solution of the Rössler system up to the first term of the Taylor series. In the perturbation technique, the Euler method can be switched with its semi-explicit counterpart.
It can be assumed that the perturbation can also result from switching methods of different orders. Therefore, as the third case in our comparison, we include Euler and Heun methods. We repeated the estimate of the period length for the sample set of initial conditions with y 0 = 0.96875 and x 0 , z 0 ∈ [−0.0315; 0.96875]. The obtained results are shown in Figure 4. The numbers of sequences with cycles are listed in Table 2.
Both methods (11) and (12) are explicit and possess the algebraic accuracy of second order. The finite-difference schemes obtained via these techniques restore the Taylor series precisely up to the second term. Therefore, this pair of methods can potentially be used together to simulate a chaotic system.
It can be assumed that the perturbation can also result from switching methods of different orders. Therefore, as the third case in our comparison, we include Euler and Heun methods.
We repeated the estimate of the period length for the sample set of initial conditions with y0 = 0.96875 and x0, z0∈ [−0.0315; 0.96875]. The obtained results are shown in Figure  4. The numbers of sequences with cycles are listed in Table 2.

Switchable Methods Number of Periodic Sequences
Two Euler-Cromer methods 21

Euler-Cromer and Euler methods 1183
Explicit midpoint and Heun methods 2608 Heun and Euler methods 599 As one can see, the efficiency of the approach based on methods of different orders was close to the technique based on switching various forms of the right-hand sides of the differential equation. Additional terms of the Taylor series, reconstructed by the Euler-Cromer method, did not provide the same perturbation when paired with the Euler method than with the semi-explicit copy with a different order of calculations. The worst performance was shown by switching between the Heun and the midpoint methods. This can be explained by the relatively low difference between the solutions produced by these Thus, we can conclude that the best perturbation is provided by switching semiexplicit methods among the approaches considered.

Conclusions and Discussion
In this study, a new approach to avoid dynamical degradation in chaotic systems was proposed. It was theoretically predicted that applying a pair of finite-difference schemes obtained by the semi-explicit integration method can help prevent chaotic degradation in the generated sequences. This phenomenon can be explained as follows: slight changes in the truncation errors and round-off errors of arithmetic operations result in simultaneous changes in the discrete system's trajectory. Experimental results confirmed the theoretical suggestions and showed that it is possible to extend the period of a chaotic sequence by switching the integration order in computations performed with short data types. The efficiency of the proposed approach was comparable to the traditional technique based on the perturbation of the bifurcation parameter. However, an indisputable advantage is the invariability of the bifurcation properties of the model during the perturbation. It can be helpful in various types of analysis that require long-term simulation and secure communication systems, where the concealment of transmission is preferred. In addition, semi-explicit integrators can be compactly implemented in embedded communication systems based on chaotic synchronization. In such systems, to avoid dynamic degradation, in the transmitter, one can switch the input values of the Euler-Cromer integrator as shown in Figure 5, obtaining different orders of calculations. In comparison, for example, from the switching between the explicit midpoint method (11) and the Heun method (12), to implement a pair of semi-explicit methods, no additional hardware resources are required since the functions of the right-hand side are calculated almost identically.
As one can see, the efficiency of the approach based on methods of different orders was close to the technique based on switching various forms of the right-hand sides of the differential equation. Additional terms of the Taylor series, reconstructed by the Euler-Cromer method, did not provide the same perturbation when paired with the Euler method than with the semi-explicit copy with a different order of calculations. The worst performance was shown by switching between the Heun and the midpoint methods. This can be explained by the relatively low difference between the solutions produced by these second-order integrators. Only the different order of algebraic operations caused the perturbation in this case.
Thus, we can conclude that the best perturbation is provided by switching semi-explicit methods among the approaches considered.

Conclusions and Discussion
In this study, a new approach to avoid dynamical degradation in chaotic systems was proposed. It was theoretically predicted that applying a pair of finite-difference schemes obtained by the semi-explicit integration method can help prevent chaotic degradation in the generated sequences. This phenomenon can be explained as follows: slight changes in the truncation errors and round-off errors of arithmetic operations result in simultaneous changes in the discrete system's trajectory. Experimental results confirmed the theoretical suggestions and showed that it is possible to extend the period of a chaotic sequence by switching the integration order in computations performed with short data types. The efficiency of the proposed approach was comparable to the traditional technique based on the perturbation of the bifurcation parameter. However, an indisputable advantage is the invariability of the bifurcation properties of the model during the perturbation. It can be helpful in various types of analysis that require long-term simulation and secure communication systems, where the concealment of transmission is preferred. In addition, semi-explicit integrators can be compactly implemented in embedded communication systems based on chaotic synchronization. In such systems, to avoid dynamic degradation, in the transmitter, one can switch the input values of the Euler-Cromer integrator as shown in Figure 5, obtaining different orders of calculations. In comparison, for example, from the switching between the explicit midpoint method (11) and the Heun method (12), to implement a pair of semi-explicit methods, no additional hardware resources are required since the functions of the right-hand side are calculated almost identically. It is of interest to extend the proposed perturbation technique using methods with higher orders of algebraic accuracy. By increasing the order of the integrator, on the one It is of interest to extend the proposed perturbation technique using methods with higher orders of algebraic accuracy. By increasing the order of the integrator, on the one hand, one can reduce the truncation error due to the reconstruction of a larger number of terms of the Taylor series. However, this process is associated with an increase in calls to the right-hand side function, which affects the number of arithmetic operations and the round-off error in one-step integration. As shown by [42], low-order methods can yield a large error while simulating with a floating-point data type. In such cases, more than half of the mantissa bits have pseudo-random properties even in solving linear ordinary differential equations. Therefore, reversing such numerical effects against cycles in chaotic sequences can be helpful in various applications. Moreover, it is interesting to test the proposed technique for simulation of other chaotic systems, which will also be the subject of our further research.
hand, one can reduce the truncation error due to the reconstruction of a larger number of terms of the Taylor series. However, this process is associated with an increase in calls to the right-hand side function, which affects the number of arithmetic operations and the round-off error in one-step integration. As shown by [42], low-order methods can yield a large error while simulating with a floating-point data type. In such cases, more than half of the mantissa bits have pseudo-random properties even in solving linear ordinary differential equations. Therefore, reversing such numerical effects against cycles in chaotic sequences can be helpful in various applications. Moreover, it is interesting to test the proposed technique for simulation of other chaotic systems, which will also be the subject of our further research. Acknowledgments: The authors are thankful to anonymous reviewers for their valuable suggestions and fruitful discussion during the review process.

Conflicts of Interest:
The authors declare no conflict of interest.