A POSSIBILITY OF ROBUST CHAOS EMERGENCE IN LORENZ-LIKE NON-AUTONOMOUS SYSTEM

Robust chaos is determined by the absence of periodic windows in bifurcation diagrams and coexisting attractors with parameter values taken from some regions of the parameter space of a dynamical system. Reliable chaos is an important characteristic of a dynamic system when it comes to its practical application. This property ensures that the chaotic behavior of the system will not deteriorate or be adversely affected by various factors. There are many methods for creating chaotic systems that are generated by adjusting the corresponding system parameters. However, most of the proposed systems are functions of well-known discrete mappings. In view of this, in this paper we consider a continuous system that illustrates some robust chaos properties.


Introduction
Chaos theory is a branch of mathematics that deals with complex systems whose behavior is highly sensitive to the initial conditions. Even small alterations can give rise to strikingly great consequences. Tools developed as part of this theory are successfully applied in various branches of science.
Chaos theory has a wide variety of applications. These include but are not limited to: modeling of population dynamics [1], modeling of social processes [2], optimization methods and control theory [3,4], cryptography [5], auxiliary methods in medicine [6,7], predicting stock market dynamics [8].
The general rule of thumb here is that if a given process is known to be deterministic and is exhibiting dynamical behavior, the chaos theory tools can be used to build a model for it. Since most processes in the known Universe possess such attributes, we can safely say that chaos theory is a universal tool that can be applied and yield satisfactory results in almost any problem.
One of the most ubiquitous tools used for modeling are systems of ordinary differential equations (ODE). Such systems allow for flexible parameterization and on-demand scaling. This means that many processes (often different in nature) can be modeled using a single ODE by tweaking its parameters.
Such tweaking of parameters, however, can lead the system into the region of its parameter space where it only exhibits chaotic behavior.
The parameter space is the space of possible parameter values that define the system's behavior. The parameters of ODE form its domain. The ranges of values of the parameters may form the axes of a plot, and particular outcomes of the model may be plotted against these axes to illustrate how different regions of the parameter space produce different types of behavior in the model. These changes of behavioral types are often accompanied by the change in the topology of the objects produced by the system.
It is natural to assume that if the system behaves differently in different points of the parameter space, there exist such points (or surfaces in the case of multidimensional parameter space) that separate one behavioral type from another. These points (surfaces) are called critical.
In 1993, Majumdar and Mitra [11] first coined the phrase robust chaos in dynamic optimization models represented by a quadratic map family. Later in 1996, a search for robust chaos in a discrete-time neural network was conducted by R. Dogaru et al. [12,13] to discover a compact set of parameters, included in a weight space, that could sustain chaotic behaviors but remain unchanged.
Such a definition implies that any changes or variations in system parameters would not result in the fragility of chaos. A practical example of robust chaos in a 2-D piecewise smooth system was also demonstrated through a current-mode controlled boost converter. A search for robust chaos generation approaches has been of considerable interest due to the suitability of robust chaos in practical applications in science and engineering, such as cryptography and secure communications [14] - [19]. Andrecut and Ali [20,21] reconstructed 2-D smooth unimodal maps via non-integer powers for robust chaos by means of mapping a critical point into an unstable fixed point that was not in the basin of attraction of a periodic attractor where, consequently, no periodic attractors occurred. G. Perez [22] has further analyzed the linear interpolation between fully chaotic logistic and quartic maps suggested by S. Thomae, and the results reveal a bifurcation diagram without any periodic windows.This idea was developed in article [23], where superpositions of logistic, sinusoidal and exponential mappings were used as discrete mappings generating robust chaos.
Let us now say that if there exists a subset ω in the parameter space Ω in which: • there are no critical points, • system exhibits chaotic behavior, then the system is said to be robustly chaotic in this subset.
The goal of the chaos robustification is to make such sets as large as possible to ensure that the chaotic behavior will not degrade due to the small changes.
Chaotically robust systems are highly resistant to degradation of chaos since there are no regions in ω where its chaotic behavior vanes. These systems are in high demand in fields where a reliable chaotic signal is required, e.g. in cryptography and secure communication technologies.
There were proposed methods for creating dynamical systems whose chaotic behavior is robust [23]. These systems, however, are discrete in nature and are obtained from the superposition of already existing maps. This can be quite limiting especially when there is a need for creating such systems in bulk. To address this issue we are going to present and explore a system of nonlinear differential equations, i.e. a continuous flow that can be used for generating an unlimited amount of discrete maps.

Exploration methodology
To perform the analysis of the non-autonomous system and demonstrate its chaotic behavior, we need to look at the instruments that are going to be used to explore the properties of the system. Let us begin by introducing some definitions.
Definition 2.1 (Criteria of chaos). A given system is said to be chaotic if: [9] 1. It is sensitive to initial conditions.
3. It has dense periodic orbits.
The definition (2.1) implies that trajectories in the chaotic system that start in close points will eventually diverge. This divergence is captured by the value of the largest Lyapunov exponent (LLE).
The behavior of a given system is defined by the values of its parameters. A combination of these values comprises a point in the system's parameter space.
Definition 2.2 (Robustness of chaos). Robust chaos is defined by the absence of periodic windows and coexisting attractors in some neighborhoods in the parameter space of a dynamical system [10].
Considering these statements, we can formulate the methods of searching and exploring the neighborhoods in the parameter space of the system where it retains its chaotic properties and is resistant to chaos degradation.
To verify that the system's behavior is chaotic at a given point in the parameter space, we are going to estimate its LLE at that point. If the exponent is greater than zero, then the system is chaotic according to the definition (2.1). We are going to use the Wolf algorithm to measure the rate at which close trajectories diverge (converge). This rate represents the LLE value. To get the overall picture of the system's behavior in a given range of parameter values, we are going to plot the estimated values of the largest Lyapunov exponent against the changing parameter.
Since we are not interested in the in-depth exploration of the behavior at each point of the given subspace, the largest Lyapunov exponent will be sufficient to draw necessary conclusions. The use of the Wolf algorithm is therefore justified.
To show that the system is robustly chaotic in the given subspace according to definition (2.2), we are going to use bifurcation diagrams. These diagrams can be used to verify the absence of the windows of periodicity and further show that the behavior of the system is chaotic throughout the given parameter subspace.

Non-autonomous system
This section presents a non-autonomous system of ordinary differential equations and explores its parameter space to find subspaces where the robustness of chaos is observed.

Derivation of the system
This subsection is based on the results presented in the article [24]. The system in question was derived from the equations that describe a convection process in a porous medium subjected to gravitational modulation. It as assumed that the fluid was incompressible and was heated from bellow. The graphical representation of the problem's domain can be seen on the Fig.3.1.
Here T denotes the temperature, v = (u, v) is the velocity, p is the pressure, κ is the coefficient of thermal diffusivity, µ is the kinematic viscosity, g is the gravity, γ is a unit vector in the downward vertical direction (in the sense of gravity), T C is the cold wall temperature, and K is the permeability. The following boundary conditions were considered Here T H is the hot wall temperature, × l and l are the horizontal and the vertical dimensions of the porous medium. The boundary conditions for the temperature correspond to the adiabatic (no-flux) condition at the lateral boundaries and fixed temperature at the lower and upper boundaries. For the velocity, its normal component at the boundary is zero. This means that the fluid does not intersect the boundary. In what follows, the characteristic temperature difference will be denoted by: To obtain the dimensionless model, new spatial variables were introduced x = x/l, y = y/l, time t = (κ/l) 2 t, velocity v(l/κ), pressure (K/µκ)p, and the frequencies σ i = (l 2 /κ)ν i . Denoting θ = (T −T c )/DT and keeping for convenience the same notation for the other variables, the following was obtained 3) here R a = gβDT C Kl/κµ is the Rayleigh number and (u, v) is the velocity vector. The parameter χ = 1/P r D stands for the inverse of Darcy-Prandtl number with P r D = P r /D a , P r = µ/κ is the Prandtl number and D a = K/l 2 is the Darcy number.
The system (3.3)-(3.6) of was supplemented by the free surface boundary conditions for the velocity, zero boundary condition for the temperature at the upper and lower boundaries of the rectangular domain, and adiabatic boundary condition at the side boundaries In order to perform a mathematical analysis of the problem, the stream function ψ was introduced The rotational operator was applied to the system of (3.4)-(3.5) in order to eliminate the pressure To apply the spectral method, a solution of the problem was found in the form of the basic solution added to a variation one. More precisely, the stream function and temperature were separated into a basic conduction part and variation convection one in the form These two representations are equivalent to a Galerkin expansion of the solution in x and y directions, truncated when i + j = 2, where i and j are the Galerkin summation indices in x-direction and in y-direction, respectively. By substituting the form of solutions (3.11) in (3.10), the following was obtained − χ ∂ ∂t + 1 π 2 1 2 + 1 A 11 sin πx sin (πy) −R a (1 + λ 1 sin (σ 1 t) + λ 2 sin(σ 2 t)) π B 11 sin πx sin (πy) = 0. the equation (3.13) can be transformed into π 2 χ dA 11 dt + π 2 A 11 + πR a (1 + λ 1 sin (σ 1 t) + λ 2 sin(σ 2 t)) B 11 × sin πx sin (πy) = 0.   Then dB 11 dt + π 2 B 11 + π A 11 + π 2 A 11 B 02 cos πx sin (πy) (3.18) The following equation was derived by multiplying (3.15) by sin (πx/ ) sin (πy) and integrating over the space domain (0, 1) × (0, ) 4 π 2 χ dA 11 dt + π 2 A 11 + πR a (1 + λ 1 sin (σ 1 t) + λ 2 sin(σ 2 t)) B 11 = 0. (3.19) Again, multiplying (3.18) by cos (πx/ ) sin (πy) and integrating over the space domain Also, multiplying (3.18) by sin (2πy) and integrating over the space domain (3.22) The first equilibrium point of this system is which is expected to be stable when R a /π 2 ς 2 < 1. However, for R a /π 2 ς 2 > 1, he following second equilibrium point is expected to be stable

The influence of modulation on the system's behavior
The parameter space of the system has seven dimensions: α, R, , σ 1,2 , λ 1,2 . We will assume that the first three parameters in this list have the following values then we are left to examine only the parameters that are associated with the influence of the gravitational modulation. The values of the R parameter is going to be used to plot results of the robustification of the system (3.33).
Let us now encapsulate the modulation function in the following way φ(t, σ 1 , σ 2 , λ 1 , λ 2 ) = λ 1 sin (σ 1 t) + λ 2 sin (σ 2 t) , (3.32) then the system (3.30) can be rewritten as followṡ Naturally, the system's equilibrium points will shift as time progresses unless the modulation function is a constant value (which can only be 0 since φ is symmetrical with respect to the t axis). Now let's explore how the stationary points positions and types depend on the values of the modulation function. The parameters of the φ are omitted for more clarity Here, we can divide the first equation by α and the third equation by 4 with an added conditions that = 0 and α = 0 (3.35) From the first and the third equations X and Z can be expressed and from this we have: By substituting X and Z in the second equation of the (3.34) with the values from (3.37) we will get the following equation (3.38) The first root of this equation in respect to Y is the obvious and trivial solution Y 0 = 0, which implies that X 0 = Z 0 = 0. Let us find the other two roots by solving for Y the expression in the parentheses (3.39) Using the discriminant method we get By substituting the (3.41) in (3.37) we get values of X and Z .
Since the modulation function is included in the denominator we have to guarantee that this denominator never equals to zero φ = −1. (3.43) From the equations (3.41-3.42) it is evident that when modulation is not introduced, the equilibrium points are always the same The modulation will cause the points to shift with a certain period and possibly change their types. It implies that the modulation can cause these points to become unstable and repel the system's trajectories. We are going to use this fact to try to robustify the system (3.33).

Robustification
In this section, we are going to present bifurcation diagrams and correspondent Lyapunov exponent estimations for the system (3.33) with various configurations of the modulation parameters.
The behavior of the system is similar to the one of the Lorenz system when modulation is not present. The bifurcation diagram and the LLE estimates for the case when φ = 0 are shown on the Fig.4.2.

Single frequency modulation
Before introducing multifrequency modulation, we are going to explore the properties of the system when single frequency modulation is introduced. The bifurcation diagram and largest Lyapunov exponent estimates can be found on the image below.
As evident from the Fig.4.3, the transition to chaos in the system occurs earlier when the system is subjected to the modulation with the specified parameters. But some occasional periodicity windows still have to be suppressed to further robustify the system.

Multifrequency modulation
The Fig.4.4 shows that the periodicity windows that are present on Fig.4.2-4.3 can be suppressed by introducing the second modulation.
The Fig.4.4 demonstrates that by tweaking the parameters of the modulation one can suppress the periodicity windows and ensure that the behavior of the system remains chaotic for a wider range of the R parameter. The same can be done for all of the other parameters. The process of robustification may continue until the desired subset of the parameter space is found. The final subset of parameters can be used to create maps via discretization of the flow generated by the system. Similar plots can be created for one of the modulation parameters. The Fig.4.5 demonstrates how one of the frequencies of the modulation affects the system's behavior. Such plots can be used to find modulation parameters that turn the system chaotic for given parameters.

Recurrence plots analysis
In the two graphs below, using recurrence analysis (see [25]- [27]), an attempt was made to find out the difference between classical chaos, which the Lorenz system demonstrates, and chaos, which system (3.33) demonstrates. The graphs are represented on For the Lorenz system (see Fig.5.6), adiabatic drift can be noted. Such behavior is typical for systems that do not have uniformity, but contain slowly varying parameters. In Fig. 5.7, this phenomenon is noted by a change in image brightness when removed in a direction perpendicular to the main diagonal. Note that for system (3.33), the same structure of the recurrence diagram (despite noticeable local differences) is preserved. for λ1 = 10, λ2 = 1 As can be seen from the comparison of Fig. 5.7 and 5.8 there is a drift in the behavior of the system. The structure of recurrence diagrams with an increase of λ 1 becomes more uniform. This condition is characteristic of non-stationary systems, but containing slowly changing parameters. In recurrence plots, this fact is noted by a characteristic change in the brightness of the image from the lower right to the upper left corner. for λ1 = 50, λ2 = 1