Dynamic Analysis and Circuit Implementation of a New 4D Lorenz-Type Hyperchaotic System

Abstract

This paper attempts to further extend the results of dynamical analysis carried out on a recent 4D Lorenz-type hyperchaotic system while exploring new analytical results concerns its local and global dynamics. In particular, the equilibrium points of the system along with solution’s continuous dependence on initial conditions are examined. Then, a detailed $Z_{\mathrm{2}}$ symmetrical Bogdanov-Takens bifurcation analysis of the hyperchaotic system is performed. Also, the possible first integrals and global invariant surfaces which exist in system’s phase space are analytically found. Theoretical results reveal the rich dynamics and the complexity of system behavior. Finally, numerical simulations and a proposed circuit implementation of the hyperchaotic system are provided to validate the present analytical study of the system.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Nonlinear dynamics analysis of various phenomena and systems of physics, engineering, biology, chemistry, economy, and industry has attracted a great interest among scientists and considered a very active area of research from 1960s in the last century till now [1–4]. There are two key reasons which interpret this great interest. The first one is that the dynamical systems tools help scientists better comprehend and analyze the varieties of nonlinear characteristics and new phenomena exhibited by systems from different disciplines. In particular, the tools of dynamical systems such as the applied bifurcation theories are successfully employed to investigate the qualitative behaviors of nonlinear systems [5–7]. This includes investigation of equilibrium points and their stability, creation, destruction and stability of periodic orbits, quasiperiodic behavior, homoclinic orbits, creation or destruction of chaotic attractors, and chaos control and synchronization. The second reason is that engineers and scientists can utilize some of the fascinating features of nonlinear dynamical systems in wide range of interesting applications.

The spellbinding chaotic dynamics, as an example, is recognized by high sensitivity to initial conditions and positive Lyapunov exponents. The generation of chaos for practical applications can be achieved by exploiting nonlinear electronic circuits, nonlinearities of laser systems with feedbacks, or via digital platforms such as FPGAs and DSPs. The noise-like behavior of chaotic systems, their wide spectrum, and the possibility of attaining chaos synchronization between two chaotic systems render them essential for cutting edge applications related to cryptography and robust physical-layer secure communication systems and ultra-fast physical random bits generation [4–11]. The other applications of dynamical systems methods and chaos theory include financial systems, mathematical biology, nonlinear circuits, nonlinear mechanical systems, plasma physics, chaos control, efficient image encryption, neuroscience research, and geophysics [12–33].

The more complex hyperchaotic system possesses at least two positive Lyapunov exponents and it has a phase space of dimension at least four. Clearly, the hyperchaotic systems have more randomness and higher unpredictability than simple chaotic systems. Therefore, hyperchaos is more preferred than simple chaos and its applications have recently become a central topic in research including chaos-based secure communications, image encryption, and cryptography

In the last two decades, some interesting high-dimensional hyperchaotic systems, in science and engineering, have been explored and their dynamics have been extensively investigated [18, 24, 27, 28, 34, 35]. In fact, it is of great importance from theoretical and practical aspects to explain complicated phenomena and internal structural characteristics of hyperchaotic systems. This research line focuses on applying codimension two or three bifurcation analysis to system under investigation. For example, analysis of Bogdanov-Takens bifurcation, degenerate Hopf bifurcation, and Heteroclinic and Homoclinic bifurcations can be undertaken. Also, it is crucial to determine whether there exists any set of parameters values for which the considered system is integrable and find the corresponding invariant surfaces if exist. As it is practically unattainable goal to apply exhaustive numerical investigation to acquire values of parameters where the hyperchaotic system is integrable, one needs a powerful analytical method that enables achieving this goal easily.

Recently, theoretical analysis, based on Fishing Principle, is applied to study some global features of a new four-dimensional Lorenz-type hyperchaotic system [36]. More specifically, conditions for existence of homoclinic orbits are obtained. This work aims at extending the aforementioned work and exploring other aspects of complicated dynamical behaviors of the 4D Lorenz-type hyperchaotic system. Analytical bifurcation structure of the model which includes $Z_{\mathrm{2}}$ symmetrical Bogdanov-Takens, Pitchfork, Andronov-Hopf bifurcation and homoclinic bifurcation are obtained. Existence, uniqueness, and continuous dependence on initial conditions for the solution of 4D hyperchaotic system are all investigated. It is important to ask whether there is any set of parameter values in which system dynamics are regular and the studied system is integrable. The different integrable cases of the system are studied and the closed-form expressions for invariant surfaces corresponding to first integrals of the system are also obtained. Moreover, the hyperchaotic system has been simulated using a proposed electronic circuit realization and numerical simulations are performed confirming the new results of theoretical analysis. Note that attractors found in dynamical system are classified as being self-excited or hidden attractors [37, 38]. The key difference between them is that the self-excited attractor has a basin of attraction excited from unstable equilibria. The hidden attractor, on the other hand, has a basin of attraction with no intersection with neighborhoods of dynamical systems equilibrium points. In this work, we focus on self-excited attractors exist in the hyperchaotic system.

The rest of the paper is structured as follows: In Section 2, we introduce the 4D Lorenz-type system and discuss the equilibrium points’ existence. A sufficient condition for continuous dependence on initial condition is determined. Phase portraits, bifurcation diagrams, and Lyapunov characteristics spectrum are obtained. In Section 3, the analysis of some possible codimension two bifurcations is performed. It is shown that the 4D Lorenz system undergoes Bogdanov-Takens bifurcation, Andronov-Hopf bifurcation, Pitchfork bifurcation, and homoclinic bifurcation. The integrability analysis of system is investigated in Section 4. A practical application to engineering will be realized by an electronic circuit in Section 5. Finally, Section 6 concludes the paper.

2. The 4D Lorenz-Type Hyperchaotic System

The following 4D hyperchaotic system of Lorenz type was presented in [36]$$\begin{array}{c}\text{(1)}& \dot{x}=a\left(y-x\right),\dot{y}=cx-dy-xz,\\ \dot{z}=-bz+xy+w,\\ \dot{w}=-rw+kz,\end{array}$$where state variables of the system are denoted by $x,y,z$, and $w$, the parameters of the system are represented by $a,b,c,d,r$, and $k$, and the dot above state variables refers to time derivative of state variables.

The following subsections examine the main properties of system (1) and provide elementary dynamical analysis of the model.

2.1. Equilibrium Points of the System

The fixed points of system (1) can be obtained when $\dot{x}=\dot{y}=\dot{z}=\dot{w}=\mathrm{0}$, i.e., through solving$$\begin{array}{c}\text{(2)}& a\left(y-x\right)=\mathrm{0},cx-dy-xz=\mathrm{0},\\ -bz+xy+w=\mathrm{0},\\ -rw+kz=\mathrm{0}.\end{array}$$Therefore, system (1) has the following equilibrium solutions:

(I) when $(d-c)(k-br)/r\le \mathrm{0}$, the fixed point is given by $$\begin{array}{}\text{(3)}& x=y=z=w=\mathrm{0},\end{array}$$

(II) when $(d-c)(k-br)/r>\mathrm{0}$, the fixed points are$$\begin{array}{c}\text{(4)}& x=y=z=w=\mathrm{0},x=y=\pm \sqrt{\frac{\left(d-c\right)\left(k-br\right)}{r}},\\ z=c-d,\\ w=\frac{k\left(c-d\right)}{r}.\end{array}$$

2.2. Existence and Uniqueness of the Solution

The hyperchaotic system (1) can be put in the form$$\begin{array}{c}\text{(5)}& \dot{\mathbf{X}}\left(\mathbf{t}\right)=\mathbf{\Phi }\left(\mathbf{X}\left(\mathbf{t}\right)\right),t\in \left(\mathrm{0},T\right],\end{array}$$where$$\begin{array}{c}\text{(6)}& \mathbf{X}=\left[\begin{array}{c}{x}_{\mathrm{1}}\\ x_{\mathrm{2}}\\ x_{\mathrm{3}}\\ x_{\mathrm{4}}\end{array}\right],{\mathbf{X}}_{\mathrm{0}}=\left[\begin{array}{c}{x}_{\mathrm{01}}\\ x_{\mathrm{02}}\\ x_{\mathrm{03}}\\ x_{\mathrm{04}}\end{array}\right],\\ \mathbf{\Phi }\left(\mathbf{X}\right)=\left[\begin{array}{c}a\left(x_{\mathrm{2}}-x_{\mathrm{1}}\right),\\ c{x}_{\mathrm{1}}-d{x}_{\mathrm{2}}-x_{\mathrm{1}}{x}_{\mathrm{3}}\\ -b{x}_{\mathrm{3}}+x_{\mathrm{1}}{x}_{\mathrm{2}}+x_{\mathrm{4}}\\ -r{x}_{\mathrm{4}}+k{x}_{\mathrm{3}}\end{array}\right],\end{array}$$and the initial conditions are given by$$\begin{array}{}\text{(7)}& \mathbf{X}\left(\mathbf{0}\right)={\mathbf{X}}_{\mathrm{0}}.\end{array}$$

For the class of continuous functions $F(t)\in C[\mathrm{0},T]$, we use the following norm in subsequent analysis:$$\begin{array}{}\text{(8)}& ‖F‖=\underset{t\in \left[\mathrm{0},T\right]}{\mathrm{s}\mathrm{u}\mathrm{p}}\left|F\left(t\right)\right|,\end{array}$$while the matrix $K=[k_{ij}[t]]$ of continuous functions employs the norm$$\begin{array}{}\text{(9)}& ‖K‖={\displaystyle \sum }_{i,j}\underset{t\in \left[\mathrm{0},T\right]}{\mathrm{s}\mathrm{u}\mathrm{p}}\left|k_{ij}\left[t\right]\right|.\end{array}$$

It is obvious that 4D system (1) is dissipative if $\nabla .\mathbf{F}(\mathbf{X})=-(a+b+d+r)<\mathrm{0}.$ Now, solution of the system is examined in specific region $\Gamma \times J$ where $J=[\mathrm{0},T]$ and$$\begin{array}{}\text{(10)}& \Gamma =\left\{\left(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}},x_{\mathrm{4}}\right):\max \left\{\left|x\right|,\left|x_{\mathrm{2}}\right|,\left|x_{\mathrm{3}}\right|\hspace{0.17em}\hspace{0.17em}\mathrm{a}\mathrm{n}\mathrm{d}\hspace{0.17em}\hspace{0.17em}\left|x_{\mathrm{4}}\right|\right\}\le M\right\},\hspace{1em}M>\mathrm{0}.\end{array}$$Parameter $M$ is utilized to lay a boundary for the phase space region where existence and uniqueness of the solution are investigated.

The solution of (5) and (7) can be represented by$$\begin{array}{}\text{(11)}& \mathbf{X}\left(\mathbf{t}\right)={\mathbf{X}}_{\mathrm{0}}+{{\displaystyle \int }}_{\mathrm{0}}^t\mathbf{F}\left(\mathbf{X}\left(\mathbf{s}\right)\right)\mathbf{d}\mathbf{s}.\end{array}$$The equivalence of the integral equation (11) and system (5)-(7) is obvious. Now, denoting the right hand side of (11) by $\mathbf{H}(\mathbf{X})$, then for ${\mathbf{X}}_{\mathrm{1}}=\left[\begin{array}{c}{x}_{\mathrm{11}}\\ x_{\mathrm{12}}\\ x_{\mathrm{13}}\\ x_{\mathrm{14}}\end{array}\right]$ and ${\mathbf{X}}_{\mathrm{2}}=\left[\begin{array}{c}{x}_{\mathrm{21}}\\ x_{\mathrm{22}}\\ x_{\mathrm{23}}\\ x_{\mathrm{24}}\end{array}\right]$ we get$$\begin{array}{}\text{(12)}& \mathbf{H}\left({\mathbf{X}}_{\mathrm{1}}\right)-\mathbf{H}\left({\mathbf{X}}_{\mathrm{2}}\right)={{\displaystyle \int }}_{\mathrm{0}}^t\left(\mathbf{\Phi }\left({\mathbf{X}}_{\mathrm{1}}\left(\mathbf{s}\right)\right)-\mathbf{\Phi }\left({\mathbf{X}}_{\mathrm{2}}\left(\mathbf{s}\right)\right)\right)\mathbf{d}\mathbf{s},\end{array}$$and therefore $$\begin{array}{}\text{(13)}& \left|\mathbf{H}\left({\mathbf{X}}_{\mathrm{1}}\right)-\mathbf{H}\left({\mathbf{X}}_{\mathrm{2}}\right)\right|\le {{\displaystyle \int }}_{\mathrm{0}}^t\left|\left(\mathbf{\Phi }\left({\mathbf{X}}_{\mathrm{1}}\left(\mathbf{s}\right)\right)-\mathbf{\Phi }\left({\mathbf{X}}_{\mathrm{2}}\left(\mathbf{s}\right)\right)\right)\right|\mathbf{d}\mathbf{s}.\end{array}$$Then, we obtain $$\begin{array}{}\text{(14)}& ‖\mathbf{H}\left({\mathbf{X}}_{\mathrm{1}}\right)-\mathbf{H}\left({\mathbf{X}}_{\mathrm{2}}\right)‖\le T\max \left\{\left|a\right|+\left|c\right|+\mathrm{2}M,\left|a\right|+\left|d\right|+M,\left|b\right|+\left|k\right|+M,\mathrm{1}+\left|r\right|\right\}‖{\mathbf{X}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{2}}‖\le K‖{\mathbf{X}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{2}}‖,\end{array}$$where$$\begin{array}{}\text{(15)}& K=T\max \left\{\left|a\right|+\left|c\right|+\mathrm{2}M,\left|a\right|+\left|d\right|+M,\left|b\right|+\left|k\right|+M,\mathrm{1}+\left|r\right|\right\}>\mathrm{0}.\end{array}$$Thus, $\mathbf{X}=\mathbf{H}(\mathbf{X})$, for $\mathrm{0}<K<\mathrm{1}$ as sufficient condition, is a contraction mapping.

Theorem 1.

Assume that $\mathrm{0}<T\max \left\{\left|a\right|+\left|c\right|+\mathrm{2}M,\left|a\right|+\left|d\right|+M,\left|b\right|+\left|k\right|+M,\mathrm{1}+\left|r\right|\right\}<\mathrm{1}$, then a unique solution of hyperchaotic system (1), starting from initial condition $\mathbf{X}(\mathbf{0})={\mathbf{X}}_{\mathrm{0}}$ in the region $\Omega \times J$, exists.

2.3. Continuous Dependence on Initial Conditions

The continuous dependence on initial conditions means that solution trajectories of the system which start close to each other still close to each other with evolution of time. This property is contrary to sensitive dependence on initial conditions which specifies chaotic dynamics. The goal of the next analysis is to find the particular parameters and time range where continuous dependence on state variables initial conditions is persevered; i.e., system (1) does not exhibit chaotic dynamics.

Assume that there are two points of initial conditions of system (5)-(7); namely, ${\mathbf{X}}_{\mathrm{01}}$ and ${\mathbf{X}}_{\mathrm{02}}$ satisfy $$\begin{array}{}\text{(16)}& ‖{\mathbf{X}}_{\mathrm{01}}-{\mathbf{X}}_{\mathrm{02}}‖\le \delta .\end{array}$$First, suppose that the condition of Theorem 1 holds. Thus $$\begin{array}{c}\text{(17)}& {\mathbf{X}}_{\mathrm{1}}\left(t\right)={\mathbf{X}}_{\mathrm{01}}+{{\displaystyle \int }}_{\mathrm{0}}^t\mathbf{H}\left({\mathbf{X}}_{\mathrm{1}}\left(\mathbf{s}\right)\right)\mathbf{d}\mathbf{s},{\mathbf{X}}_{\mathrm{2}}\left(t\right)={\mathbf{X}}_{\mathrm{02}}+{{\displaystyle \int }}_{\mathrm{0}}^t\mathbf{H}\left({\mathbf{X}}_{\mathrm{2}}\left(\mathbf{s}\right)\right)\mathbf{d}\mathbf{s},\end{array}$$and also we obtain $$\begin{array}{}\text{(18)}& ‖{\mathbf{X}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{2}}‖\le ‖{\mathbf{X}}_{\mathrm{01}}-{\mathbf{X}}_{\mathrm{02}}‖+K‖{\mathbf{X}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{2}}‖,\end{array}$$and $$\begin{array}{}\text{(19)}& \left(\mathrm{1}-K\right)‖{\mathbf{X}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{2}}‖\le ‖{\mathbf{X}}_{\mathrm{01}}-{\mathbf{X}}_{\mathrm{02}}‖,\end{array}$$where $$\begin{array}{}\text{(20)}& \mathrm{0}<K<\mathrm{1}\end{array}$$is defined by (15). Finally, we get $$\begin{array}{}\text{(21)}& ‖{\mathbf{X}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{2}}‖\le ϵ,\end{array}$$where $ϵ=\delta /(\mathrm{1}-K)$. Thus, we can formulate the following theorem.

Theorem 2.

The solution of hyperchaotic system (1) exhibits continuous dependence on initial conditions if system (1) satisfies the condition of Theorem 1. More specifically, $\forall$$ϵ>\mathrm{0}$$\exists$$\delta (ϵ)=(\mathrm{1}-K)ϵ>\mathrm{0}$ such that whenever $‖{\mathbf{X}}_{\mathrm{01}}-{\mathbf{X}}_{\mathrm{02}}‖\le \delta$ it implies that $‖{\mathbf{X}}_{\mathrm{1}}-{\mathbf{X}}_{\mathrm{2}}‖\le ϵ$.

3. $Z_2$ Symmetry Bogdanov-Takens Bifurcation

Now, we study the case where $c=c^{\ast }=-a$ and $d=d^{\ast }=-a$ that implies Jacobian matrix has two real zero eigenvalues, ${\lambda }_{\mathrm{1,2}}=\mathrm{0}$, and two negative eigenvalues if either $b^{\mathrm{2}}-\mathrm{2}br+\mathrm{4}k+r^{\mathrm{2}}<b+r$ or $b^{\mathrm{2}}-\mathrm{2}br+\mathrm{4}k+r^{\mathrm{2}}<\mathrm{0}$ evaluated at the origin. The following coordinates’ transformation is applied to (1) in order to put the system in standard form$$\begin{array}{}\text{(22)}& \left(\begin{array}{c}{x}_{\mathrm{1}}\\ x_{\mathrm{2}}\\ x_{\mathrm{3}}\\ x_{\mathrm{4}}\end{array}\right)=\left(\begin{array}{cccc}\mathrm{1}& -\frac{\mathrm{1}}{a}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \frac{-b+r-\varpi }{\mathrm{2}k}& \frac{-b+r+\sqrt{b^{\mathrm{2}}-\mathrm{2}rb+r^{\mathrm{2}}+\mathrm{4}k}}{\mathrm{2}k}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{1}\end{array}\right)\left(\begin{array}{c}{u}_{\mathrm{1}}\\ u_{\mathrm{2}}\\ u_{\mathrm{3}}\\ u_{\mathrm{4}}\end{array}\right),\end{array}$$which yields$$\begin{array}{}\text{(23)}& \left(\begin{array}{c}{\dot{u}}_{\mathrm{1}}\\ {\dot{u}}_{\mathrm{2}}\\ {\dot{u}}_{\mathrm{3}}\\ {\dot{u}}_{\mathrm{4}}\end{array}\right)=\left(\begin{array}{c}{ϝ}_{\mathrm{1}}\\ {ϝ}_{\mathrm{2}}\\ {ϝ}_{\mathrm{3}}\\ {ϝ}_{\mathrm{4}}\end{array}\right),\end{array}$$where$$\begin{array}{c}\text{(24)}& {ϝ}_{\mathrm{1}}=\frac{y_{\mathrm{4}}{y}_{\mathrm{2}}\varpi }{\mathrm{2}ak}-\frac{y_{\mathrm{3}}{y}_{\mathrm{2}}\varpi }{\mathrm{2}ak}-\frac{b{y}_{\mathrm{3}}{y}_{\mathrm{2}}}{\mathrm{2}ak}-\frac{b{y}_{\mathrm{4}}{y}_{\mathrm{2}}}{\mathrm{2}ak}+\frac{r{y}_{\mathrm{3}}{y}_{\mathrm{2}}}{\mathrm{2}ak}+\frac{r{y}_{\mathrm{4}}{y}_{\mathrm{2}}}{\mathrm{2}ak}+\frac{y_{\mathrm{1}}{y}_{\mathrm{3}}\varpi }{\mathrm{2}k}-\frac{y_{\mathrm{1}}{y}_{\mathrm{4}}\varpi }{\mathrm{2}k}+\frac{b{y}_{\mathrm{1}}{y}_{\mathrm{3}}}{\mathrm{2}k}+\frac{b{y}_{\mathrm{1}}{y}_{\mathrm{4}}}{\mathrm{2}k}-\frac{r{y}_{\mathrm{1}}{y}_{\mathrm{3}}}{\mathrm{2}k}-\frac{r{y}_{\mathrm{1}}{y}_{\mathrm{4}}}{\mathrm{2}k}+y_{\mathrm{2}},{ϝ}_{\mathrm{2}}=\frac{a{y}_{\mathrm{1}}{y}_{\mathrm{3}}\varpi }{\mathrm{2}k}-\frac{a{y}_{\mathrm{1}}{y}_{\mathrm{4}}\varpi }{\mathrm{2}k}+\frac{ab{y}_{\mathrm{1}}{y}_{\mathrm{3}}}{\mathrm{2}k}+\frac{ab{y}_{\mathrm{1}}{y}_{\mathrm{4}}}{\mathrm{2}k}-\frac{ar{y}_{\mathrm{1}}{y}_{\mathrm{3}}}{\mathrm{2}k}-\frac{ar{y}_{\mathrm{1}}{y}_{\mathrm{4}}}{\mathrm{2}k}-\frac{y_{\mathrm{2}}{y}_{\mathrm{3}}\varpi }{\mathrm{2}k}+\frac{y_{\mathrm{2}}{y}_{\mathrm{4}}\varpi }{\mathrm{2}k}-\frac{b{y}_{\mathrm{2}}{y}_{\mathrm{3}}}{\mathrm{2}k}-\frac{b{y}_{\mathrm{2}}{y}_{\mathrm{4}}}{\mathrm{2}k}+\frac{r{y}_{\mathrm{2}}{y}_{\mathrm{3}}}{\mathrm{2}k}+\frac{r{y}_{\mathrm{2}}{y}_{\mathrm{4}}}{\mathrm{2}k},\\ {ϝ}_{\mathrm{3}}=\frac{k{y}_{\mathrm{1}}{y}_{\mathrm{2}}}{a\varpi }-\frac{b^{\mathrm{2}}{y}_{\mathrm{3}}}{\mathrm{4}\varpi }-\frac{b^{\mathrm{2}}{y}_{\mathrm{4}}}{\mathrm{4}\varpi }+\frac{br{y}_{\mathrm{3}}}{\mathrm{2}\varpi }+\frac{br{y}_{\mathrm{4}}}{\mathrm{2}\varpi }-\frac{\mathrm{1}}{\mathrm{4}}{y}_{\mathrm{3}}\varpi +\frac{\mathrm{1}}{\mathrm{4}}{y}_{\mathrm{4}}\varpi -\frac{k{y}_{\mathrm{1}}^{\mathrm{2}}}{\varpi }-\frac{k{y}_{\mathrm{3}}}{\varpi }-\frac{k{y}_{\mathrm{4}}}{\varpi }-\frac{r^{\mathrm{2}}{y}_{\mathrm{3}}}{\mathrm{4}\varpi }-\frac{r^{\mathrm{2}}{y}_{\mathrm{4}}}{\mathrm{4}\varpi }-\frac{b{y}_{\mathrm{3}}}{\mathrm{2}}-\frac{r{y}_{\mathrm{3}}}{\mathrm{2}},\\ {ϝ}_{\mathrm{4}}=-\frac{k{y}_{\mathrm{1}}{y}_{\mathrm{2}}}{a\varpi }+\frac{b^{\mathrm{2}}{y}_{\mathrm{3}}}{\mathrm{4}\varpi }+\frac{b^{\mathrm{2}}{y}_{\mathrm{4}}}{\mathrm{4}\varpi }-\frac{br{y}_{\mathrm{3}}}{\mathrm{2}\varpi }-\frac{br{y}_{\mathrm{4}}}{\mathrm{2}\varpi }+\frac{k{y}_{\mathrm{1}}^{\mathrm{2}}}{\varpi }-\frac{\mathrm{1}}{\mathrm{4}}{y}_{\mathrm{3}}\varpi +\frac{k{y}_{\mathrm{3}}}{\varpi }+\frac{r^{\mathrm{2}}{y}_{\mathrm{3}}}{\mathrm{4}\varpi }+\frac{\mathrm{1}}{\mathrm{4}}{y}_{\mathrm{4}}\varpi +\frac{k{y}_{\mathrm{4}}}{\varpi }+\frac{r^{\mathrm{2}}{y}_{\mathrm{4}}}{\mathrm{4}\varpi }-\frac{b{y}_{\mathrm{4}}}{\mathrm{2}}-\frac{r{y}_{\mathrm{4}}}{\mathrm{2}},\end{array}$$such that $\varpi =\sqrt{b^{\mathrm{2}}-\mathrm{2}br+\mathrm{4}k+r^{\mathrm{2}}}.$

The center manifold is assumed in the form of second-order polynomial for sufficiently small $u_{\mathrm{1}}$ and $u_{\mathrm{2}}$ as$$\begin{array}{}\text{(25a)}& u_{\mathrm{3}}={\alpha }_{\mathrm{20}}{u}_{\mathrm{1}}^{\mathrm{2}}+{\alpha }_{\mathrm{11}}{u}_{\mathrm{1}}{u}_{\mathrm{2}}+{\alpha }_{\mathrm{02}}{u}_{\mathrm{2}}^{\mathrm{2}},\text{(25b)}& u_{\mathrm{4}}={\beta }_{\mathrm{20}}{u}_{\mathrm{1}}^{\mathrm{2}}+{\beta }_{\mathrm{11}}{u}_{\mathrm{1}}{u}_{\mathrm{2}}+{\beta }_{\mathrm{02}}{u}_{\mathrm{2}}^{\mathrm{2}}.\end{array}$$then substituting from ((25a)-(25b)) and the first two equations of (23) into the last two equation of (23). By comparing the coefficient of $u_{\mathrm{2}}^i{u}_{\mathrm{3}}^j$, $i+j=\mathrm{2}$ in both sides of last two equations of (23) after substitution in ((25a)-(25b)), the values of center manifold coefficients are obtained as follows: $$\begin{array}{c}\text{(26)}& {\alpha }_{\mathrm{20}}=-\frac{\mathrm{2}k}{b^{\mathrm{2}}+b\left(\varpi -\mathrm{2}r\right)+r\left(\varpi +r\right)+\mathrm{4}k},{\alpha }_{\mathrm{11}}=\frac{\mathrm{2}k\left(\mathrm{4}a\varpi +b^{\mathrm{2}}+b\left(\varpi -\mathrm{2}r\right)+r\left(\varpi +r\right)+\mathrm{4}k\right)}{a{\left(b^{\mathrm{2}}+b\left(\varpi -\mathrm{2}r\right)+r\left(\varpi +r\right)+\mathrm{4}k\right)}^{\mathrm{2}}},\\ {\alpha }_{\mathrm{02}}=-\frac{\mathrm{4}k\left({\left(b-r\right)}^{\mathrm{2}}+\mathrm{4}k\right)\left(\mathrm{4}a+\varpi +b+r\right)}{a{\left(b^{\mathrm{2}}+b\left(\sqrt{{\left(b-r\right)}^{\mathrm{2}}+\mathrm{4}k}-\mathrm{2}r\right)+r\left(\varpi +r\right)+\mathrm{4}k\right)}^{\mathrm{3}}},\\ {\beta }_{\mathrm{20}}=-\frac{k\left(\varpi +b+r\right)}{\mathrm{2}\varpi \left(k-br\right)},\\ {\beta }_{\mathrm{11}}=\frac{k\left(-\mathrm{4}a+\varpi +b+r\right)}{\mathrm{2}a\varpi \left(k-br\right)},\\ {\beta }_{\mathrm{02}}=\frac{k\left(a\varpi +a\left(-\left(b+r\right)\right)-br+k\right)}{a\varpi {\left(k-br\right)}^{\mathrm{2}}}.\end{array}$$ The dynamics of the system on the center manifold is then described by $$\begin{array}{}\text{(27)}& \left(\begin{array}{c}{\dot{u}}_{\mathrm{1}}\\ {\dot{u}}_{\mathrm{2}}\end{array}\right)=\left(\begin{array}{c}{u}_{\mathrm{2}}\\ \mathrm{0}\end{array}\right)+\frac{\mathrm{1}}{\Psi }\left(\begin{array}{c}G\left(u_{\mathrm{2}},u_{\mathrm{3}}\right)\\ aG\left(u_{\mathrm{2}},u_{\mathrm{3}}\right)\end{array}\right),\end{array}$$where$$\begin{array}{c}\text{(28)}& \Psi =\mathrm{2}\left(\frac{b^{\mathrm{2}}}{\mathrm{4}\varpi }-\frac{br}{\mathrm{2}\varpi }+\frac{\mathrm{1}}{\mathrm{4}}\varpi +\frac{k}{\varpi }+\frac{r^{\mathrm{2}}}{\mathrm{4}\varpi }+\frac{b}{\mathrm{2}}+\frac{r}{\mathrm{2}}\right),G\left(u_{\mathrm{2}},u_{\mathrm{3}}\right)=\frac{\mathrm{2}\varpi u_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}{r}^{\mathrm{5}}}{a\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\varpi u_{\mathrm{1}}^{\mathrm{3}}{r}^{\mathrm{5}}}{\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\varpi u_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{5}}}{a^{\mathrm{2}}\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}+\frac{\mathrm{6}\varpi u_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{4}}}{a\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}+\frac{\mathrm{6}{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{4}}}{a{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }-\frac{\mathrm{4}\varpi u_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}{r}^{\mathrm{4}}}{\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\mathrm{2}\varpi u_{\mathrm{2}}^{\mathrm{3}}{r}^{\mathrm{4}}}{a^{\mathrm{2}}\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\mathrm{4}{u}_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}{r}^{\mathrm{4}}}{{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }-\frac{\mathrm{2}{u}_{\mathrm{2}}^{\mathrm{3}}{r}^{\mathrm{4}}}{a^{\mathrm{2}}{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{2}b{u}_{\mathrm{2}}^{\mathrm{3}}{r}^{\mathrm{3}}}{a^{\mathrm{2}}{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{4}\varpi u_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{3}}}{\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}+\frac{\mathrm{4}{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{3}}}{{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{4}k\varpi u_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}{r}^{\mathrm{3}}}{a\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}+\frac{\mathrm{4}b{u}_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}{r}^{\mathrm{3}}}{{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }-\frac{\mathrm{2}k\varpi u_{\mathrm{1}}^{\mathrm{3}}{r}^{\mathrm{3}}}{\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\mathrm{4}\varpi u_{\mathrm{2}}^{\mathrm{3}}{r}^{\mathrm{3}}}{a\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\mathrm{2}k\varpi u_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{3}}}{a^{\mathrm{2}}\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\mathrm{4}{u}_{\mathrm{2}}^{\mathrm{3}}{r}^{\mathrm{3}}}{a{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }-\frac{\mathrm{6}b{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{3}}}{a{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{4}b{u}_{\mathrm{2}}^{\mathrm{3}}{r}^{\mathrm{2}}}{a{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{2}k{u}_{\mathrm{2}}^{\mathrm{3}}{r}^{\mathrm{2}}}{a^{\mathrm{2}}{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{6}k\varpi u_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{2}}}{a\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}+\frac{\mathrm{4}k{u}_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}{r}^{\mathrm{2}}}{{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }-\frac{\mathrm{4}k\varpi u_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}{r}^{\mathrm{2}}}{\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\mathrm{2}k\varpi u_{\mathrm{2}}^{\mathrm{3}}{r}^{\mathrm{2}}}{a^{\mathrm{2}}\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\mathrm{4}b{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{2}}}{{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }-\frac{\mathrm{6}k{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}{r}^{\mathrm{2}}}{a{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{4}k\varpi u_{\mathrm{2}}^{\mathrm{3}}r}{a\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}+\frac{\mathrm{20}k{u}_{\mathrm{2}}^{\mathrm{3}}r}{a{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{2}bk{u}_{\mathrm{2}}^{\mathrm{3}}r}{a^{\mathrm{2}}{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{2}\varpi u_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}r}{{\left(br-k\right)}^{\mathrm{2}}}+\frac{\mathrm{2}{k}^{\mathrm{2}}\varpi u_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}r}{a\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}+\frac{\mathrm{4}bk{u}_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}r}{{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }-\frac{\mathrm{2}\varpi u_{\mathrm{2}}^{\mathrm{3}}r}{a{\left(br-k\right)}^{\mathrm{2}}}-\frac{k^{\mathrm{2}}\varpi u_{\mathrm{1}}^{\mathrm{3}}r}{\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\mathrm{4}k\varpi u_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}r}{\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{k^{\mathrm{2}}\varpi u_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}r}{a^{\mathrm{2}}\left(br-k\right){\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}}-\frac{\mathrm{20}k{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}r}{{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\sqrt{b^{\mathrm{2}}-\mathrm{2}rb+r^{\mathrm{2}}+\mathrm{4}k}}-\frac{\mathrm{6}bk{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}r}{a{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{4}{k}^{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{3}}}{a^{\mathrm{2}}{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{4}bk{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}}{{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{\mathrm{8}{k}^{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}}{{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }+\frac{-a^{\mathrm{2}}{r}^{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{3}}-a^{\mathrm{2}}k{u}_{\mathrm{1}}^{\mathrm{3}}+\mathrm{2}a{r}^{\mathrm{2}}{u}_{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}ak{u}_{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{2}}+\mathrm{2}{a}^{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}}{u}_{\mathrm{1}}-r^{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{2}}{u}_{\mathrm{1}}-k{u}_{\mathrm{2}}^{\mathrm{2}}{u}_{\mathrm{1}}-\mathrm{2}a{u}_{\mathrm{2}}^{\mathrm{3}}}{a^{\mathrm{2}}\left(br-k\right)}-\frac{u_{\mathrm{1}}{\left(a{u}_{\mathrm{1}}-u_{\mathrm{2}}\right)}^{\mathrm{2}}}{a^{\mathrm{2}}}-\frac{\mathrm{2}\left(-r^{\mathrm{2}}{u}_{\mathrm{2}}^{\mathrm{3}}-k{u}_{\mathrm{2}}^{\mathrm{3}}+a{r}^{\mathrm{2}}{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}+ak{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}\right)}{a{\left(br-k\right)}^{\mathrm{2}}}-\frac{\mathrm{4}bk{u}_{\mathrm{2}}^{\mathrm{3}}}{a{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }-\frac{\mathrm{12}{k}^{\mathrm{2}}{u}_{\mathrm{1}}{u}_{\mathrm{2}}^{\mathrm{2}}}{a{\left(r^{\mathrm{2}}+k\right)}^{\mathrm{2}}\varpi }.\end{array}$$

Conditional normal form at bifurcation values of parameters can be obtained by following Kuznetsov approach which implies that the simplified system which represents the conditional normal form at bifurcation value is the following system: $$\begin{array}{c}\text{(29)}& {\dot{v}}_{\mathrm{1}}=v_{\mathrm{2}},{\dot{v}}_{\mathrm{2}}=\frac{ar}{\mathrm{1}-br}{v}_{\mathrm{1}}^{\mathrm{3}}+\frac{r\left(a\left(-\mathrm{3}{b}^{\mathrm{2}}+\mathrm{6}br-\mathrm{3}{r}^{\mathrm{2}}-\mathrm{12}\right)-\mathrm{3}{b}^{\mathrm{2}}+\mathrm{6}br-\mathrm{3}{r}^{\mathrm{2}}-\mathrm{12}\right)}{\mathrm{1}{b}^{\mathrm{3}}r+b^{\mathrm{2}}\left(-\mathrm{2}{r}^{\mathrm{2}}-\mathrm{1}\right)+br\left(\mathrm{1}{r}^{\mathrm{2}}+\mathrm{6}\right)-\mathrm{1}{r}^{\mathrm{2}}-\mathrm{4}}{v}_{\mathrm{1}}^{\mathrm{2}}{v}_{\mathrm{2}}\end{array}$$

The next step in this analysis is to obtain the universal unfolding of the BT bifurcation. In order to put the system in an appropriate form, the bifurcation parameters should be perturbed around bifurcation value, such that $c=-a+{\Delta }_{\mathrm{1}}$ and $d=-a+{\Delta }_{\mathrm{2}}.$ The following system is attained$$\begin{array}{c}\text{(30)}& {\dot{x}}_{\mathrm{1}}=a\left(-x_{\mathrm{1}}+x_{\mathrm{2}}\right),{\dot{x}}_{\mathrm{2}}=\left(-a+{\Delta }_{\mathrm{1}}\right)x_{\mathrm{2}}-\left(-a+{\Delta }_{\mathrm{2}}\right)x_{\mathrm{2}}-x_{\mathrm{1}}{x}_{\mathrm{3}},\\ {\dot{x}}_{\mathrm{3}}=x_{\mathrm{4}}+x_{\mathrm{1}}{x}_{\mathrm{2}}-b{x}_{\mathrm{3}},\\ {\dot{x}}_{\mathrm{4}}=-r{x}_{\mathrm{4}}+k{x}_{\mathrm{3}},\\ {\dot{\Delta }}_{\mathrm{1}}=\mathrm{0},\\ {\dot{\Delta }}_{\mathrm{2}}=\mathrm{0}.\end{array}$$ Now, applying transformation (22) to get$$\begin{array}{c}\text{(31)}& \left(\begin{array}{c}{\dot{u}}_{\mathrm{1}}\\ {\dot{u}}_{\mathrm{2}}\\ {\dot{u}}_{\mathrm{3}}\\ {\dot{u}}_{\mathrm{4}}\end{array}\right)=\left(\begin{array}{c}{\widetilde{ϝ}}_{\mathrm{1}}\\ {\widetilde{ϝ}}_{\mathrm{2}}\\ {\widetilde{ϝ}}_{\mathrm{3}}\\ {\widetilde{ϝ}}_{\mathrm{4}}\end{array}\right),\text{(32)}& {\dot{\Delta }}_{\mathrm{1}}=\mathrm{0},{\dot{\Delta }}_{\mathrm{2}}=\mathrm{0},\end{array}$$where $$\begin{array}{c}\text{(33)}& {\widetilde{ϝ}}_{\mathrm{1}}=\frac{u_{\mathrm{4}}{u}_{\mathrm{2}}\varpi }{\mathrm{2}ak}-\frac{u_{\mathrm{3}}{u}_{\mathrm{2}}\varpi }{\mathrm{2}ak}-\frac{b{u}_{\mathrm{3}}{u}_{\mathrm{2}}}{\mathrm{2}ak}-\frac{b{u}_{\mathrm{4}}{u}_{\mathrm{2}}}{\mathrm{2}ak}+\frac{r{u}_{\mathrm{3}}{u}_{\mathrm{2}}}{\mathrm{2}ak}+\frac{r{u}_{\mathrm{4}}{u}_{\mathrm{2}}}{\mathrm{2}ak}-\frac{{\Delta }_{\mathrm{1}}{u}_{\mathrm{2}}}{a}+\frac{u_{\mathrm{1}}{u}_{\mathrm{3}}\varpi }{\mathrm{2}k}-\frac{u_{\mathrm{1}}{u}_{\mathrm{4}}\varpi }{\mathrm{2}k}+\frac{b{u}_{\mathrm{1}}{u}_{\mathrm{3}}}{\mathrm{2}k}+\frac{b{u}_{\mathrm{1}}{u}_{\mathrm{4}}}{\mathrm{2}k}-\frac{r{u}_{\mathrm{1}}{u}_{\mathrm{3}}}{\mathrm{2}k}-\frac{r{u}_{\mathrm{1}}{u}_{\mathrm{4}}}{\mathrm{2}k}+{\Delta }_{\mathrm{1}}{u}_{\mathrm{1}}-{\Delta }_{\mathrm{2}}{u}_{\mathrm{1}}+u_{\mathrm{2}},{\widetilde{ϝ}}_{\mathrm{2}}=\frac{a{u}_{\mathrm{1}}{u}_{\mathrm{3}}\varpi }{\mathrm{2}k}-\frac{a{u}_{\mathrm{1}}{u}_{\mathrm{4}}\varpi }{\mathrm{2}k}+\frac{ab{u}_{\mathrm{1}}{u}_{\mathrm{3}}}{\mathrm{2}k}+\frac{ab{u}_{\mathrm{1}}{u}_{\mathrm{4}}}{\mathrm{2}k}-\frac{ar{u}_{\mathrm{1}}{u}_{\mathrm{3}}}{\mathrm{2}k}-\frac{ar{u}_{\mathrm{1}}{u}_{\mathrm{4}}}{\mathrm{2}k}+a{\Delta }_{\mathrm{1}}{u}_{\mathrm{1}}-a{\Delta }_{\mathrm{2}}{u}_{\mathrm{1}}-\frac{u_{\mathrm{2}}{u}_{\mathrm{3}}\varpi }{\mathrm{2}k}+\frac{u_{\mathrm{2}}{u}_{\mathrm{4}}\varpi }{\mathrm{2}k}-\frac{b{u}_{\mathrm{2}}{u}_{\mathrm{3}}}{\mathrm{2}k}-\frac{b{u}_{\mathrm{2}}{u}_{\mathrm{4}}}{\mathrm{2}k}+\frac{r{u}_{\mathrm{2}}{u}_{\mathrm{3}}}{\mathrm{2}k}+\frac{r{u}_{\mathrm{2}}{u}_{\mathrm{4}}}{\mathrm{2}k}-{\Delta }_{\mathrm{1}}{u}_{\mathrm{2}},\\ {\widetilde{ϝ}}_{\mathrm{3}}=\frac{k{u}_{\mathrm{1}}{u}_{\mathrm{2}}}{a\varpi }-\frac{b^{\mathrm{2}}{u}_{\mathrm{3}}}{\mathrm{4}\varpi }-\frac{b^{\mathrm{2}}{u}_{\mathrm{4}}}{\mathrm{4}\varpi }+\frac{br{u}_{\mathrm{3}}}{\mathrm{2}\varpi }+\frac{br{u}_{\mathrm{4}}}{\mathrm{2}\varpi }-\frac{\mathrm{1}}{\mathrm{4}}{u}_{\mathrm{3}}\varpi +\frac{\mathrm{1}}{\mathrm{4}}{u}_{\mathrm{4}}\varpi -\frac{k{u}_{\mathrm{1}}^{\mathrm{2}}}{\varpi }-\frac{k{u}_{\mathrm{3}}}{\varpi }-\frac{k{u}_{\mathrm{4}}}{\varpi }-\frac{r^{\mathrm{2}}{u}_{\mathrm{3}}}{\mathrm{4}\varpi }-\frac{r^{\mathrm{2}}{u}_{\mathrm{4}}}{\mathrm{4}\varpi }-\frac{b{u}_{\mathrm{3}}}{\mathrm{2}}-\frac{r{u}_{\mathrm{3}}}{\mathrm{2}},\\ {\widetilde{ϝ}}_{\mathrm{4}}=-\frac{k{u}_{\mathrm{1}}{u}_{\mathrm{2}}}{a\varpi }+\frac{b^{\mathrm{2}}{u}_{\mathrm{3}}}{\mathrm{4}\varpi }+\frac{b^{\mathrm{2}}{u}_{\mathrm{4}}}{\mathrm{4}\varpi }-\frac{br{u}_{\mathrm{3}}}{\mathrm{2}\varpi }-\frac{br{u}_{\mathrm{4}}}{\mathrm{2}\varpi }+\frac{k{u}_{\mathrm{1}}^{\mathrm{2}}}{\varpi }-\frac{\mathrm{1}}{\mathrm{4}}{u}_{\mathrm{3}}\varpi +\frac{k{u}_{\mathrm{3}}}{\varpi }+\frac{r^{\mathrm{2}}{u}_{\mathrm{3}}}{\mathrm{4}\varpi }+\frac{\mathrm{1}}{\mathrm{4}}{u}_{\mathrm{4}}\varpi +\frac{k{u}_{\mathrm{4}}}{\varpi }+\frac{r^{\mathrm{2}}{u}_{\mathrm{4}}}{\mathrm{4}\varpi }-\frac{b{u}_{\mathrm{4}}}{\mathrm{2}}-\frac{r{u}_{\mathrm{4}}}{\mathrm{2}}.\end{array}$$

The center manifold is assumed as follows: $$\begin{array}{c}\text{(34)}& u_{\mathrm{3}}={\displaystyle \sum }_{i,j,k,l=\mathrm{0}}^{i+j+k+l=\mathrm{2}}{\tilde{\alpha }}_{ijkl}{u}_{\mathrm{1}}^i{u}_{\mathrm{2}}^j{\Delta }_{\mathrm{1}}^k{\Delta }_{\mathrm{2}}^l,u_{\mathrm{4}}={\displaystyle \sum }_{i,j,k,l=\mathrm{0}}^{i+j+k+l=\mathrm{2}}{\tilde{\beta }}_{ijkl}{u}_{\mathrm{1}}^i{u}_{\mathrm{2}}^j{\Delta }_{\mathrm{1}}^k{\Delta }_{\mathrm{2}}^l,\end{array}$$where parameters ${\tilde{\alpha }}_{ijkl}$ and ${\tilde{\beta }}_{ijkl}$ can be attained through the same way employed in the first part of this subsection. The following reduced system on center manifold is obtained after some simplifications $$\begin{array}{c}\text{(35)}& \left(\begin{array}{c}{\dot{u}}_{\mathrm{1}}\\ {\dot{u}}_{\mathrm{2}}\end{array}\right)=A\left(\begin{array}{c}{u}_{\mathrm{1}}\\ u_{\mathrm{2}}\end{array}\right)+\left(\begin{array}{c}{\tilde{G}}_{\mathrm{1}}\left(u_{\mathrm{1}},u_{\mathrm{2}},{\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)\\ {\tilde{G}}_{\mathrm{2}}\left(u_{\mathrm{1}},u_{\mathrm{2}},{\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)\end{array}\right),A=\left(\begin{array}{cc}\mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}\end{array}\right),\end{array}$$where (36) G ~ 1 u 1 , u 2 , Δ 1 , Δ 2 = u 1 2 u 2 b 3 a b 2 + ϖ ∗ - 2 r b + r r + ϖ ∗ + 16 2 + ϖ ∗ u 1 2 u 2 b 2 a b 2 + ϖ ∗ - 2 r b + r r + ϖ ∗ + 16 2 - 3 r u 1 2 u 2 b 2 a b 2 + ϖ ∗ - 2 r b + r r + ϖ ∗ + 16 2 + u 1 3 b 2 b r - 4 + b 2 - 2 r b + r 2 + 16 u 1 2 u 2 b 2 a ϖ ∗ b r - 4 + 4 ϖ ∗ u 1 2 u 2 b b 2 + ϖ ∗ - 2 r b + r r + ϖ ∗ + 16 2 + 3 r 2 u 1 2 u 2 b a b 2 + ϖ ∗ - 2 r b + r r + ϖ ∗ + 16 2 + 16 u 1 2 u 2 b a b 2 + ϖ ∗ - 2 r b + r r + ϖ ∗ + 16 2 - u 1 2 u 2 b 2 a b r - 4 - b 2 - 2 r b + r 2 + 16 u 1 3 b 2 ϖ ∗ b r - 4 + r b 2 - 2 r b + r 2 + 16 u 1 2 u 2 2 a ϖ ∗ b r - 4 + r u 1 2 u 2 2 a b r - 4 - r u 1 3 2 b r - 4 - r b 2 - 2 r b + r 2 + 16 u 1 3 2 ϖ ∗ b r - 4 - 4 ϖ ∗ r u 1 2 u 2 b 2 + ϖ ∗ - 2 r b + r r + ϖ ∗ + 16 2 - 16 r u 1 2 u 2 a b 2 + ϖ ∗ - 2 r b + r r + ϖ ∗ + 16 2 , G ~ 2 u 1 , u 2 , Δ 1 , Δ 2 = - a b u 1 3 b 2 - 2 b r + r 2 + 16 2 ϖ ∗ b r - 4 - a r u 1 3 b 2 - 2 b r + r 2 + 16 2 ϖ ∗ b r - 4 + 4 a b u 1 2 u 2 ϖ ∗ b 2 + b ϖ ∗ - 2 r + r ϖ ∗ + r + 16 2 + a b u 1 3 2 b r - 4 - a r u 1 3 2 b r - 4 + b u 1 2 u 2 b 2 - 2 b r + r 2 + 16 2 ϖ ∗ b r - 4 + r u 1 2 u 2 b 2 - 2 b r + r 2 + 16 2 ϖ ∗ b r - 4 + b 2 u 1 2 u 2 ϖ ∗ b 2 + b ϖ ∗ - 2 r + r ϖ ∗ + r + 16 2 + 16 b u 1 2 u 2 b 2 + b ϖ ∗ - 2 r + r ϖ ∗ + r + 16 2 - 16 r u 1 2 u 2 b 2 + b ϖ ∗ - 2 r + r ϖ ∗ + r + 16 2 + b 3 u 1 2 u 2 b 2 + b ϖ ∗ - 2 r + r ϖ ∗ + r + 16 2 - b u 1 2 u 2 2 b r - 4 + r u 1 2 u 2 2 b r - 4 ,such that ${\varpi }^{\ast }=\sqrt{(b-r{)}^{\mathrm{2}}+\mathrm{16}}$.

Now, it necessary to define the next vector valued functions$$\begin{array}{c}\text{(37)}& C\left(s,v,w\right)={\displaystyle \sum }_{i,j,k=\mathrm{1}}^{\mathrm{3}}{\left.\frac{{\partial }^{\mathrm{3}}\tilde{G}\left(u_{\mathrm{1}},u_{\mathrm{2}},{\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)}{\partial u_i\partial u_j\partial u_k}\right|}_{\begin{array}{c}{u}_{\mathrm{1}}=u_{\mathrm{2}}=\mathrm{0}\\ {\Delta }_{\mathrm{1}}={\Delta }_{\mathrm{2}}=\mathrm{0}\end{array}}{s}_i{v}_j{w}_k,A_{\mathrm{1}}\left(s,\Lambda \right)={\displaystyle \sum }_{i=\mathrm{1}}^{\mathrm{2}}{\displaystyle \sum }_{j=\mathrm{1}}^{\mathrm{2}}{\left.\frac{{\partial }^{\mathrm{2}}\tilde{G}\left(u_{\mathrm{1}},u_{\mathrm{2}},{\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)}{\partial u_i\partial {\Delta }_j}\right|}_{\begin{array}{c}{u}_{\mathrm{2}}=u_{\mathrm{3}}=\mathrm{0}\\ {\Delta }_{\mathrm{1}}={\Delta }_{\mathrm{2}}=\mathrm{0}\end{array}}{s}_i{\Lambda }_j,\\ A_{\mathrm{2}}\left(s,\Lambda ,\Omega \right)={\displaystyle \sum }_{i=\mathrm{1}}^{\mathrm{2}}\hspace{0.17em}{\displaystyle \sum }_{j,k=\mathrm{1}}^{\mathrm{2}}{\left.\frac{{\partial }^{\mathrm{2}}\tilde{G}\left(u_{\mathrm{2}},u_{\mathrm{3}},{\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)}{\partial u_i\partial {\Delta }_j\partial {\Delta }_k}\right|}_{\begin{array}{c}{u}_{\mathrm{2}}=u_{\mathrm{3}}=\mathrm{0}\\ {\Delta }_{\mathrm{1}}={\Delta }_{\mathrm{2}}=\mathrm{0}\end{array}}{s}_i{\Lambda }_j{\Omega }_k,\end{array}$$where $$\begin{array}{}\text{(38)}& \tilde{G}\left(u_{\mathrm{2}},u_{\mathrm{3}},{\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)=\left(\begin{array}{c}{\tilde{G}}_{\mathrm{1}}\left(u_{\mathrm{1}},u_{\mathrm{2}},{\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)\\ {\tilde{G}}_{\mathrm{2}}\left(u_{\mathrm{1}},u_{\mathrm{2}},{\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)\end{array}\right).\end{array}$$The exact expression for these vector valued functions is omitted for brevity. Now, following the procedure described in [39], we assume that the following relation holds between unfolding parameters $\eta =({\eta }_{\mathrm{1}}$, ${\eta }_{\mathrm{2}}{)}^T$ and perturbation values $\Delta =({\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}{)}^T$$$\begin{array}{}\text{(39)}& {\Delta }_i={\rho }_i{\eta }_{\mathrm{1}}+{\upsilon }_i{\eta }_{\mathrm{2}}+{\gamma }_i{\eta }_{\mathrm{1}}^{\mathrm{2}}+{\delta }_i{\eta }_{\mathrm{1}}{\eta }_{\mathrm{2}}+{\zeta }_i{\eta }_{\mathrm{2}}^{\mathrm{2}}+O\left({‖\eta ‖}^{\mathrm{3}}\right),\hspace{1em}i=\mathrm{1,2}\end{array}$$

or equivalently$$\begin{array}{}\text{(40)}& {\eta }_i={\tilde{\rho }}_i{\Delta }_{\mathrm{1}}+{\tilde{\upsilon }}_i{\Delta }_{\mathrm{2}}+{\tilde{\gamma }}_i{\Delta }_{\mathrm{1}}^{\mathrm{2}}+{\tilde{\delta }}_i{\Delta }_{\mathrm{1}}{\Delta }_{\mathrm{2}}+{\tilde{\zeta }}_i{\Delta }_{\mathrm{2}}^{\mathrm{2}}+O\left({‖\Delta ‖}^{\mathrm{3}}\right),\hspace{1em}i=\mathrm{1,2}.\end{array}$$

The two real linearly independent generalized eigenvectors $q_{\mathrm{1}}$ and $q_{\mathrm{2}}$ that satisfy $A{q}_{\mathrm{1}}=\mathrm{0}$ and $A{q}_{\mathrm{2}}=q_{\mathrm{1}}$ can be found as $$\begin{array}{c}\text{(41)}& q_{\mathrm{1}}=\left[\begin{array}{c}\mathrm{1}\\ \mathrm{0}\end{array}\right],q_{\mathrm{2}}=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{1}\end{array}\right].\end{array}$$Also, for the transposed matrix $A^T$, there exist eigenvectors $p_{\mathrm{1}}$ and $p_{\mathrm{2}}$ that have the properties $$\begin{array}{c}\text{(42)}& A^T{p}_{\mathrm{2}}=\mathrm{0},A^T{p}_{\mathrm{1}}=p_{\mathrm{2}},\\ p_{\mathrm{1}}^T{q}_{\mathrm{1}}=p_{\mathrm{2}}^T{q}_{\mathrm{2}}=\mathrm{1},\\ p_{\mathrm{1}}^T{q}_{\mathrm{2}}=p_{\mathrm{2}}^T{q}_{\mathrm{1}}=\mathrm{0},\end{array}$$ and they are computed as $$\begin{array}{c}\text{(43)}& p_{\mathrm{1}}=\left[\begin{array}{c}\mathrm{1}\\ \mathrm{0}\end{array}\right],p_{\mathrm{2}}=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{1}\end{array}\right].\end{array}$$

After some calculations, the linear part coefficients of (39) can be evaluated as $$\begin{array}{c}\text{(44)}& {\rho }_{\mathrm{1}}=\frac{\mathrm{1}}{a},{\rho }_{\mathrm{2}}=\mathrm{0},\\ {\upsilon }_{\mathrm{1}}=-\mathrm{1},\\ {\upsilon }_{\mathrm{2}}=-\mathrm{1}.\end{array}$$ Then, the transversality condition is verified by evaluating $$\begin{array}{}\text{(45)}& \det \left(\begin{array}{cc}\frac{\mathrm{1}}{a}& -\mathrm{1}\\ \mathrm{0}& -\mathrm{1}\end{array}\right)=-\frac{\mathrm{1}}{a}\ne \mathrm{0},\end{array}$$which implies that transversality condition is fulfilled. The other coefficients of (39) can be evaluated by solving a set of systems of linear equations described in [39] which results in $$\begin{array}{c}\text{(46)}& {\gamma }_{\mathrm{1}}={\gamma }_{\mathrm{2}}=\mathrm{0},{\delta }_{\mathrm{1}}={\delta }_{\mathrm{2}}=\mathrm{0},\\ {\zeta }_{\mathrm{1}}={\zeta }_{\mathrm{2}}=\mathrm{0}.\end{array}$$By substituting from (39) into (40) and comparing the coefficients of ${\eta }_{\mathrm{1}}^i{\eta }_{\mathrm{2}}^j$, $i+j=\mathrm{1,2}$, on both sides of (40), the unfolding parameters ${\eta }_{\mathrm{1}}({\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}})$ and ${\eta }_{\mathrm{2}}({\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}})$ can be obtained in the following expression $$\begin{array}{c}\text{(47)}& {\eta }_{\mathrm{1}}\left({\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)=a{\Delta }_{\mathrm{1}}-a{\Delta }_{\mathrm{2}},{\eta }_{\mathrm{2}}\left({\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\right)=-a{\Delta }_{\mathrm{2}}.\end{array}$$

The system representing the universal unfolding of BT bifurcation is thus written as$$\begin{array}{c}\text{(48)}& {\dot{u}}_{\mathrm{1}}=u_{\mathrm{2}},{\dot{u}}_{\mathrm{2}}=\left(a{\Delta }_{\mathrm{1}}-a{\Delta }_{\mathrm{2}}\right)u_{\mathrm{1}}-a{\Delta }_{\mathrm{2}}{u}_{\mathrm{2}}+\frac{ar}{\mathrm{1}-br}{u}_{\mathrm{1}}^{\mathrm{3}}+\frac{r\left(a\left(-\mathrm{3}{b}^{\mathrm{2}}+\mathrm{6}br-\mathrm{3}{r}^{\mathrm{2}}-\mathrm{12}\right)-\mathrm{3}{b}^{\mathrm{2}}+\mathrm{6}br-\mathrm{3}{r}^{\mathrm{2}}-\mathrm{12}\right)}{\mathrm{1}{b}^{\mathrm{3}}r+b^{\mathrm{2}}\left(-\mathrm{2}{r}^{\mathrm{2}}-\mathrm{1}\right)+br\left(\mathrm{1}{r}^{\mathrm{2}}+\mathrm{6}\right)-\mathrm{1}{r}^{\mathrm{2}}-\mathrm{4}}{u}_{\mathrm{1}}^{\mathrm{2}}{u}_{\mathrm{2}}.\end{array}$$Let ${\vartheta }_{\mathrm{30}}=ar/\left(\mathrm{1}-br\right)$ and ${\vartheta }_{\mathrm{21}}=r\left(a\left(-\mathrm{3}{b}^{\mathrm{2}}+\mathrm{6}br-\mathrm{3}{r}^{\mathrm{2}}-\mathrm{12}\right)-\mathrm{3}{b}^{\mathrm{2}}+\mathrm{6}br-\mathrm{3}{r}^{\mathrm{2}}-\mathrm{12}\right)/\left(\mathrm{1}{b}^{\mathrm{3}}r+b^{\mathrm{2}}\left(-\mathrm{2}{r}^{\mathrm{2}}-\mathrm{1}\right)+br\left(\mathrm{1}{r}^{\mathrm{2}}+\mathrm{6}\right)-\mathrm{1}{r}^{\mathrm{2}}-\mathrm{4}\right)$, then apply the following rescaling of coordinates and time to further simplify system (48) $$\begin{array}{c}\text{(49)}& u_{\mathrm{1}}=\frac{s_{\mathrm{1}}\sqrt{\left|{\vartheta }_{\mathrm{30}}\right|}}{\left|{\vartheta }_{\mathrm{21}}\right|}{\xi }_{\mathrm{1}},u_{\mathrm{2}}=-\frac{s_{\mathrm{1}}{s}_{\mathrm{2}}{\left|{\vartheta }_{\mathrm{30}}\right|}^{\mathrm{3}/\mathrm{2}}}{{\left|{\vartheta }_{\mathrm{21}}\right|}^{\mathrm{2}}}{\xi }_{\mathrm{2}},\\ t=-\frac{{\vartheta }_{\mathrm{21}}}{\left|{\vartheta }_{\mathrm{30}}\right|}\tau .\end{array}$$ After simplification, we can obtain$$\begin{array}{c}\text{(50)}& {\dot{\xi }}_{\mathrm{1}}={\xi }_{\mathrm{2}},{\dot{\xi }}_{\mathrm{2}}={\mu }_{\mathrm{1}}{\xi }_{\mathrm{1}}+{\mu }_{\mathrm{2}}{\xi }_{\mathrm{2}}+s_{\mathrm{1}}{\xi }_{\mathrm{1}}^{\mathrm{3}}-{\xi }_{\mathrm{1}}^{\mathrm{2}}{\xi }_{\mathrm{2}},\end{array}$$where derivative is considered with respect to $\tau .$

The bifurcation structure of (50) is summarized in the next theorem [39–41].

Theorem 3.

The bifurcation diagram of (50) in ${\mu }_{\mathrm{1}}-{\mu }_{\mathrm{2}}$ plane involves the following bifurcation curves:

(i) $PB=\{({\mu }_{\mathrm{1}},{\mu }_{\mathrm{2}}):{\mu }_{\mathrm{1}}=\mathrm{0},{\mu }_{\mathrm{2}}\ne \mathrm{0}\}$,

(ii) $H{B}_{\mathrm{1}}=\{({\mu }_{\mathrm{1}},{\mu }_{\mathrm{2}}):{\mu }_{\mathrm{2}}=\mathrm{0},{\mu }_{\mathrm{1}}<\mathrm{0}\}$,

(iii) $H{B}_{\mathrm{2}}=\{({\mu }_{\mathrm{1}},{\mu }_{\mathrm{2}}):{\mu }_{\mathrm{2}}={\mu }_{\mathrm{1}}+O({\mu }_{\mathrm{1}}^{\mathrm{2}}),{\mu }_{\mathrm{1}}>\mathrm{0}\}$,

(iv) $H_{L}B=\{({\mu }_{\mathrm{1}},{\mu }_{\mathrm{2}}):{\mu }_{\mathrm{2}}=\left(\mathrm{4}/\mathrm{5}\right){\mu }_{\mathrm{1}}+O({\mu }_{\mathrm{1}}^{\mathrm{3}/\mathrm{2}}),{\mu }_{\mathrm{1}}>\mathrm{0}\}$,

(v) $SNB=\{({\mu }_{\mathrm{1}},{\mu }_{\mathrm{2}}):{\eta }_{\mathrm{2}}=\mathrm{0.752}{\mu }_{\mathrm{1}}+O({\mu }_{\mathrm{1}}^{\mathrm{3}/\mathrm{2}}),{\mu }_{\mathrm{1}}>\mathrm{0}\}$, where ${\mu }_{\mathrm{1}}=\mathrm{9}(a{\Delta }_{\mathrm{1}}+a^{\mathrm{2}}({\Delta }_{\mathrm{1}}^{\mathrm{2}}+{\Delta }_{\mathrm{1}}{\Delta }_{\mathrm{2}}))$ and ${\mu }_{\mathrm{2}}=\mathrm{3}(a({\Delta }_{\mathrm{1}}+{\Delta }_{\mathrm{2}})+a^{\mathrm{2}}(\mathrm{2}{\Delta }_{\mathrm{1}}^{\mathrm{2}}+\mathrm{3}{\Delta }_{\mathrm{1}}{\Delta }_{\mathrm{2}}+{\Delta }_{\mathrm{2}}^{\mathrm{2}}))$ where $PB$ represents a pitchfork bifurcation of system (50), $H{B}_{\mathrm{1}}$ denotes Hopf bifurcation of the origin equilibrium point, $H{B}_{\mathrm{2}}$ refers to the Hopf bifurcation of the nontrivial equilibrium points, $H_{L}B$ corresponds to Homoclinic connection to the origin, and finally $SNB$ denotes the saddle-node bifurcation of symmetric periodic orbits.

The phase portraits of system (1) are obtained via numerical simulations and illustrate examples of different dynamics that can be exhibited by the hyperchaotic system. In particular, Figure 1 verifies the occurrences of Hopf and pitchfork bifurcations of system (2) whereas Figure 2 shows 2D and 3D portraits of hyperchaotic attractor of the system.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

Bifurcation diagrams are used to provide useful illustrations for variations of a particular system dynamics with changes in its parameters. Figure 3 shows examples of bifurcation diagrams of system (1) with respect to different parameters in the system. Figure 4 depicts the associated Lyapunov exponent plots to the cases presented in Figure 3.

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

4. First Integrals of the System

It is well known that, if the analyzed system has chaotic dynamics, the behaviors of the system are irregular behaviors. So, it is impossible to predict the behavior of the system in a long time scale. Therefore, it is important to ask whether there is any set of parameter values in which system dynamics are regular and that the studied system is integrable. In this section, we provide the question of global analytic integrability investigated for hyperchaotic system (1). In particular, we provide conditions of analytic integrability in hyperchaotic system (1) and obtain the possible first integrals of it.

Theorem 4.

System (1) has the following first integrals:

(1)

$I_{\mathrm{1}}(x,y,z,w,t)=w^{\mathrm{2}}-k(y^{\mathrm{2}}+z^{\mathrm{2}})$ when $c=b=d=r=\mathrm{0}$,

(2)

$I_{\mathrm{2}}(x,y,z,w,t)=e^{\mathrm{2}bt}\left(b{w}^{\mathrm{2}}-kb\left(y^{\mathrm{2}}+z^{\mathrm{2}}\right)+ck{x}^{\mathrm{2}}\right)$ when $a=d=r=b$,

(3)

$I_{\mathrm{3}}(x,y,z,w,t)=e^{\mathrm{2}bt}\left(\mathrm{2}{c}_{\mathrm{3}}\left(w^{\mathrm{2}}+cb{x}^{\mathrm{2}}-b^{\mathrm{2}}\left(y^{\mathrm{2}}+z^{\mathrm{2}}\right)\right)+c_{\mathrm{2}}\left(-\mathrm{2}bz+\mathrm{2}w+x^{\mathrm{2}}\right)\right)$ when $a=c=d=r=b$ and $k=b^{\mathrm{2}}$,

(4)

$I_{\mathrm{4}}(x,y,z,w,t)=e^{\mathrm{2}at}\left(x^{\mathrm{2}}(r-\mathrm{2}a)+\mathrm{2}az(\mathrm{2}a-r)-\mathrm{2}aw\right)$ when $k=(\mathrm{2}a-b)(\mathrm{2}a-r)$,

(5)

$I_{\mathrm{5}}(x,y,z,w,t)=e^{\mathrm{2}bt}\left(\mathrm{2}cb(bz-w)-b^{\mathrm{2}}\left(y^{\mathrm{2}}+z^{\mathrm{2}}\right)+w^{\mathrm{2}}\right)$ when $c=d=r=b$ and $k=b^{\mathrm{2}}$,

where $c_{\mathrm{2}}$ and $c_{\mathrm{3}}$ are arbitrary constants.

Proof.

The function $I(x(t),y(t),z(t),w(t),t)$ is said to be the first integral of system (1) if it satisfies the condition$$\begin{array}{}\text{(51)}& \frac{dI}{dx}=I_t+\dot{x}{I}_x+\dot{y}{I}_y+\dot{z}{I}_z+\dot{w}{I}_w=\mathrm{0},\end{array}$$Substituting (1) into (51), we obtain$$\begin{array}{}\text{(52)}& I_t+a\left(y-x\right)I_x+\left(cx-dy-xz\right)I_y+\left(-bz+xy+w\right)I_z+\left(-rw+kz\right)I_w=\mathrm{0}.\end{array}$$Assume that the solution of (52) takes the form$$\begin{array}{}\text{(53)}& I\left(x,y,z,w,t\right)=R\left(x,y,z,w\right)e^{st}.\end{array}$$where $s$ is a constant. Substituting (53) into (52), we obtain$$\begin{array}{}\text{(54)}& sR+a\left(y-x\right)R_x+\left(cx-dy-xz\right)R_y+\left(-bz+xy+w\right)R_z+\left(-rw+kz\right)R_w=\mathrm{0},\end{array}$$Let the solution of (54) take the form$$\begin{array}{}\text{(55)}& R\left(x,y,z,w\right)=A\left(x\right)+B\left(y\right)+T\left(z\right)+Q\left(w\right).\end{array}$$Substituting (55) into (54), we obtain$$\begin{array}{}\text{(56)}& s\left(A\left(x\right)+B\left(y\right)+T\left(z\right)+Q\left(w\right)\right)+a\left(y-x\right)A^{\prime }\left(x\right)+\left(cx-dy-xz\right)B^{\prime }\left(y\right)+\left(-bz+xy+w\right)T^{\prime }\left(z\right)+\left(-rw+kz\right)Q^{\prime }\left(w\right)=\mathrm{0}.\end{array}$$By derivative (56) twice with respect to $w$, we obtain$$\begin{array}{}\text{(57)}& kz{Q}^{\prime \prime \prime }\left(w\right)-rw{Q}^{\prime \prime \prime }\left(w\right)+\left(s-\mathrm{2}r\right)Q^{\prime \prime }\left(w\right)=\mathrm{0}.\end{array}$$Equating the coefficient of $z$ with zero and solving the obtained equation, we obtain$$\begin{array}{}\text{(58)}& Q\left(w\right)=c_{\mathrm{1}}+c_{\mathrm{2}}w+c_{\mathrm{3}}{w}^{\mathrm{2}},\end{array}$$where $c_{\mathrm{1}}$ are arbitrary constant. Substituting (58) into (56), we obtain$$\begin{array}{}\text{(59)}& sT\left(z\right)+xy{T}^{\prime }\left(z\right)+a\left(y-x\right)A^{\prime }\left(x\right)+sA\left(x\right)-bz{T}^{\prime }\left(z\right)+B^{\prime }\left(y\right)\left(cx-dy-xz\right)+sB\left(y\right)+c_{\mathrm{1}}s+w\left(-c_{\mathrm{2}}r+c_{\mathrm{2}}s+\mathrm{2}\mathrm{c}\mathrm{3}kz+T^{\prime }\left(z\right)\right)+c_{\mathrm{2}}kz+w^{\mathrm{2}}\left(c_{\mathrm{3}}s-\mathrm{2}{c}_{\mathrm{3}}r\right)=\mathrm{0}.\end{array}$$Equating the coefficient of $w$ with zero and solving the obtained equation, we obtain$$\begin{array}{}\text{(60)}& T\left(z\right)=c_{\mathrm{2}}\left(r-s\right)z-c_{\mathrm{3}}k{z}^{\mathrm{2}}+c_{\mathrm{4}},\end{array}$$where $c_{\mathrm{4}}$ are arbitrary constant. Substituting (60) into (59), we obtain$$\begin{array}{}\text{(61)}& c_{\mathrm{4}}s+a\left(y-x\right)A^{\prime }\left(x\right)+s\left(A\left(x\right)+B\left(y\right)\right)+z\left(c_{\mathrm{2}}\left(\left(b-s\right)\left(s-r\right)+k\right)\right)+c_{\mathrm{3}}k{z}^{\mathrm{2}}\left(\mathrm{2}b-s\right)+B^{\prime }\left(y\right)\left(cx-dy\right)+xz\left(-B^{\prime }\left(y\right)-\mathrm{2}{c}_{\mathrm{3}}ky\right)+c_{\mathrm{1}}s+c_{\mathrm{2}}xy\left(r-s\right)+c_{\mathrm{3}}{w}^{\mathrm{2}}\left(s-\mathrm{2}r\right)=\mathrm{0}\end{array}$$Equating the coefficient of $xz$ with zero and solving the obtained equation, we obtain$$\begin{array}{}\text{(62)}& B\left(y\right)=c_{\mathrm{5}}-c_{\mathrm{3}}k{y}^{\mathrm{2}},\end{array}$$where $c_{\mathrm{5}}$ are arbitrary constant. Substituting (62) into (61), we obtain$$\begin{array}{}\text{(63)}& y\left(a{A}^{\prime }\left(x\right)+x\left(c_{\mathrm{2}}\left(r-s\right)-\mathrm{2}c{c}_{\mathrm{3}}k\right)\right)-ax{A}^{\prime }\left(x\right)+sA\left(x\right)+c_{\mathrm{2}}z\left(\left(b-s\right)\left(s-r\right)+k\right)+c_{\mathrm{3}}k{z}^{\mathrm{2}}\left(\mathrm{2}b-s\right)+c_{\mathrm{1}}s+c_{\mathrm{3}}k{y}^{\mathrm{2}}\left(\mathrm{2}d-s\right)-\mathrm{2}{c}_{\mathrm{3}}r{w}^{\mathrm{2}}+c_{\mathrm{3}}s{w}^{\mathrm{2}}+c_{\mathrm{4}}s+c_{\mathrm{5}}s=\mathrm{0}.\end{array}$$Equating the coefficient of $y$ with zero and solving the obtained equation, we obtain$$\begin{array}{}\text{(64)}& A\left(x\right)=\frac{\mathrm{1}}{\mathrm{2}a}{x}^{\mathrm{2}}\left(\mathrm{2}c{c}_{\mathrm{3}}k+c_{\mathrm{2}}\left(s-r\right)\right)+c_{\mathrm{6}},\end{array}$$where $c_{\mathrm{6}}$ are arbitrary constant. Substituting (64) into (63), we obtain$$\begin{array}{}\text{(65)}& \frac{\mathrm{1}}{\mathrm{2}a}{x}^{\mathrm{2}}\left(s-\mathrm{2}a\right)\left(\mathrm{2}c{c}_{\mathrm{3}}k+c_{\mathrm{2}}\left(s-r\right)\right)+c_{\mathrm{2}}z\left(\left(b-s\right)\left(s-r\right)+k\right)+c_{\mathrm{3}}k{z}^{\mathrm{2}}\left(\mathrm{2}b-s\right)+s\left(c_{\mathrm{1}}+c_{\mathrm{4}}+c_{\mathrm{5}}+c_{\mathrm{6}}\right)+c_{\mathrm{3}}k{y}^{\mathrm{2}}\left(\mathrm{2}d-s\right)+c_{\mathrm{3}}{w}^{\mathrm{2}}\left(s-\mathrm{2}r\right)\end{array}$$Equation (65) satisfies when

(1)

$s=\mathrm{0},b=c=d=r=\mathrm{0}$, and $c_{\mathrm{2}}=\mathrm{0}$,

(2)

$s=\mathrm{2}b,a=d=r=b,c_{\mathrm{6}}=-c_{\mathrm{1}}-c_{\mathrm{4}}-c_{\mathrm{5}}$, and $c_{\mathrm{2}}=\mathrm{0}$,

(3)

$s=\mathrm{2}b,a=d=r=b,c_{\mathrm{6}}=-c_{\mathrm{1}}-c_{\mathrm{4}}-c_{\mathrm{5}}$, and $k=b^{\mathrm{2}}$,

(4)

$s=\mathrm{2}a,k=(\mathrm{2}a-b)(\mathrm{2}a-r),c_{\mathrm{6}}=-c_{\mathrm{1}}-c_{\mathrm{4}}-c_{\mathrm{5}}$, and $c_{\mathrm{3}}=\mathrm{0}$,

(5)

$s=\mathrm{2}b,d=r=b,c_{\mathrm{6}}=-c_{\mathrm{1}}-c_{\mathrm{4}}-c_{\mathrm{5}},k=b^{\mathrm{2}}$, and $c_{\mathrm{3}}=-c_{\mathrm{2}}/\mathrm{2}bc$.

So, when we substitute (58), (60), (62), and (64) into (53), we obtain the following first integrals of system (1), respectively:

(1)

$I_{\mathrm{1}}(x,y,z,w,t)=w^{\mathrm{2}}-k(y^{\mathrm{2}}+z^{\mathrm{2}})$,

(2)

$I_{\mathrm{2}}(x,y,z,w,t)=e^{\mathrm{2}bt}\left(b{w}^{\mathrm{2}}-kb\left(y^{\mathrm{2}}+z^{\mathrm{2}}\right)+ck{x}^{\mathrm{2}}\right)$,

(3)

$I_{\mathrm{3}}(x,y,z,w,t)=e^{\mathrm{2}bt}\left(\mathrm{2}{c}_{\mathrm{3}}\left(w^{\mathrm{2}}+cb{x}^{\mathrm{2}}-b^{\mathrm{2}}\left(y^{\mathrm{2}}+z^{\mathrm{2}}\right)\right)+c_{\mathrm{2}}\left(-\mathrm{2}bz+\mathrm{2}w+x^{\mathrm{2}}\right)\right)$,

(4)

$I_{\mathrm{4}}(x,y,z,w,t)=e^{\mathrm{2}at}\left(x^{\mathrm{2}}(r-\mathrm{2}a)+\mathrm{2}az(\mathrm{2}a-r)-\mathrm{2}aw\right)$,

(5)

$I_{\mathrm{5}}(x,y,z,w,t)=e^{\mathrm{2}bt}\left(\mathrm{2}cb(bz-w)-b^{\mathrm{2}}\left(y^{\mathrm{2}}+z^{\mathrm{2}}\right)+w^{\mathrm{2}}\right).$

Hence, the proof is completed.

4.1. First Integrals Containing the Origin Equilibrium Point

In this subsection, the invariant surfaces in phase space of hyperchaotic system (1) which contain the origin equilibrium point are attained. If we substitute the equilibrium points into the first integrals of the system (1), we obtain$$\begin{array}{}\text{(66)}& I_{\mathrm{1}}=w^{\mathrm{2}}-k\left(y^{\mathrm{2}}+z^{\mathrm{2}}\right)=\mathrm{0},\text{(67)}& I_{\mathrm{2}}=y^{\mathrm{2}}+z^{\mathrm{2}}-\frac{c}{b}{x}^{\mathrm{2}}-\frac{\mathrm{1}}{k}{w}^{\mathrm{2}}=\mathrm{0},\text{(68)}& I_{\mathrm{3}}=\mathrm{2}{c}_{\mathrm{3}}\left(w^{\mathrm{2}}+cb{x}^{\mathrm{2}}-b^{\mathrm{2}}\left(y^{\mathrm{2}}+z^{\mathrm{2}}\right)\right)+c_{\mathrm{2}}\left(-\mathrm{2}bz+\mathrm{2}w+x^{\mathrm{2}}\right)=\mathrm{0},\text{(69)}& I_{\mathrm{4}}=x^{\mathrm{2}}\left(r-\mathrm{2}a\right)+\mathrm{2}az\left(\mathrm{2}a-r\right)-\mathrm{2}aw=\mathrm{0},\text{(70)}& I_{\mathrm{5}}=\mathrm{2}cb\left(bz-w\right)-b^{\mathrm{2}}\left(y^{\mathrm{2}}+z^{\mathrm{2}}\right)+w^{\mathrm{2}}=\mathrm{0}.\end{array}$$Hence, we will classify the solutions (66)-(70) as follows:

(1)

For (66), when $k>\mathrm{0}$, (66) becomes circular cone (see Figure 5).

(2)

For (67), there are some cases:

(i)

when $c=\mathrm{0}$ and $k>\mathrm{0}$, (67) becomes circular cone (see Figure 6(a)),

(ii)

when $c/b>\mathrm{0}$, $k>\mathrm{0}$, and $w=\beta$, (67) becomes hyperboloid of one sheet (see Figure 6(b)),

(iii)

when $c/b<\mathrm{0}$, $k>\mathrm{0}$, and $w=\beta$, (67) becomes ellipsoid (see Figure 6(c)),

(iv)

when $c=b$, $k>\mathrm{0}$, and $w=\beta$, (67) becomes sphere (see Figure 6(d)),

(v)

when $c/b>\mathrm{0}$, $k>\mathrm{0}$, and $w=\beta$, (67) becomes hyperboloid of two sheets (see Figure 6(e)).

(3)

Let $c_{\mathrm{2}}=\mathrm{2}{c}_{\mathrm{3}}$; (68) becomes$$\begin{array}{}\text{(71)}& y^{\mathrm{2}}+{\left(\frac{\mathrm{1}}{b}+z\right)}^{\mathrm{2}}-\frac{\left(bc+\mathrm{1}\right)}{b^{\mathrm{2}}}{x}^{\mathrm{2}}-\frac{\mathrm{1}}{b^{\mathrm{2}}}{\left(w+\mathrm{1}\right)}^{\mathrm{2}}=\mathrm{0}.\end{array}$$

 

which is classified as follows:

(i)

when $bc=-\mathrm{1}$, (71) becomes circular cone (see Figure 7(a)),

(ii)

when $bc>-\mathrm{1}$ and $w=\beta$, (71) becomes hyperboloid of one sheet (see Figure 7(b)),

(iii)

when $bc<-\mathrm{1}$ and $w=\beta$, (71) becomes ellipsoid (see Figure 7(c)),

(iv)

when $c=-b-\mathrm{1}/b$ and $w=\beta$, (71) becomes sphere (see Figure 7(d)).

(4)

Equation (69) represents parabola (see Figure 8).

(5)

Equation (70)can be rewritten in the form$$\begin{array}{}\text{(72)}& y^{\mathrm{2}}+{\left(z-c\right)}^{\mathrm{2}}-\frac{{\left(w-bc\right)}^{\mathrm{2}}}{b^{\mathrm{2}}}=\mathrm{0},\end{array}$$

 

which represents circular cone (see Figure 9).

[figure omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

$\beta$ is a constant.

5. Circuit Implementation

The circuit implementation of hyperchaotic system (1) can be designed using operational amplifier, analog multipliers, resistors, and capacitors [17, 18]. This hyperchaotic nonlinear circuit is equivalent to system (1) and enables the useful experimental study of complex dynamics exhibited by system (1) and verify the results of theoretical analysis carried out upon it. The Multisim package is used to carry out circuit simulations and depict the circuit outputs on the oscilloscope. The nonlinear circuit implementation of hyperchaotic system can be used as a central in a variety of advanced applications, e.g., chaos based physical-layer secure communications and physical random number generators which outperform pseudo random number generators. Figure 10 shows schematics of proposed the circuit realization of hyperchaotic system (1) for the values of parameters given by $a=\mathrm{21.7},b=\mathrm{7.3},c=\mathrm{6.6},d=-\mathrm{2},r=\mathrm{0.1}$, and $k=\mathrm{9.5}.$

[figures omitted; refer to PDF]

Regarding digital realization of chaotic circuits, it is important to note that dynamics degradation of digital chaos occurs due to impossibility of the chaotic map to reach an ideal chaotic state under the limitations of digital domain with finite-precision. Various methods are proposed to counteract dynamics degradation; see, for example, [42] and references therein.

6. Conclusion

In this work, dynamical and bifurcation analysis are applied to explore local and global characteristics of 4D Lorenz-type hyperchaotic system. Through investigation of codimension two local bifurcation of the hyperchaotic system, it is found that the 4D hyperchaotic system undergoes various types of bifurcations including BT bifurcation, Hopf bifurcation, pitchfork bifurcation, and homoclinic bifurcation. Moreover, the paper successfully acquires, for first time, the closed forms of first integrals and regular invariant surfaces in solution space, which pass through the origin equilibrium points, and relates these forms to the parameters of the system. Thus, it predicts some analytical solutions of the system and estimates the associated parameters values. Finally, a circuit realization for 4D hyperchaotic system is proposed as an engineering application and numerical examples are presented to verify theoretical results. Future work can include analytical and numerical investigations of hidden attractors, if exist, in hyperchaotic 4D Lorenz system. Also, other possible first integrals and invariant structures of the state space of the system may be explored in future work.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to extend their sincere appreciation to the Deanship of Scientific Research at King Saud University for funding this research (Group no. RG-1438-046).

References

[1] E. M. Izhikevich, Dynamical Systems in Neuroscience: the Geometry of Excitability And Bursting, 2007.

[2] P. Stavroulakis, Chaos Applications in Telecommunications,DOI: 10.1201/9780203025314, 2006.

[3] S. H. Strogatz, Nonlinear dynamics and chaos with applications to physics. Biology, chemistry, and engineering, 2001.

[4] P. N. V. Tu, Dynamical Systems: An Introduction with Applications in Economics And Biology,DOI: 10.1007/978-3-662-02779-0, 2012.

[5] M. Han, P. Yu, Normal Forms, Melnikov Functions and Bifurcations of Limit Cycles, vol. 181,DOI: 10.1007/978-1-4471-2918-9, 2012.

[6] J. D. Meiss, Differential Dynamical Systems, vol. 14,DOI: 10.1137/1.9780898718232, 2007.

[7] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, vol. 2, 2003.

[8] A. Elsonbaty, S. F. Hegazy, S. S. A. Obayya, "Simultaneous concealment of time delay signature in chaotic nanolaser with hybrid feedback," Optics and Lasers in Engineering, vol. 107, pp. 342-351, DOI: 10.1016/j.optlaseng.2018.03.027, 2018.

[9] A. Elsonbaty, S. F. Hegazy, S. S. A. Obayya, "Numerical analysis of ultrafast physical random number generator using dual-channel optical chaos," Optical Engineering, vol. 55 no. 9, 2016.

[10] A. Elsonbaty, S. F. Hegazy, S. S. A. Obayya, "Simultaneous suppression of time-delay signature in intensity and phase of dual-channel chaos communication," IEEE Journal of Quantum Electronics, vol. 51,DOI: 10.1109/JQE.2015.2466176, 2015.

[11] L. Kocarev, S. Lian, Chaos-Based Cryptography,DOI: 10.1007/978-3-642-20542-2, 2011.

[12] A. A. Elsadany, A. M. Yousef, A. Elsonbaty, "Further analytical bifurcation analysis and applications of coupled logistic maps," Applied Mathematics and Computation, vol. 338, pp. 314-336, DOI: 10.1016/j.amc.2018.06.008, 2018.

[13] A. A. Elsadany, A. Elsonbaty, H. N. Agiza, "Qualitative dynamical analysis of chaotic plasma perturbations model," Communications in Nonlinear Science and Numerical Simulation, vol. 59, pp. 409-423, DOI: 10.1016/j.cnsns.2017.11.020, 2018.

[14] A. Al-khedhairi, A. A. Elsadany, A. Elsonbaty, A. G. Abdelwahab, "Dynamical study of a chaotic predator-prey model with an omnivore," The European Physical Journal Plus, vol. 133 no. 29, 2018.

[15] A. Elsonbaty, A. A. Elsadany, "Bifurcation analysis of chaotic geomagnetic field model," Chaos, Solitons & Fractals, vol. 103, pp. 325-335, DOI: 10.1016/j.chaos.2017.06.022, 2017.

[16] A. Elsonbaty, A. M. El-Sayed, "Analytical study of global bifurcations, stabilization and chaos synchronization of jerk system with multiple attractors," Nonlinear Dynamics, vol. 90 no. 4, pp. 2637-2655, DOI: 10.1007/s11071-017-3828-z, 2017.

[17] A. R. Elsonbaty, A. M. El-Sayed, "Further nonlinear dynamical analysis of simple jerk system with multiple attractors," Nonlinear Dynamics, vol. 87 no. 2, pp. 1169-1186, DOI: 10.1007/s11071-016-3108-3, 2017.

[18] A. M. El-Sayed, H. M. Nour, A. Elsaid, A. E. Matouk, A. Elsonbaty, "Dynamical behaviors, circuit realization, chaos control, and synchronization of a new fractional order hyperchaotic system," Applied Mathematical Modelling: Simulation and Computation for Engineering and Environmental Systems, vol. 40 no. 5-6, pp. 3516-3534, DOI: 10.1016/j.apm.2015.10.010, 2016.

[19] Q. Wang, S. Yu, C. Li, "Theoretical design and FPGA-based implementation of higher-dimensional digital chaotic systems," IEEE Transactions on Circuits and Systems I-Regular Papers, vol. 63 no. 3, pp. 401-412, DOI: 10.1109/tcsi.2016.2515398, 2016.

[20] Q. Hong, Q. Xie, Y. Shen, X. Wang, "Generating multi-double-scroll attractors via nonautonomous approach," Chaos: An Interdisciplinary Journal of Nonlinear Science, vol. 26 no. 8,DOI: 10.1063/1.4959538, 2016.

[21] C. Li, "Cracking a hierarchical chaotic image encryption algorithm based on permutation," Signal Processing, vol. 118, pp. 203-210, DOI: 10.1016/j.sigpro.2015.07.008, 2016.

[22] E. Harjanto, J. M. Tuwankotta, "Bifurcation of periodic solution in a Predator–Prey type of systems with non-monotonic response function and periodic perturbation," International Journal of Non-Linear Mechanics, vol. 85, pp. 188-196, DOI: 10.1016/j.ijnonlinmec.2016.06.011, 2016.

[23] L. Barreira, C. Valls, J. Llibre, "Integrability and limit cycles of the Moon-Rand system," International Journal of Non-Linear Mechanics, vol. 69, pp. 129-136, DOI: 10.1016/j.ijnonlinmec.2014.11.029, 2015.

[24] J. Kengne, "Coexistence of chaos with hyperchaos, period-3 doubling bifurcation, and transient chaos in the hyperchaotic oscillator with gyrators," International Journal of Bifurcation and Chaos, vol. 25 no. 4,DOI: 10.1142/S0218127415500522, 2015.

[25] A. P. Kuznetsov, S. P. Kuznetsov, E. Mosekilde, N. V. Stankevich, "Co-existing hidden attractors in a radio-physical oscillator system," Journal of Physics A: Mathematical and Theoretical, vol. 48 no. 12,DOI: 10.1088/1751-8113/48/12/125101, 2015.

[26] J. Kengne, J. C. Chedjou, T. Fonzin Fozin, K. Kyamakya, G. Kenne, "On the analysis of semiconductor diode-based chaotic and hyperchaotic generators---a case study," Nonlinear Dynamics, vol. 77 no. 1-2, pp. 373-386, DOI: 10.1007/s11071-014-1301-9, 2014.

[27] A. M. El-Sayed, H. M. Nour, A. Elsaid, A. E. Matouk, A. Elsonbaty, "Circuit realization, bifurcations, chaos and hyperchaos in a new 4D system," Applied Mathematics and Computation, vol. 239, pp. 333-345, DOI: 10.1016/j.amc.2014.04.109, 2014.

[28] H. M. Nour, A. Elsaid, A. Elsonbaty, "Circuit realization, chaos synchronization and estimation of parameters of a hyperchaotic system with unknown parameters," Journal of the Egyptian Mathematical Society, vol. 22 no. 3, pp. 550-557, DOI: 10.1016/j.joems.2013.11.007, 2014.

[29] A. M. El-Sayed, A. Elsaid, H. M. Nour, A. Elsonbaty, "Synchronization of different dimensional chaotic systems with time varying parameters, disturbances and input nonlinearities," Journal of Applied Analysis and Computation, vol. 4 no. 4, pp. 323-338, 2014.

[30] A. M. El-Sayed, A. Elsaid, H. M. Nour, A. Elsonbaty, "Dynamical behavior, chaos control and synchronization of a memristor-based ADVP circuit," Communications in Nonlinear Science and Numerical Simulation, vol. 18 no. 1, pp. 148-170, DOI: 10.1016/j.cnsns.2012.06.011, 2013.

[31] C. Li, Y. Liu, T. Xie, M. Z. Q. Chen, "Breaking a novel image encryption scheme based on improved hyperchaotic sequences," Nonlinear Dynamics, vol. 73 no. 3, pp. 2083-2089, DOI: 10.1007/s11071-013-0924-6, 2013.

[32] B. Muthuswamy, L. O. Chua, "Simplest chaotic circuit," International Journal of Bifurcation and Chaos, vol. 20 no. 5, pp. 1567-1580, DOI: 10.1142/S0218127410027076, 2010.

[33] A. H. Nayfeh, B. Balachandran, Applied Nonlinear Dynamics: Analytical, Computational and Experimental Methods,DOI: 10.1002/9783527617548, 2008.

[34] N. A. A. Fataf, S. K. Palit, S. Mukherjee, M. R. M. Said, D. H. Son, S. Banerjee, "Communication scheme using a hyperchaotic semiconductor laser model: Chaos shift key revisited," The European Physical Journal Plus, vol. 132 no. 492, 2017.

[35] S. Vaidyanathan, A. Akgul, S. Kaçar, U. Çavuşoğlu, "A new 4-D chaotic hyperjerk system, its synchronization, circuit design and applications in RNG, image encryption and chaos-based steganography," The European Physical Journal Plus, vol. 133 no. 46, 2018.

[36] Y. Chen, "The existence of homoclinic orbits in a 4D Lorenz-type hyperchaotic system," Nonlinear Dynamics, vol. 87 no. 3, pp. 1445-1452, DOI: 10.1007/s11071-016-3126-1, 2017.

[37] G. A. Leonov, N. V. Kuznetsov, V. I. Vagaitsev, "Hidden attractor in smooth Chua systems," Physica D: Nonlinear Phenomena, vol. 241 no. 18, pp. 1482-1486, DOI: 10.1016/j.physd.2012.05.016, 2012.

[38] C. Li, J. C. Sprott, "Coexisting hidden attractors in a 4-D simplified lorenz system," International Journal of Bifurcation and Chaos, vol. 24 no. 3,DOI: 10.1142/S0218127414500345, 2014.

[39] G. Peng, Y. Jiang, "Computation of universal unfolding of the double-zero bifurcation in symmetric systems by a homological method," Journal of Difference Equations and Applications, vol. 19 no. 9, pp. 1501-1512, DOI: 10.1080/10236198.2012.761980, 2013.

[40] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory, 1997.

[41] J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcation of Vector Fields, 1983.

[42] C. Li, B. Feng, S. Li, J. Kurths, G. Chen, "Dynamic analysis of digital chaotic maps via state-mapping networks," IEEE Transactions on Circuits and Systems I: Regular Papers,DOI: 10.1109/TCSI.2018.2888688, 2018.