Dynamics, Periodic Orbit Analysis, and Circuit Implementation of a New Chaotic System with Hidden Attractor

Abstract

Hidden attractors are associated with multistability phenomena, which have considerable application prospects in engineering. By modifying a simple three-dimensional continuous quadratic dynamical system, this paper reports a new autonomous chaotic system with two stable node-foci that can generate double-wing hidden chaotic attractors. We discuss the rich dynamics of the proposed system, which have some interesting characteristics for different parameters and initial conditions, through the use of dynamic analysis tools such as the phase portrait, Lyapunov exponent spectrum, and bifurcation diagram. The topological classification of the periodic orbits of the system is investigated by a recently devised variational method. Symbolic dynamics of four and six letters are successfully established under two sets of system parameters, including hidden and self-excited chaotic attractors. The system is implemented by a corresponding analog electronic circuit to verify its realizability.

1. Introduction

Chaos theory, which is regarded as the third scientific theory revolution in the 20th century, has been extensively and intensively studied since the meteorologist Lorenz discovered chaotic phenomena for three-dimensional (3D) autonomous quadratic systems in 1963 [1]. As chaotic states in nonlinear dynamic systems are extremely sensitive to initial values, a large amount of research work shows that chaos is closely related to engineering technology, with wide application in fields such as circuit control [2], image encryption [3], secure communications [4], and neural networks [5].

Many chaotic systems have been constructed [6,7,8,9] that include both self-excited and hidden attractors [10]. Self-excited attractors have a basin of attraction related to the unstable equilibrium, whereas those of hidden attractors do not intersect with small neighborhoods of any equilibria [11,12]. Most well-known dynamical systems have self-excited chaotic attractors [13,14,15]. As hidden attractors cannot be calculated from the initial conditions in the neighborhood of the equilibrium point, they were not introduced until recently, and there are some new studies on how to locate them [16,17]. Hidden attractors have attracted great interest in recent years due to their considerable importance in both theory and engineering, because they allow unexpected and potentially catastrophic responses to structural disturbances such as to bridges or aircraft wings. It has been shown that attractors in a dynamical system with stable equilibria [18,19,20], an infinite number of equilibria [21,22,23], or no equilibrium points [24,25,26,27,28,29] are hidden attractors. These are also represented in a 3D continuous dynamical system with only one unstable node as the equilibrium point [30].

Hidden attractors have been broadly investigated in the literature. Wang and Chen constructed a chaotic system with only one stable equilibrium via a constant control parameter added to the Sprott E system [31]. Wei found a new chaotic system with no equilibrium by adding a simple constant to the Sprott D system [32]. Two modified Sprott systems that have only stable node-focus points with hidden chaotic attractors were analyzed [33]. The generalized Sprott C system with two stable equilibrium points was proposed [34], and its chaotic and complex dynamic behaviors in the parametric space were investigated. Hidden chaos and hyperchaos have been found in the jerk system [35,36,37], meminductor-based chaotic system [38], extended Rikitake system [39], 4D Rabinovich system [40], 4D modified Lorenz-Stenflo system [41], and 5D homopolar disc dynamo [42]. Self-excited and hidden attractors can also be generated in a modified Chua’s circuit and in a 3D memristive Hindmarsh–Rose neuron model [43,44]. Moreover, some 3D dynamical systems with three different families of hidden attractors have been discovered [45,46,47]. A 4D autonomous chaotic system that has two types of hidden attractors with a line of equilibria or no equilibria was derived [48]. A 5D chaotic system with both hidden attractors and extreme multistability was introduced [49], and coexisting self-excited and hidden attractors in a Lorenz-like system with two equilibria were found [50].

This paper proposes a hidden chaotic attractor system with two stable fixed points. With the change of parameters, its complex dynamical behaviors are analyzed using multiple dynamical tools, such as phase portraits, time sequence, power spectrum, and Lyapunov exponents. We establish two symbolic dynamics in the system, and classify the unstable periodic orbits embedded in hidden and self-excited chaotic attractors topologically for two sets of parameters. The electronic circuit of the system is designed and simulated by Multisim software, which proves the existence of chaos. Compared to the above contributions in the literature, the novelty of the work lies in the unstable periodic orbits of the new system showing a complexity of significant differences with different parameter values. For the system with hidden chaotic attractors, which is determined by four parameters, the complexity is relatively simple; however, the system with self-excited chaotic attractors, which contains only two parameters, unexpectedly has more complex dynamics.

The rest of this paper is arranged as follows. Section 2 introduces the mathematical model of the system, and its nonlinear dynamic characteristics are investigated. In Section 3, observation of chaotic and complex dynamics in the system is implemented by varying different parameters. To locate the unstable cycles in the system, we review the variational method in Section 4, which can be effectively utilized in calculations. We systematically calculate all short unstable cycles of the new system under two parameters. To establish appropriate symbolic dynamics, one needs four letters, and the other needs six. Section 5 presents a circuit implementation of the system to validate its feasibility. Section 6 discusses our conclusions.

2. The New System and Its Dynamic Characteristics

Inspired by the chaotic system proposed by [51], a new system can be easily constructed by adding a nonlinear term of cross-product $kxz$ to the first equation,

$$\begin{array}{ccc}\hfill \frac{dx}{dt}& =& a(y-x)+kxz\hfill \\ \hfill \frac{dy}{dt}& =& -cy-xz\hfill \\ \hfill \frac{dz}{dt}& =& -b+xy,\hfill \end{array}$$

(1)

where $x,y,$ and z are state variables, and a, b, c, and k are the control parameters. Note that adding a cross-product nonlinear term is not a general method to realize chaos in a 3D quadratic system. In addition, the proposed system (1) has three nonlinear terms, where notably each equation has one single cross-product term, so it certainly does not belong to algebraic simple chaotic flows, but is suitable for practical implementation as an electronic circuit. When the parameters of system (1) are assigned as $(a,b,c,k)=(10,100,11.2,-0.2),$ and the initial values $(x_0,y_0,z_0)$ are set as $(1,1,1)$, a fourth-order Runge–Kutta method is adopted in the numerical integration, which reveals the chaotic behaviors characterized by strange attractors, as shown in Figure 1a–c, and the power spectrum with continuous broadband characteristics (Figure 1d) verifies the emergence of chaos. Correspondingly, the Lyapunov exponents are calculated based on the Wolf algorithm [52], which gives $L{E}_1=0.7457,L{E}_2=-0.0057,L{E}_3=-26.8144$ (see Figure 2). The positive Lyapunov exponent indicates that the phase volume of the system is expanding and folding in a certain direction, which means that the system is in a chaotic state. The Kaplan–Yorke dimension is $D_{KY}=2+(L{E}_1+L{E}_2)/\left|L{E}_3\right|=2.0276$, which also verifies the chaoticity of system (1).

The new system has the following fundamental dynamic properties:

(1) System (1) is rotationally symmetric versus the z-axis, which is invariant under the coordinate transformation $(x,y,z)\to (-x,-y,z)$. Any attractors are either a symmetric pair or symmetric under a 180-degree rotation around the z-axis;

(2) Since the divergence of system (1) is

$$\nabla ·V=\frac{\partial \stackrel{·}{x}}{\partial x}+\frac{\partial \stackrel{·}{y}}{\partial y}+\frac{\partial \stackrel{·}{z}}{\partial z}=-a+kz-c,$$

(2)

under the condition $-a+kz-c<0,$ it is dissipative and can converge to a set of zero measure in exponential form,

$$\frac{dV}{dt}=e^{-a+kz-c};$$

(3)

(3) System (1) possesses two equilibrium points:

$$\begin{array}{ccc}& & E_1:(-\sqrt{\frac{ab-bck}{a}},-\frac{b}{\sqrt{\frac{ab-bck}{a}}},\frac{ac}{ck-a}),\hfill \\ & & E_2:(\sqrt{\frac{ab-bck}{a}},\frac{b}{\sqrt{\frac{ab-bck}{a}}},\frac{ac}{ck-a}).\hfill \end{array}$$

(4)

Linearizing the system gives the Jacobian matrix

$$A=\left(\begin{array}{ccc}-a+kz& a& kx\\ -z& -c& -x\\ y& x& 0\end{array}\right).$$

(5)

For the parameters $(a,b,c,k)=(10,100,11.2,-0.2),$ the Jacobian eigenvalues at $E_1$ and $E_2$ can be calculated by solving the corresponding characteristic equations, and they have the same values: ${\lambda }_1=-18.7413,$${\lambda }_{2,3}=-0.314\pm 11.424i.$ From the eigenvalues, we can see that $E_1$ and $E_2$ are both stable node-focus points.

Since the new system can generate strange attractors, it is implied that system (1) under current parameters has hidden chaotic attractors. Figure 3a displays the coexistence of chaotic motion and stable node-foci in 3D phase space. We can see clearly that the trajectory starting from initial conditions $I_1=(1,1,1)$ becomes a disordered state; however, orbits starting from initial conditions $I_2=(-20,-10,-10)$ and $I_3=$$(20,10,-10)$ spirally converge to $E_1$ and $E_2$, respectively. Figure 3b displays the time-domain waveform diagram for initial conditions $I_1,I_2$, and $I_3$; an apparently chaotic waveform of $x\left(t\right)$ illustrates that hidden chaotic attractors exist in system (1). To avoid transient chaos, we also confirmed the existence of a hidden chaotic attractor, since the orbit remains on it for $t={10}^6$.

The graphics of the basin of attraction, which is defined as the initial condition set that the orbits converge to a given attractor, can clearly exhibit the initial point distributions of different attractors. To further check whether the chaotic attractor in Figure 3 is hidden, a section $z=-9.1503$, including $E_1$ and $E_2$, is selected, and the initial condition regions of coexisting attractors are explored, as shown in Figure 4. Three types of basins of attraction on the cross section are colored yellow, blue, and red. Yellow areas with black stripes represent the basin of attraction of a chaotic attractor and the Poincaré section concerning the chaotic attractor with two wings, while blue and red areas denote that the movement from these initial conditions will converge to equilibria $E_1$ and $E_2$, respectively. From Figure 4, we find that the basin of attraction has the expected symmetric similarity and a smooth boundary. Moreover, in view of the topology of basins in Figure 4, the basins of attraction of chaotic attractors do not intersect with small neighborhoods of stable equilibrium points $E_1$ and $E_2$, which also illustrates that there is a hidden chaotic attractor in system (1).

3. Chaotic and Complex Dynamics in New System

The system parameters can significantly influence the dynamics, and the qualitative or topological variety in the behavior of dynamic systems means that a bifurcation occurs [53]. We discuss the chaotic and complex dynamics of the proposed system (1) through varying the parameters, taking the initial values as $(x_0,y_0,z_0)=(1,1,1)$. The bifurcation diagram, largest Lyapunov exponent, and division diagram are adopted as tools to observe the impacts of parameters.

3.1. Fix $a=10,c=11.2,k=-0.2$ and Vary b

To investigate the effect of the parameters on the dynamics of system (1), we first take parameters $(a,c,k)=(10,11.2,-0.2)$ and vary $b\in [10,140]$. When altering b, the system shows many complex dynamic behaviors that can be explored in the parameter space. The typical Benettin method is employed to calculate the maximum Lyapunov exponent spectrum, and the corresponding bifurcation diagram versus parameter b is displayed in Figure 5. Obviously, a positive maximum Lyapunov exponent implies that the system is chaotic over a wide range of parameters. It is clear that the bifurcation diagram with the variation of parameter b matches well with the largest Lyapunov exponent spectrum. We can see that the state of system (1) becomes chaotic through pitchfork and period-doubling bifurcations, then periodic, and chaotic again. The bifurcation diagram in Figure 5b demonstrates that the system evolves smoothly from a periodic solution to a chaotic region through a typical period-doubling route; hence, no clear boundary exists between a periodic phase portrait and chaos. When b is in the interval $[125,140]$, the largest Lyapunov exponent quickly becomes negative, and the corresponding bifurcation diagram suddenly changes to no cutoff points, which both mean that the trajectory of system (1) finally converges to a fixed point. More details are presented through the 3D projections of phase portraits of system (1) at different b values, as shown in Figure 6.

3.2. Fix $b=100,c=11.2,k=-0.2$ and Vary a

Now, we fix $(b,c,k)=(100,11.2,-0.2)$ and vary $a\in [6,20]$ and explore the dynamical evolution of system (1). Figure 7a shows the largest Lyapunov exponent spectrum with respect to a, and Figure 7b displays the bifurcation diagram of the whole evolution process. As can be seen from Figure 7, these results are consistent with each other and demonstrate that the dynamical behaviors vary when a undergoes change. Obviously, when a is in the interval $[6,9.7]$, the system converges to one stable equilibrium, as shown in Figure 8a, where $a=7$. When a$\in (9.7,12.1]$, the system has a chaotic status. Near $a=12.2$, the largest Lyapunov exponent is about zero, which implies that system (1) is periodic in a small parameter range, as demonstrated in Figure 8b, where $a=12.2.$ However, the system becomes chaotic again when a is in the interval $[12.6,13.7]$ (see Figure 8c). Then, when $a>13.7$, the system experiences an inverse period-doubling bifurcation process with the increase of a, and eventually becomes periodic again. Figure 8d displays the 3D phase diagram for $a=20$.

3.3. Fix $a=10,b=100,k=-0.2$ and Vary c

Here, we fix the parameters $a=10,b=100,k=-0.2$, and vary c. The Lyapunov exponent spectrum and bifurcation diagram presented in Figure 9 reveal that the limit cycle, chaos, and equilibrium point appear alternately as c increases from −30 to 20. We can see that the system shows complex dynamic behavior in this region, and generates chaos via period-doubling bifurcation. Moreover, periodic windows exist in such a parameter region. An attractor of system (1) becomes a limit cycle from chaos through a process of reverse period-doubling bifurcations, then becomes chaotic again through period-doubling bifurcations, and finally converges to a stable equilibrium point. With initial values $(x_0,y_0,z_0)=(-1,-1,1)$, the system will undergo the same bifurcation process, and the orbits in the phase space will eventually converge toward another stable fixed point.

When the parameters $a=10,b=100,k=-0.2$, the emergence of a hidden chaotic attractor is dependent on the value of $c>0$. When we take $c\in [-10.7,0]$, the system generates a self-excited chaotic attractor. It is worth noting that the periodic regions do not share the same dynamical characteristics; diverse limit cycles appear in four periodic regions, as shown in Figure 10.

3.4. Fix $a=10,b=100,c=11.2$ and Vary k

The Lyapunov exponent spectrum and bifurcation diagram shown in Figure 11 reveal that equilibrium point and chaotic orbit appear alternately with k increasing gradually from $-2$ to $0.9$. When we take parameters $(a,b,c,k)=(10,100,11.2,0)$, system (1) becomes the system with a hidden chaotic attractor in Ref. [51]. As shown in Figure 11, when k becomes slightly positive or negative, chaotic attractors can still be generated. However, the type of chaotic attractors will depend on the positive or negative value of k; these are the new structures that have emerged. As required by the Routh–Hurwitz stability criterion, when we take parameters $(a,b,c)=(10,100,11.2)$, hidden chaotic attractors can exist if the following inequalities are satisfied: $a+c-\frac{ack}{ck-a}>0$, which requires $k<0$ or $k>1.69$. Hence, when k is slightly positive, the chaotic attractor is self-excited, while when k is slightly negative, it generates a hidden chaotic attractor with two stable node-foci. It is noteworthy that, with k increasing in the range $k>0.9$, the orbit finally leads to infinity.

3.5. Fix $a=10,b=100$ and Vary k and c

We draw a division diagram to capture different kinds of dynamical modes of system (1) with respect to parameters k–c. Varying k and c within the region of $k\in [-4,0]$, $c\in [0,25]$, by calculating the largest Lyapunov exponents, we obtain a pseudo-colored map on a $300\times 250$ grid of parameters $(k,c)$ (see Figure 12a). Colors correspond to magnitudes of the largest Lyapunov exponents, where green and blue imply equilibrium, yellow indicates a limit cycle, and red represents a state of chaos. It can be observed from Figure 12a that the dynamical mode of system (1) evolves as k and c change. To more clearly show the evolution of chaos, we fix $k=-0.2$, take $c\in [5,15]$, and plot a vertical line A–B–C, presented in Figure 12a for $A=5,B=10,C=15$. The rich dynamics of the evolution process in the division diagram are shown in Figure 9. Starting with the periodic region A, as c increases, the chaos degenerates through period-doubling bifurcations in line A–B, and the system abruptly changes to one equilibrium through the chaotic status in line B–C. Similarly, we plot a horizontal line, D–E, fix c at 11.2, and take $k\in [-2,-1]$. The system changes from one stable equilibrium to a chaotic state at $k=-1.5$.

Taking initial conditions $(x_0,y_0,z_0)=(1,10,1)$, the numerical results in Figure 12b illustrate that the other initial conditions have an impact on the division diagram, leading to expansion of the regions of stable equilibrium. The original chaotic state regions become an equilibrium state, indicating the existence of a hidden attractor in the corresponding parameters. If the regions remain in a chaotic state, a self-excited attractor might exist under the corresponding parameters, as the regions may become one equilibrium at other initial values.

Inspired by the k–c division diagram, it is clear that the parameter values $k=0$ and $c=0$ at the lower-right corner are dark red, which means the chaos is most complex. As will be discussed shortly, system (1) has a self-excited chaotic attractor at parameter values $(a,b,c,k)=(10,64,0,0)$.

4. Diverse Symbolic Dynamics for Unstable Periodic Orbits

We employ the variational method for the cycle search in system (1) and establish appropriate symbolic dynamics for the found periodic orbits. We first introduce the variational method, which can be effectively used in the calculations. After that, we aim to accurately find the surrounding mode of the orbit in system (1), and develop a universal approach for the symbolic encodings of cycles. We select two sets of parameters, one corresponding to the hidden chaotic attractor, and the other to the self-excited chaotic attractor. The symbolic encoding method based on orbit topology will enable us to analyze periodic orbits by establishing diverse symbolic dynamics.

As shown in Figure 1, the strange attractor of system (1) is composed of numerous unstable periodic orbits. The 3D continuous flow can be transformed to a 2D discrete mapping by an appropriate Poincaré section. The idea is to select a section properly in a high-dimensional phase space, on which a pair of conjugate variables are fixed; then, the information about the motion characteristics can be obtained by observing the intersection points of the motion trajectory and cross section. Figure 13a shows the first return map of system (1) for $(a,b,c,k)=(10,100,11.2,-0.2)$. When we choose a special Poincaré section $z=-9.1503$, the initial values are $[1,1,1],$ where a dense point with a four-branch structure is presented under these parameters, which indicates the necessity to encode all short cycles by symbolic dynamics with four letters. For the parameters $(a,b,c,k)=(10,64,0,0)$, which also correspond to a chaotic state, the first return map with a Poincaré section $z=0$ with the same initial values is shown in Figure 13b. We can see more branches in this case, which means that more symbols are needed to encode the periodic orbits, and demonstrates better complexity in the topological structure of periodic orbits. To the best of our knowledge, investigations of such complex unstable cycles in the chaotic attractor have rarely been reported. The Lyapunov exponents under the parameters $(a,b,c,k)=(10,64,0,0)$ are also calculated, which gives $L{E}_1=1.4456$, $L{E}_2=0.0017$, and $L{E}_3=-11.4473$ (see Figure 14). Correspondingly, the Kaplan–Yorke dimension is $D_{KY}=2.1264$. Compared with the largest Lyapunov exponent, i.e., 0.7457, under the parameters $(a,b,c,k)=(10,100,11.2,-0.2)$, the largest Lyapunov exponent becomes larger, which indicates that the chaotic characteristics of the system are more complex for the parameters $(a,b,c,k)=(10,64,0,0)$.

4.1. Variational Method

Periodic orbits play important roles in physical and engineering applications. It is relatively easy to locate the unstable cycles in low-dimensional chaotic systems in general [54]. When we locate them in a high-dimensional state space, because the topological structure of the dynamical system is difficult to perceive, even if points on the cycle are guessed, the shooting method may fail. This problem can be solved by initializing a complete orbit with similar topology, and making it gradually evolve into a real cycle. This is also the basic idea for the variational method to calculate unstable periodic orbits of dynamical systems. For the calculations of unstable cycles, a discretization equation was derived as [55]

$$\left(\begin{array}{cc}\stackrel{\wedge }{A}& -\stackrel{\wedge }{v}\\ \stackrel{\wedge }{a}& 0\end{array}\right)\left(\genfrac{}{}{0pt}{}{\delta \stackrel{\sim }{x}}{\delta \lambda }\right)=\delta \tau \left(\genfrac{}{}{0pt}{}{\lambda \stackrel{\wedge }{v}-\stackrel{\wedge }{\stackrel{\sim }{v}}}{0}\right),$$

(6)

where $\tau$ is the virtual time related to iteration times. We use $\lambda$ to adjust the period, which has the relationship with the period $T=2\pi \lambda$ when the periodic orbit converges. We want to match the vector fields $\stackrel{\wedge }{v}={(v_1,v_2,\cdots ,v_N)}^t,\stackrel{\wedge }{\stackrel{\sim }{v}}={(\stackrel{\sim }{v_1},\stackrel{\sim }{v_2},\cdots \stackrel{\sim }{v_N})}^t$ everywhere along the loop, and v and $\stackrel{\sim }{v}$ represent the flow velocity and loop velocity vector, respectively. $\stackrel{\wedge }{a}$ is an $Nd$-dimensional row vector that restricts coordinate alterations. $\stackrel{\wedge }{A}=\stackrel{\wedge }{D}-\lambda \mathrm{diag}[A_1,A_2,\cdots ,A_N],$ where $A_{ij}=\frac{\partial v_i}{\partial x_j}$, is the gradient matrix of the velocity field, and the five-point approximation matrix is

$$\stackrel{\wedge }{D}=\frac{N}{24\pi }\left(\begin{array}{ccccccccccc}0& 8& -1& & & & & & & 1& -8\\ -8& 0& 8& -1& & & & & & & 1\\ 1& -8& 0& 8& -1& & & & & & \\ & & & & & ···& & & & & \\ & & & & & & 1& -8& 0& 8& -1\\ -1& & & & & & & 1& -8& 0& 8\\ 8& -1& & & & & & & 1& -8& 0\end{array}\right),$$

(7)

where each matrix element is $d\times d$ dimensional, and blanks are filled with zeros. The $\left[\right(Nd+1)\times (Nd+1\left)\right]$ matrix on the left side of Equation (6) must be inverted to solve for $\delta \stackrel{\sim }{x}$ and $\delta \lambda$, and the banded lower-upper decomposition method for accelerated computing and the Woodbury formula are adopted. In addition, because the virtual time steps are sometimes not significant, we may choose larger time steps for numerical integration, so as to effectively search the real periodic orbit.

To utilize the variational method, as a first step, a loop guess which is fit to the cycle calculations is a prerequisite, and the loop guess can be initialized in many ways in the numerical calculations [55]. For example, we can use a fast Fourier transform of a nearly closed orbital fragment obtained from the numerical integration, keep only the lowest-frequency components, and use a reverse fast Fourier transform back to the phase space, in which emerges a glossy loop guess that can be used to initialize. We can also easily construct the initial loop guess by utilizing the homotopy evolution method [56]. We mention other initialization methods below.

The flexibility of the variational method for cycle searching has been verified by many examples [57,58,59], including conservative systems and low- or high-dimensional dissipative systems. The method can also locate other invariant sets in dynamical systems with proper modification [60,61]. The method cannot only find cycles with fixed parameters, but can be used to study the deformations of cycles when changing some parameters, i.e., to investigate the bifurcation behaviors of a dynamical system [62,63]. Hence, this approach can be used to study the generation or disappearance of periodic orbits and the change of cycle stability.

4.2. Unstable Cycles Embedded in Hidden Chaotic Attractor for $(a,b,c,k)=(10,100,11.2,-0.2)$

The unstable periodic orbits in system (1) when $(a,b,c,k)=(10,100,11.2,-0.2)$ are investigated based on the variational method. In the process, establishing appropriate symbolic dynamics is important for locating all short cycles without missing any [64]. To obtain the topological shape of the periodic orbit to be calculated, we perform numerical simulations, intercepting part of the simple orbital fragment to construct the initial loop guess. Several short periodic orbits with uncomplicated topological structures are found, as shown in Figure 15. Figure 15a shows a periodic orbit that revolves one turn around the left equilibrium $E_1$ with an elliptical shape, which has a relatively small extension in the z-axis with shortest period $T=0.635920.$ We mark it as cycle 0. Figure 15b shows a cycle that rotates once around the right equilibrium $E_2$ with an elliptical shape, and mark it as cycle 1. It can be seen that the two periodic orbits are symmetric to each other. Similarly, we mark the cycle with a wing shape rotating around the fixed point on the left once as cycle 2, as shown in Figure 15c, and its symmetric cycle with reasonably large extension in the z orientation is denoted as cycle 3 (see Figure 15d). The above four cycles can be regarded as the building blocks, and other periodic orbits can be calculated systematically by the symbolic dynamics of four letters.

There are four situations in which an orbit revolves one turn both around the left and right fixed points, and they are the cycles with topological length 2, as listed in Figure 16a–d. The rotationally symmetric property of system (1) implies the exchange symmetry 0 and 1 or 2 and 3 of the symbol sequence. Consequently, it is shown that the symmetry partner of cycle 12 is 03, and they have the same period. Cycles 01 or 23 are conjugated with themselves, so that no other orbit has the same period. Figure 16e–h display four cycles with topological length 3. Utilizing symbolic dynamics, we can calculate the cycles up to any topological length, i.e., we first construct the loop guess of the corresponding symbol sequence, and use the variational technique to verify its existence. Altogether, we found 20 periodic orbits with topological length 3, as listed in

Table 1. Twenty unstable periodic orbits embedded in hidden chaotic attractor of system (1) for $(a,b,c,k)=(10,100,11.2,-0.2)$, showing topological length, itinerary p, period $T_p$, and three coordinates of a point on the periodic orbit.

Length	p	${\mathit{T}}_{\mathit{p}}$	x	y	z
1	0	0.635920	−7.028076	0.430355	1.092913
	1	0.635920	7.028076	−0.430355	1.092913
	2	1.192933	2.544434	12.123766	20.650672
	3	1.192933	−2.544434	−12.123766	20.650672
2	12	1.752388	2.807407	6.313538	12.146665
	03	1.752388	−2.807407	−6.313538	12.146665
	01	1.467965	0.100280	1.271667	−7.073590
	23	2.383824	−14.090307	17.922194	33.086632
3	001	2.174153	−7.214950	4.081725	7.082844
	011	2.174153	7.214950	−4.081725	7.082844
	003	2.361334	−1.162938	−0.669254	−14.983886
	112	2.361334	1.162938	0.669254	−14.983886
	132	2.940945	−2.570910	−0.291669	11.591391
	023	2.940945	2.570910	0.291669	11.591391
	021	2.554559	−0.016908	−0.016629	−31.681837
	013	2.554559	0.016908	0.016629	−31.681837
	033	2.946229	0.291076	0.540818	−60.155264
	122	2.946229	−0.291076	−0.540818	−60.155264
	223	3.954898	−5.509519	−11.111451076	−71.906347
	233	3.954898	5.509519	11.111451076	−71.906347

.

4.3. Unstable Periodic Orbits Embedded in Self-Excited Chaotic Attractor for $(a,b,c,k)=(10,64,0,0)$

If we make the third and fourth parameters in the chosen set $(a,b,c,k)$ zero, according to the Routh–Hurwitz stability criterion, $a^2{b}^2<0$ must be satisfied, so there is no solution, which means that there will be no hidden chaotic attractors in the system. Only when c and k are not zero is it possible to satisfy the Routh–Hurwitz stability criterion and a hidden chaotic attractor with stable equilibrium points can exist. When we take another set of parameters, $(a,b,c,k)=(10,64,0,0)$, system (1) becomes a dynamical system with only five terms and also exhibits the existence of a chaotic state. As with some simple chaotic flows, namely the Sprott system, listed in Ref. [15], system (1) is simpler, but it has more complex dynamics, which is worthy of further research. For these parameter values, the two equilibrium points yield $E_1=(-\sqrt{b},-\sqrt{b},0)$ and $E_2=(\sqrt{b},\sqrt{b},0),$ and three eigenvalues of the fixed points $E_{1,2}$ are ${\lambda }_{1,2}=1.4034\pm 9.8983i$ and ${\lambda }_3=-12.8068$. Since $E_{1,2}$ are two saddle-foci, system (1) has a self-excited chaotic attractor under current parameters.

We locate the unstable periodic orbits via the variational method, which brings great convenience. Here, we numerically integrate Equation (1) and extract approximate closed trajectories with different shapes, then artificially connect them so as to initialize the search. Figure 17 displays our calculated results for some short periodic orbits. We can also record the cycle swirling around the left fixed point $E_1$ once with a wing shape in Figure 17a by cycle 2, and its symmetric partner in Figure 17b with cycle 3. However, we did not find cycle 0 and cycle 1 to exist. The cycles in Figure 17c,d both have a knot, which indicates that they have a self-linking number of 1, which can conveniently be calculated [65]. We mark them as cycle 4 and cycle 5, respectively. Noting that they are symmetric to each other, the commutative symmetry 4 and 5 of the symbol sequence is satisfied.

We can use the above cycles as building blocks to find more complicated cycles, and the symbolic dynamics is established successfully. Figure 18 shows part of the found cycles together with equilibria $E_1$ and $E_2$. According to the symbolic dynamics, we find 41 unstable cycles with topological lengths up to 3, and sort them in

Table 2. Forty-one unstable periodic orbits embedded in self-excited chaotic attractor of system (1) for $(a,b,c,k)=(10,64,0,0)$.

Length	p	${\mathit{T}}_{\mathit{p}}$	Self-Linking	Length	p	${\mathit{T}}_{\mathit{p}}$	Self-Linking	p	${\mathit{T}}_{\mathit{p}}$	Self-Linking
1	2	1.016946	0	3	223	2.994130	0	031	2.447450	2
	3	1.016946	0		233	2.994130	0	012	2.447450	2
2	01	1.358438	1		033	2.609712	2	132	2.505368	0
	23	1.965825	1		122	2.609712	2	023	2.505368	0
	12	1.587528	1		021	2.323226	0
	03	1.587528	1		013	2.323226	0
1	4	1.312552	1		445	4.235720	3	354	3.955079	2
	5	1.312552	1		455	4.235720	3	234	3.263831	1
2	24	2.289914	2		344	3.667897	1	325	3.263831	1
	25	2.354458	0		255	3.667897	1	225	3.367249	1
	34	2.354458	0		335	3.312270	1	334	3.367249	1
	35	2.289914	2		224	3.312270	1	254	3.606269	3
	45	2.642183	1		244	3.600833	3	345	3.606269	3
3	235	3.349139	1		355	3.600833	3
	324	3.349139	1		245	3.955079	2

. A total of 14 cycles are pruned, e.g., cycles 02, 002, and 123. Compared with the two sets of parameters, for the same periodic orbits, the cycles embedded in the hidden chaotic attractor have longer periods than those embedded in the self-excited chaotic attractor. The unstable cycles of system (1), as discussed under current parameters, must invoke symbolic dynamics for six letters, which is usually complicated. The topological classification approach used here indicates its flexibility. Additionally, although the symbolic dynamics of six letters can produce many symbol sequences within the topological length of 3, it is found that the 2 and 3 building blocks can be combined with all the other building blocks, while the 0 and 1 building blocks cannot be combined with the 4 or 5 building blocks. These behaviors are surely unusual. Therefore, the number of cycles actually allowed by the symbol sequence is greatly reduced. Whether this empirical pruning rule is applicable to longer periodic orbits is an open problem worthy of further investigation.

Regarding the system’s overall dynamical complexity from the two sets of parameters chosen in the study, we can draw the following conclusions:

(1) The proposed system (1) with two parameters has more complex dynamics than the system with four parameters.

(2) The system with a self-excited attractor has more complex dynamics than the system with a hidden attractor.

(3) The system with periodic orbits containing building blocks of self-linking number 1 has more complex dynamics than that containing building blocks of self-linking number 0.

5. Circuit Simulation

We discuss the circuit implementation of system (1) to confirm the realizability of the mathematical model. Because all the values of state variables $(x,y,z)$ in system (1) are out of the dynamic range, they should be scaled down to avoid problems during simulation. We set the amplitude scaling factor to 10, where $X=\frac{1}{10}x$, $Y=\frac{1}{10}y$, and $Z=\frac{1}{10}z$. The time scale factor is set to ${\tau }_0=\frac{1}{R_0{C}_0}=2500$ to better match the system, a new time variable $\tau$ is defined instead of t, and $t={\tau }_0\tau .$ As a result, system (1) after scale transformation is described as

$$\begin{array}{ccc}\hfill R_0{C}_0\stackrel{·}{X}& =& a(Y-X)+10kXZ\hfill \\ \hfill R_0{C}_0\stackrel{·}{Y}& =& -cY-10XZ\hfill \\ \hfill R_0{C}_0\stackrel{·}{Z}& =& -\frac{b}{10}+10XY,\hfill \end{array}$$

(8)

where $a=10,b=100,c=11.2,$ and $k=-0.2$.

The proposed circuit design is depicted in Figure 19, in which X, Y, and Z are the voltages at the outputs of operational amplifiers U2, U5, and U8, respectively. The circuit consists of nine AD811AN operational amplifiers, whose supply voltage is ±18 V; three multipliers with an output coefficient of 0.1; three capacitors; and 19 resistors. Based on Kirchhoff’s law, we can get

$$\begin{array}{ccc}\hfill \stackrel{·}{X}& =& \frac{R_4}{R_2{R}_5{C}_1}Y-\frac{R_4}{R_1{R}_5{C}_1}X-\frac{R_4}{R_3{R}_5{C}_1}0.1XZ\hfill \\ \hfill \stackrel{·}{Y}& =& -\frac{R_{10}}{R_8{R}_{11}{C}_2}Y-\frac{R_{10}}{R_9{R}_{11}{C}_2}0.1XZ\hfill \\ \hfill \stackrel{·}{Z}& =& \frac{R_{16}}{R_{15}{R}_{17}{C}_3}{V}_1+\frac{R_{16}}{R_{14}{R}_{17}{C}_3}0.1XY.\hfill \end{array}$$

(9)

Comparing Equation (8) with Equation (9), we select all the capacitors $C_i$ = 40 nF $(i=1,2,3)$ and $V1$ = −1 V. The resistors $R_9=R_{14}$ = 1 kΩ, $R_3$ = 5 kΩ, $R_8$ = 8.93 kΩ, $R_i$ = 10 kΩ $(i=1,2,5,6,7,11,12,13,15,17,18,19)$, and $R_i$ = 100 kΩ $(i=4,10,16).$ We used NI Multisim 14.0 to simulate the circuit, as shown in Figure 20. It can be seen that the circuit well emulates the proposed system, which is in good agreement with the numerical results in Figure 1. Therefore, we can conclude that system (1) can be realized in physical experiments.

6. Conclusions and Discussion

We constructed a new 3D autonomous chaotic system with coexisting self-excited and hidden attractors, in which the generation of different types of attractors depends on its parameters. The complex dynamics of the system were analyzed by different tools, and it was proved to be chaotic in the sense of having a fractional Kaplan–Yorke dimension, a phase portrait of a strange attractor, and a period-doubling route to chaos. Moreover, an applicable generic procedure for topological classification of unstable cycles in the proposed system was addressed. Guessing an entire orbit, we utilized the variational method for the calculation of cycles, and the initial conjecture loop could be gradually evolved into a real cycle. Diverse symbolic dynamics based on orbital topology was successfully established in the phase space, including four and six letters, corresponding to hidden and self-excited chaotic attractors. Periodic orbits up to certain topological lengths were found accordingly, which indicates the utility of the topological classification approach in the periodic orbit taxonomy. A Multisim circuit simulation of the system was implemented to further verify the mathematical model.

The symbolic encoding method employed here could also be applied to discrete dynamical systems, such as the memristive Rulkov neuron model [66], discrete memristor hyperchaotic maps [67], 2D memristive hyperchaotic maps [68], and 2D sine map [69]. Quotienting symmetries of a given dynamical system prior to the symbolic dynamics analysis is an attractive research direction. Symmetry reduction could not only reduce the multiple-letter symbolic encodings of periodic orbits to a single letter, but could visualize self-linking in the symmetry-reduced state space, which requires further investigation. The new system still contains rich and complex dynamic behavior, and its topology requires comprehensive and deep exploration. Moreover, as the system proposed in Ref. [51], the newly proposed system (1) is a mathematical model at present and does not correspond to any physical phenomena. It is hoped that more detailed theoretical analysis and application investigations will be carried out in the future.