Local Dynamics Analysis of a Four-Dimensional Hyperchaotic Lorenz System

Yu Ren^1,2, Liangqiang Zhou^1,2*

¹School of Mathematics, Nanjing University of Aeronautics and Astronautics, Nanjing, China.
²Key Laboratory of Mathematical Modelling and High Performance Computing of Air Vehicles (NUAA), Nanjing, China.

DOI: 10.4236/jamp.2023.1110182                

Abstract

The local dynamical behaviors of a four-dimensional hyperchaotic Lorenz system, including stability and bifurcations, are investigated in this paper by analytical and numerical methods. The equilibriums and their stability under different parameter conditions are analyzed by applying Routh-Hurwitz criterion. The results indicate that the system may exist one, three and five equilibrium points for different system parameters. Based on the central manifold theorem and normal form theorem, the pitchfork bifurcation and Hopf bifurcation are studied respectively. By using the Hopf bifurcation theorem and calculating the first Lyapunov coefficient, the Hopf bifurcation of this system is obtained as supercritical for certain parameters. Finally, the results of theoretical parts are verified by some numerical simulations.

Keywords

Hyperchaotic Lorenz System, Stability, Pitchfork Bifurcation, Hopf Bifurcation, Central Manifold

1. Introduction

In 1963, Lorenz, an American meteorologist, discovered the first chaotic attractor while studying the atmospheric motion, which started the study of chaos. The Lorenz model also has become a classical model in the field of chaos.

Based on Lorenz system, many scholars proposed Lorenz-like systems, which are widely used in image encryption, secure communication and many other fields. Perevoznikov et al. [1] considered two methods of stability analysis of systems described by dynamical equations through Lorenz and Rossler Model Problems. Liu et al. [2] gave new insights into a new Lorenz-like chaotic system from the literature. Wu et al. [3] investigated the stability of equilibrium points and Hopf bifurcation of a Lorenz-like system with the normal form theory. The conditions guaranteeing the Hopf bifurcation were derived. Hu et al. [4] studied stability and Hopf bifurcation of a new four-dimensional quadratic autonomous system both analytically and numerically. Aziz et al. [5] designed a new electronic circuit as an engineering application on a four-dimensional chaotic system. Fu et al. [6] combined a three-dimensional autonomous system with a complex network, and used the chaotic synchronization of this system to apply it to confidential communication. Wang et al. [7] suggested a brand-new Julia-fractal system in three-dimensional to increase the randomness of chaotic sequences and widen the key space and then a fresh quantum circuit for Fibonacci scrambling was created.

The bifurcation problem also has a profound application background. Zhang et al. [8] investigated the bifurcation characteristics of the active magnetic bearing-rotor system subjected to the external excitation analytically. Wang et al. [9] presented a two-dimensional simplified Hodgkin-Huxley model under exposure to electric fields, and the Hopf bifurcation of this model was analyzed both analytically and numerically. Wen et al. [10] obtained the Hopf bifurcation associated with positive non-homogeneous steady states of a reaction-diffusion-advection logistic model, with two non-local delayed density-dependent terms and zero-Dirichlet boundary conditions. Li et al. [11] studied a four-dimensional inertial two-nervous system with delay and obtained the critical value of zero-Hopf bifurcation by analyzing the distribution of eigenvalues.

The following manuscript is organized as follows: a four-dimensional hyperchaotic Lorenz system is presented in the next section. In Section 3, the equilibrium points and their stability are studied. The pitchfork bifurcation and Hopf bifurcation are analyzed in Section 4. In the last section, numerical simulations are presented to verify the theoretical results.

2. Dynamic Modeling

Wang et al. [12] proposed a four-dimensional hyperchaotic system by adding the nonlinear controller to Lorenz system, and analyzed its basic dynamic characteristics. On the basis of these works, Xie et al. [13] proposed a plaintext related color image encryption algorithm based on multiple chaotic maps to solve the problem of color image encryption and the local dynamics analysis of this system is carried out. The specific model is shown in the following:

$\{\begin{array}{l}\stackrel{\dot{}}{x}=a\left(y-x\right)+w\hfill \\ \stackrel{\dot{}}{y}=bx-y-xz\hfill \\ \stackrel{\dot{}}{z}=xy-cz\hfill \\ \stackrel{\dot{}}{w}=-yz+dw\hfill \end{array}$ (1)

Where ${\left[x\mathrm{,}y\mathrm{,}z\mathrm{,}w\right]}^{\text{T}}$ are the state variables, $a\mathrm{,}b\mathrm{,}c\mathrm{,}d$ are controllable parameters. When $a=10$ , $b=28$ , $c=\frac{8}{3}$ , $-1.52\le d\le -0.06$ , the system is in a hyperchaotic state [13] . Its hyperchaotic state can be observed by Figure 1.

Figure 1. The projections of system (1) when $a=10$ , $b=28$ , $c=\frac{8}{3}$ , $d=-1$ . (a) The three-dimensional projection of x-z-w plane; (b) The two-dimensional projection of x-y plane; (c) The two-dimensional projection of y-z plane; (d) The two-dimensional projection of x-z plane.

3. Equilibrium Analysis

In this section, the stability of the equilibrium points of this system will be analyzed.

Theorem 1. Let $b\ne 0$ , $c\ne 0$ , $\Delta =b^2{c}^2\left(b^2+a^2{c}^2{d}^2+2ab{c}^{2}d-4a{c}^{2}d\right)$ , $t_1=\frac{\left(-2acd+abcd+b^{2}c\right)+\sqrt{\Delta }}{2ad}$ , $t_2=\frac{\left(-2acd+abcd+b^{2}c\right)-\sqrt{\Delta }}{2ad}$ , the system has only one equilibrium point $S_0=\left(\mathrm{0,0,0,0}\right)$ as long as it meets any of the following conditions: 1) $\Delta <0$ i.e. $b^2+a^2{c}^2{d}^2+2ab{c}^{2}d-4a{c}^{2}d<0$ ; 2) $\Delta =0$ and also $\frac{-2acd+abcd+b^{2}c}{2ad}<0$ ; 3) $\Delta >0$ and also $t_1,t_2<0$ . And the system has three equilibrium points as long as it meets any of the following conditions: 1) $\Delta =0$ ; 2) $\Delta >0$ and also $t_1>0$ , $t_2<0$ ; 3) $\Delta >0$ and also $t_1<0$ , $t_2>0$ . The system has five equilibrium points as long as it meets the condition: both $\Delta >0$ and $t_1>0$ , $t_2>0$ .

Proof 1. Solving the equations to find the equilibrium points of system (1):

$\{\begin{array}{l}a\left(y-x\right)+w=0\hfill \\ bx-y-xz=0\hfill \\ xy-cz=0\hfill \\ -yz+dw=0\hfill \end{array}$ (2)

when $x=0$ , it’s easy to obtain the equilibrium point $S_0=\left(\mathrm{0,0,0,0}\right)$ ; and when $x\ne 0$ ,we have

$ad{\left(x^2+c\right)}^2-b^{2}c{x}^2-abcd\left(x^2+c\right)=0$ (3)

let $t_1=\frac{\left(-2acd+abcd+b^{2}c\right)+\sqrt{\Delta }}{2ad}$ , $t_2=\frac{\left(-2acd+abcd+b^{2}c\right)-\sqrt{\Delta }}{2ad}$ , where $\Delta =b^2{c}^2\left(b^2+a^2{c}^2{d}^2+2ab{c}^{2}d-4a{c}^{2}d\right)$ .

When $\Delta <0$ the system has only one equilibrium point S₀.

When $\Delta =0$ and also $\frac{-2acd+abcd+b^{2}c}{2ad}<0$ , the system has equilibrium point S₀. And when $\Delta =0$ and the $\frac{-2acd+abcd+b^{2}c}{2ad}>0$ , the system has three equilibrium points as following: S₀; $S_1=\left(x_1\mathrm{,}\frac{bc{x}_1}{x_1^2+c}\mathrm{,}\frac{b{x}_1^2}{x_1^2+c}\mathrm{,}a{x}_1-\frac{abc{x}_1}{x_1^2+c}\right)$ ; $S_2=\left(x_2\mathrm{,}\frac{bc{x}_2}{x_2^2+c}\mathrm{,}\frac{b{x}_2^2}{x_2^2+c}\mathrm{,}a{x}_2-\frac{abc{x}_2}{x_2^2+c}\right)$ . Where $x_1=\sqrt{\frac{-2acd+abcd+b^{2}c}{2ad}}$ , $x_2=-\sqrt{\frac{-2acd+abcd+b^{2}c}{2ad}}$ .

When $\Delta >0,t_1>0,t_2<0$ or $\Delta >0,t_1<0,t_2>0$ , the system has three equilibrium points. The equilibrium points in the former case are $S_0$ ; $M_1=\left(m_1,\frac{bc{m}_1}{m_1^2+c},\frac{b{m}_1^2}{m_1^2+c},a{m}_1-\frac{abc{m}_1}{m_1^2+c}\right)$ , $m_1=\sqrt{t_1}$ ; $M_2=\left(m_2\mathrm{,}\frac{bc{m}_2}{m_2^2+c}\mathrm{,}\frac{b{m}_2^2}{m_2^2+c}\mathrm{,}a{m}_2-\frac{abc{m}_2}{m_2^2+c}\right)$ , $m_2=-\sqrt{t_1}$ respectively. And the equilibrium points in the other case are S₀; $N_1=\left(n_1\mathrm{,}\frac{bc{n}_1}{n_1^2+c}\mathrm{,}\frac{b{n}_1^2}{n_1^2+c}\mathrm{,}a{n}_1-\frac{abc{n}_1}{n_1^2+c}\right)$ , $n_1=\sqrt{t_2}$ ; $N_2=\left(n_2\mathrm{,}\frac{bc{n}_2}{n_2^2+c}\mathrm{,}\frac{b{n}_2^2}{n_2^2+c}\mathrm{,}a{n}_2-\frac{abc{n}_2}{n_2^2+c}\right)$ , $n_2=-\sqrt{t_2}$ respectively.

When $\Delta >0,t_1>0,t_2>0$ , the system has five equilibrium points. The equilibrium points are S₀, M₁, M₂, N₁, N₂. Therefore, Theorem 1 is proved.

Then the stability will be discussed. Choosing case $\Delta >0,t_1>0,t_2>0$ as an example and the analysis process for other situations are similar. Based on Theorem 1, the system has equilibrium points S₀, M₁, M₂, N₁, N₂ when $\Delta >0,t_1>0,t_2>0$ . Since the system has symmetry, the stability of M₂, N₂ are same as the M₁, N₁. The Routh-Hurwitz criterion is applied to analyze stability of S₀, M₁, N₁, and the following theorems are obtained.

Theorem 2. When ${\left(a+1\right)}^2-4\left(a-ab\right)<0$ , $c>0$ , $d<0$ , $a>-1$ ; or ${\left(a+1\right)}^2-4\left(a-ab\right)\ge 0$ and $c>0$ , $d<0$ , $\frac{-\left(a+1\right)\pm \sqrt{{\left(a+1\right)}^2-4\left(a-ab\right)}}{2}<0$ the equilibrium point $S_0=\left(\mathrm{0,0,0,0}\right)$ is asymptotic stable.

Proof 2. The Jacobian matrix of system (1) at S₀ is as following:

$J_{S_0}=\left(\begin{array}{cccc}-a& a& 0& 1\\ b& -1& 0& 0\\ 0& 0& -c& 0\\ 0& 0& 0& d\end{array}\right)$

It’s characteristic polynomial is $P_1\left(\lambda \right)=\left(\lambda +c\right)\left(\lambda -d\right)\left[{\lambda }^2+\left(a+1\right)\lambda +\left(a-ab\right)\right]=0$ . And characteristic roots are ${\lambda }_1=-c$ , ${\lambda }_2=d$ , ${\lambda }_{3,4}=\frac{-\left(a+1\right)\pm \sqrt{{\left(a+1\right)}^2-4\left(a-ab\right)}}{2}$ .

All characteristic roots have negative real parts if these parameters satisfy ${\left(a+1\right)}^2-4\left(a-ab\right)<0$ , $c>0$ , $d<0$ , $a>-1$ . Applying the Routh-Hurwitz criterion [14] we can know S₀ is asymptotically stable. Similarly, when the parameters simultaneously satisfy ${\left(a+1\right)}^2-4\left(a-ab\right)\ge 0$ and $c>0$ , $d<0$ , $\frac{-\left(a+1\right)\pm \sqrt{{\left(a+1\right)}^2-4\left(a-ab\right)}}{2}<0$ , the equilibrium point is also asymptotically stable.

Theorem 3. When $a_1>0$ , $a_1{a}_2-a_3>0$ , $a_1{a}_2{a}_3-a_1^2{a}_4-a_3^2>0$ , $a_4\left(a_1{a}_2{a}_3-a_1^2{a}_4-a_3^2\right)>0$ , the equilibrium point $M_1=\left(m_1\mathrm{,}\frac{bc{m}_1}{m_1^2+c}\mathrm{,}\frac{b{m}_1^2}{m_1^2+c}\mathrm{,}a{m}_1-\frac{abc{m}_1}{m_1^2+c}\right)$ is asymptotically stable.

Proof 3. The Jacobian matrix of system (1) at M₁ is as following:

$J_{M_1}=\left(\begin{array}{cccc}-a& a& 0& 1\\ b-\frac{b{t}_1}{t_1+c}& -1& -\sqrt{t_1}& 0\\ \frac{bc\sqrt{t_1}}{t_1+c}& \sqrt{t_1}& -c& 0\\ 0& -\frac{b{t}_1}{t_1+c}& -\frac{bc\sqrt{t_1}}{t_1+c}& d\end{array}\right)$

The characteristic polynomial of this matrix is

$P_2\left(\lambda \right)={\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4$

where

$\begin{array}{l}{a}_0=1,a_1=a+c-d+1,a_2=a+c+ac+t_1-\frac{abc}{t_1+c}-ad-cd-d,\\ a_3=ac+a{t}_1-\frac{ab{c}^2}{t_1+c}+\frac{abc{t}_1}{t_1+c}-ad-cd-acd-d{t}_1+\frac{abcd}{t_1+c}+\frac{b^{2}c{t}_1+b^2{c}^2{t}_1}{{(t_1+c)}^2},\\ a_4=\frac{3b^2{c}^2{t}_1-b^{2}c{t}_1^2}{{\left(t_1+c\right)}^2}-acd-ad{t}_1+\frac{ab{c}^{2}d-abcd{t}_1}{t_1+c}.\end{array}$

Calculating the following determinants: ${\Delta }_1>0$ , ${\Delta }_2>0$ , ${\Delta }_3>0$ , ${\Delta }_4>0$ . According to Routh-Hurwitz critetion, the equilibrium point M₁ is asymptotically stable.

The analysis process of equilibrium N₁ is the same as that of M₁. Simply replace the value of t₁ with t₂ in the coordinate of equation M₁, then the stability of point N₁ can be obtained.

4. Bifurcation Analysis

4.1. Pitchfork Bifurcation

In this section, the pitchfork bifurcation of equilibrium point S₀ is investigated. b is set to the bifurcation parameter, which leads to the following theorem.

Theorem 4. Assuming the parameters satisfy conditions that $a\ne -1$ , $c\ne 0$ , $d\ne 0$ , the pitchfork bifurcation occurs when the parameter b crosses the critical value $b=1$ .

Proof 4. The Jacobian matrix of S₀ is

$J_{S_0{|}_{b=1}}=\left(\begin{array}{cccc}-a& a& 0& 1\\ 1& -1& 0& 0\\ 0& 0& -c& 0\\ 0& 0& 0& d\end{array}\right)$

It’s characteristic polynomial is $P_3\left(\lambda \right)=\left(\lambda +c\right)\left(\lambda -d\right)\left[{\lambda }^2+\left(a+1\right)\lambda \right]=0$ and the characteristic roots are ${\lambda }_1=0$ , ${\lambda }_2=-\left(a+1\right)$ , ${\lambda }_3=-c$ , ${\lambda }_4=d$ , respectively. In other words, the system has one zero eigenvalue and the others with non-zero real parts. The corresponding eigenvectors of these four eigenvalues are calculated as follows:

${\xi }_1=\left(\begin{array}{c}1\\ 1\\ 0\\ 0\end{array}\right)$ , ${\xi }_2=\left(\begin{array}{c}-a\\ 1\\ 0\\ 0\end{array}\right)$ , ${\xi }_3=\left(\begin{array}{c}0\\ 0\\ 1\\ 0\end{array}\right)$ , ${\xi }_4=\left(\begin{array}{c}d+1\\ 1\\ 0\\ d^2+d+ad\end{array}\right)$

Let

$T=\left({\xi }_1,{\xi }_2,{\xi }_3,{\xi }_4\right)=\left(\begin{array}{cccc}1& -a& 0& d+1\\ 1& 1& 0& 1\\ 0& 0& 1& 0\\ 0& 0& 0& d^2+d+ad\end{array}\right)$

Using the following transformation

$\left(\begin{array}{c}x\\ y\\ z\\ w\end{array}\right)=T\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right)$

Then system (1) becomes

$\left(\begin{array}{c}{\stackrel{\dot{}}{x}}_1\\ {\stackrel{\dot{}}{x}}_2\\ {\stackrel{\dot{}}{x}}_3\\ {\stackrel{\dot{}}{x}}_4\end{array}\right)=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& -\left(a+1\right)& 0& 0\\ 0& 0& -c& 0\\ 0& 0& 0& d\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right)+\left(\begin{array}{c}{g}_1\\ g_2\\ g_3\\ g_4\end{array}\right)$ (4)

where

$\begin{array}{c}{g}_1=\frac{1}{1+a}\left[\left(a^2+a\right)x_2+\left(a^2+d\right)x_4\right]+\frac{a}{1+a}[\left(b-1\right)x_1-\left(ab+1\right)x_2\\ +\left(bd+b-1\right)x_4-x_1{x}_3+a{x}_2{x}_3-\left(d+1\right)x_3{x}_4]\\ -\frac{a+d+1}{\mathrm{<mo\; stretchy="false">\; (\; </mo>}a+\mathrm{1\; <mo\; stretchy="false">\; )\; </mo>\; <mo\; stretchy="false">\; (\; </mo>}{d}^2+ad+d\mathrm{<mo\; stretchy="false">\; )\; </mo>}}\left[-x_1{x}_3-x_2{x}_3-x_3{x}_4+\left(d^3+a{d}^2+d^2\right)x_4\right]\mathrm{,}\end{array}$

$\begin{array}{c}{g}_2=-\frac{1}{1+a}\left[\left(a^2+d\right)x_4\right]+\frac{1}{1+a}[\left(b-1\right)x_1-\left(ab+1\right)x_2\\ +\left(bd+b-1\right)x_4-x_1{x}_3+a{x}_2{x}_3-\left(d+1\right)x_3{x}_4]\\ +\frac{d}{\left(a+1\right)\left(d^2+ad+d\right)}\left[-x_1{x}_3-x_2{x}_3-x_3{x}_4+\left(d^3+a{d}^2+d^2\right)x_4\right],\end{array}$

$g_3=x_1^2-a{x}_2^2+\left(d+1\right)x_4^2+\left(1-a\right)x_1{x}_2+\left(d+2\right)x_1{x}_4+\left(d+1-a\right)x_2{x}_4\mathrm{,}$

$g_4=\frac{1}{d^2+ad+d}\left(-x_1{x}_3-x_2{x}_3-x_3{x}_4\right)\mathrm{.}$

Next, the central manifold theorem [15] [16] is used to analyze the stability of point S₀ when $b=1$ . Constructing the following local central manifold represented by the two-parameter family x₁ and b. $\delta$ and $\bar{\delta }$ are sufficiently small.

$\begin{array}{l}{W}^c\left(S_0\right)=\{\left(x_1,x_2,x_3,x_4,b\right)\in {ℝ}^5|x_2=h_1\left(x_1,b\right),x_3=h_2\left(x_1,b\right),x_4=h_3\left(x_1,b\right),\\ \left|x_1\right|<\delta ,\left|b\right|<\bar{\delta },h_i\left(0,0\right)=0,D{h}_i\left(0,0\right)=0,i=1,2,3\}\end{array}$

Suppose $h_1\mathrm{,}{h}_2\mathrm{,}{h}_3$ have expressions of the following forms:

$x_2=h_1\left(x_1,b\right)=a_1{x}_1+a_2{x}_{1}b+a_{3}b+\cdots$ (5)

$x_3=h_2\left(x_1,b\right)=b_1{x}_1^2+b_2{x}_{1}b+b_3{b}^2+\cdots$ (6)

$x_4=h_3\left(x_1,b\right)=c_1{x}_1^3+c_2{x}_1^{2}b+c_3{x}_1{b}^2+\cdots$ (7)

According to the chain rule for derivatives, we have

$Dh{g}_1=Bh+g$ (8)

where

$h=\left(\begin{array}{c}{h}_1\\ h_2\\ h_3\end{array}\right),g=\left(\begin{array}{c}{g}_2\\ g_3\\ g_4\end{array}\right),B=\left(\begin{array}{ccc}-\left(a+1\right)& 0& 0\\ 0& -c& 0\\ 0& 0& d\end{array}\right)$

Substituting Equations (4)-(6) to Equation (7), and comparing the coefficients of the corresponding terms with equal exponents, we can calculate the expressions of $h_i$ as follows:

$x_2=-\frac{1}{a^2+a+1}{x}_1+\frac{{\left(a+1\right)}^2}{{\left(a^2+a+1\right)}^2}{x}_{1}b+\cdots$ (9)

$x_3=\frac{a+1}{ac-2a+c}{x}_1^2+\cdots$ (10)

$x_4=\frac{{\left(a+1\right)}^2}{\left(d^2+ad+d\right)\left(ac-2a+c\right)\left(3a+ad+d\right)}{x}_1^3+\cdots$ (11)

Finally, substituting expressions (8)-(10) to Equation (3), we can obtain the vector field of the system on the central manifold.

$(\begin{array}{l}{\stackrel{\dot{}}{x}}_1=[\left(\frac{a^2+d+abd+ab-a}{1+a}-\frac{\left(a+d+1\right)\left(d^3+a{d}^2+d^3\right)}{\left(a+1\right)\left(d^2+ad+d\right)}\right)c_1\\ +\left(\frac{a+d+1}{\left(a+1\right)\left(d^2+ad+d\right)}-\frac{a}{1+a}\right)b_1+\left(\frac{a^2}{1+a}+\frac{a+d+1}{\left(a+1\right)\left(d^2+ad+d\right)}\right)a_1{b}_1]x_1^3\\ +\left[\left(\frac{a^2-a^{2}b}{1+a}\right)a_1+\frac{ab-a}{a+1}\right]x_1+\left(\frac{a^2-a^{2}b}{1+a}\right)a_2{x}_{1}b+\cdots \\ \stackrel{\dot{}}{b}=0\end{array}$ (12)

$\begin{array}{l}{a}_1=-\frac{1}{a^2+a+1};a_2=\frac{{\left(a+1\right)}^2}{{\left(a^2+a+1\right)}^2};b_1=\frac{1+a}{ac-2a+c};\\ c_1=\frac{{\left(1+a\right)}^2}{\left(d^2+ad+d\right)\left(ac-2a+c\right)\left(3a+ad+d\right)}\mathrm{.}\end{array}$

The equilibrium point of this system in this expression is $S_p=\left(\mathrm{0,1}\right)$ . To satisfy the condition that pitchfork bifurcation occurs, letting $F\left(x\mathrm{,}b\right)={\stackrel{\dot{}}{x}}_1$ , we calculate the value at $b=1$ of the following equations:

$F\left(0,0\right)=0,{\frac{\partial F}{\partial x_1}|}_{\left(0,1\right)}=0,{\frac{\partial F}{\partial b}|}_{\left(0,1\right)}=0,{\frac{{\partial }^{2}F}{\partial x_1^2}|}_{\left(0,1\right)}=0;$

${\frac{{\partial }^{2}F}{\partial x_1\partial b}|}_{\left(0,1\right)}=-\frac{a^2\left(a_1+a_2\right)}{1+a}\ne 0,$

$\begin{array}{c}{\frac{{\partial }^{3}F}{\partial x_1^3}|}_{\left(0,1\right)}=\left[\frac{a^2+d+abd}{a+1}-\frac{\left(a+d+1\right)\left(d^3+a{d}^2+d^3\right)}{\left(a+1\right)\left(d^2+ad+d\right)}\right]6c_1\\ +\left(\frac{a+d+1}{\left(a+1\right)\left(d^2+ad+d\right)}-\frac{a}{1+a}\right)6b_1\end{array}$

$\begin{array}{c}+\left(\frac{a^2}{1+a}+\frac{a+d+1}{\left(a+1\right)\left(d^2+ad+d\right)}\right)6a_1{b}_1\\ \ne 0.\end{array}$

According to the pitchfork bifurcation theory [17] , when parameter b crosses the critical value $b=1$ , the pitchfork bifurcation occurs at the equilibrium point S_p.

4.2. Hopf Bifurcation

This section focuses on Hopf bifurcation with a as the bifurcation parameter.

Theorem 5. Assume $b>1,c>0,d<0$ , when the parameter a crosses the critical value −1, the Hopf bifurcation will occur at the point S₀.

Proof 5. The characteristic polynomial of Jacobian matrix of the system at S₀ is

$P_4\left(\lambda \right)=\left(\lambda +c\right)\left(\lambda -d\right)\left({\lambda }^2+\left(a+1\right)\lambda +\left(a-ab\right)\right)=0.$

when $b>1,c>0,d<0$ , the eigenvalues are ${\lambda }_1=-c$ , ${\lambda }_2=d$ , ${\lambda }_{\mathrm{3,4}}=\pm \sqrt{b-1}\text{i}$ . From calculating, we obtain ${\frac{\text{d}}{\text{d}a}Re\left(\lambda \left(a\right)\right)|}_{a=-1}=-\frac{1}{2}\ne 0$ , which satisfy the two conditions for generating Hopf bifurcation. According to Hopf bifurcation theory the system generates Hopf bifurcation at S₀.

Normal form theory [15] [16] [18] and central manifold theorem are used for further analysis. According to the conclusion of Theorem 5, the specific values of the relevant parameters are set to $a=-1,b=2,c=1,d=-2$ to satisfy the conditions for the Hopf bifurcation to occur. The corresponding eigenvalues of these parameters are ${\lambda }_1=-\mathrm{1,}{\lambda }_2=-\mathrm{2,}{\lambda }_3=\text{i}\mathrm{,}{\lambda }_4=-\text{i}$ respectively. Converting the system (1) to the following normal form: $\stackrel{\dot{}}{x}=\mathcal{A}x+F\left(x\right)$ where

$\mathcal{A}=\left[\begin{array}{cccc}1& -1& 0& 1\\ 2& -1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -2\end{array}\right]$

The $F\left(x\right)$ has Taylor expansion near the S₀: $F\left(x\right)=\frac{1}{2}\mathcal{B}\left(x\mathrm{,}y\right)+\frac{1}{6}\mathcal{C}\left(x\mathrm{,}y\mathrm{,}z\right)+\mathcal{O}\left({\Vert X\Vert }^4\right)$ where

$\mathcal{B}\left(x\mathrm{,}y\right)=\left[\begin{array}{c}0\\ -2xz\\ 2xy\\ -2yz\end{array}\right]\mathrm{,}\mathcal{C}\left(x\mathrm{,}y\mathrm{,}z\right)=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right]$

Letting $q\in {ℂ}^n$ is the complex eigenvector which is corresponding to the eigenvalue ${\lambda }_3=\text{i}$ , $p\in {ℂ}^n$ is adjoint eigenvector which meets ${\mathcal{A}}^{\text{T}}p=-\text{i}p$ and the standardized condition $\langle p\mathrm{,}q\rangle =1$ . From calculation we have the value of q and p:

$q={\left(\mathrm{1,1}-\text{i}\mathrm{,0,0}\right)}^{\text{T}}\mathrm{,}p={\left(\frac{1}{2}+\frac{1}{2}\text{i}\mathrm{,}-\frac{1}{2}\text{i}\mathrm{,0,}\frac{1}{10}+\frac{3}{10}\text{i}\right)}^{\text{T}}$ (13)

The critical real eigenspace $T^c$ corresponds to ${\lambda }_3\mathrm{,}{\lambda }_4$ and $T^{su}$ corresponds to ${\lambda }_1\mathrm{,}{\lambda }_2$ are both two-dimensional and is spanned by $\left\{Req\mathrm{,}Imq\right\}$ . For any $x\in {ℝ}^n$ , $x=zq+z\bar{q}+y$ , where $z\in {ℂ}^1$ , $zq+\bar{z}\bar{q}\in T^c$ , $y\in T^{su}$ . The complex variable z is the coordinate on $T_c$ , having

$(\begin{array}{l}z=\langle p\mathrm{,}x\rangle \\ y=x-\langle p\mathrm{,}x\rangle q-\langle \bar{p}\mathrm{,}x\rangle \bar{q}\end{array}$ (14)

Therefore, system (33) has the following form in the coordinate system (36):

$(\begin{array}{l}\stackrel{\dot{}}{z}=\text{i}z+\frac{1}{2}{G}_{20}{z}^2+G_{11}z\bar{z}+\frac{1}{2}{G}_{02}{\bar{z}}^2+\frac{1}{2}{G}_{21}{z}^2\bar{z}+z\langle p\mathrm{,}\mathcal{B}\left(q\mathrm{,}y\right)\rangle +\bar{z}\langle p\mathrm{,}\mathcal{B}\left(\bar{q}\mathrm{,}y\right)\rangle +\cdots \\ \stackrel{\dot{}}{y}=\mathcal{A}y+\frac{1}{2}{H}_{20}{z}^2+H_{11}z\bar{z}+\frac{1}{2}{H}_{02}{\bar{z}}^2+\cdots \end{array}$ (15)

where

$G_{20}=\langle p\mathrm{,}\mathcal{B}\left(q\mathrm{,}q\right)\rangle \mathrm{,}{G}_{11}=\langle p\mathrm{,}\mathcal{B}\left(q\mathrm{,}\bar{q}\right)\rangle \mathrm{,}{G}_{02}=\langle p\mathrm{,}\mathcal{B}\left(\bar{q}\mathrm{,}\bar{q}\right)\rangle$

$(\begin{array}{l}{H}_{20}=\mathcal{B}\left(q\mathrm{,}q\right)-\langle p\mathrm{,}\mathcal{B}\left(q\mathrm{,}q\right)\rangle q-\langle \bar{p}\mathrm{,}\mathcal{B}\left(q\mathrm{,}q\right)\rangle \bar{q}\\ H_{11}=\mathcal{B}\left(q\mathrm{,}\bar{q}\right)-\langle p\mathrm{,}\mathcal{B}\left(q\mathrm{,}\bar{q}\right)\rangle q-\langle \bar{p}\mathrm{,}\mathcal{B}\left(q\mathrm{,}\bar{q}\right)\rangle \bar{q}\end{array}$

The central manifold $W_c$ has an expression:

$y=V\left(z,\bar{z}\right)=\frac{1}{2}{w}_{20}{z}^2+w_{11}z\bar{z}+\frac{1}{2}{w}_{02}{\bar{z}}^2+o{\left|z\right|}^3.$

We also have $\langle p\mathrm{,}{w}_{ij}\rangle =0$ , by combining the quadratic forms with respect to z and $\bar{z}$ in the invariance condition in $W_c$ , we can obtain $\stackrel{\dot{}}{y}=V_z\stackrel{\dot{}}{z}+V_{\bar{z}}\stackrel{\dot{}}{\bar{z}}$ , $w_{ij}\in {ℂ}^4$ satisfy the following system of linear equations

$(\begin{array}{l}\left(2\text{i}{\omega }_0{I}_n-A\right)w_{20}=H_{20}\\ -A{w}_{11}=H_{11}\\ \left(-2\text{i}{\omega }_0{I}_n-A\right)w_{02}=H_{02}\end{array}$ (16)

According to all the above formulas, the following results can be obtained:

$\mathcal{B}\left(q\mathrm{,}q\right)={\left(\mathrm{0,0,2}-2\text{i}\mathrm{,0}\right)}^{\text{T}}\mathrm{,}\mathcal{B}\left(q\mathrm{,}\bar{q}\right)={\left(\mathrm{0,0,2,0}\right)}^{\text{T}}\mathrm{,}\mathcal{B}\left(\bar{q}\mathrm{,}\bar{q}\right)={\left(\mathrm{0,0,2}+2\text{i}\mathrm{,0}\right)}^{\text{T}}$

$G_{20}=0,G_{11}=0,G_{02}=0$

$H_{20}={\left(\mathrm{0,0,2}-2\text{i}\mathrm{,0}\right)}^{\text{T}}\mathrm{,}{H}_{02}={\bar{H}}_{20}={\left(\mathrm{0,0,2}+2\text{i}\mathrm{,0}\right)}^{\text{T}}\mathrm{,}{H}_{11}={\left(\mathrm{0,0,2,0}\right)}^{\text{T}}$

$w_{20}={\left(\mathrm{0,,0,}-\frac{2}{5}-\frac{6}{5}\text{i}\mathrm{,0}\right)}^{\text{T}}\mathrm{,}{w}_{11}={\left(\mathrm{0,0,2,0}\right)}^{\text{T}}\mathrm{,}{w}_{02}={\left(\mathrm{0,0,}-\frac{2}{5}+\frac{6}{5}\text{i}\mathrm{,0}\right)}^{\text{T}}\mathrm{.}$

Then the constraint equation of the original system on $W^c$ is

$\begin{array}{c}\stackrel{\dot{}}{z}=\text{i}z+\frac{1}{2}{G}_{20}{z}^2+G_{11}z\bar{z}+\frac{1}{2}{G}_{02}{\bar{z}}^2\\ +\frac{1}{2}\left(G_{21}+2\langle p\mathrm{,}\mathcal{B}\left(q\mathrm{,}{w}_{11}\right)\rangle +\langle p\mathrm{,}\mathcal{B}\left(\bar{q}\mathrm{,}{w}_{20}\right)\rangle \right)z^2\bar{z}+\cdots \end{array}$

Calculating the first Lyapunov coefficient

$\begin{array}{c}{l}_1\left(0\right)=\frac{1}{2{\omega }_0}Re[\langle p\mathrm{,}\mathcal{C}\left(q\mathrm{,}q\mathrm{,}\bar{q}\right)\rangle -2\langle p\mathrm{,}\mathcal{B}\left(q\mathrm{,}{A}^{-1}\mathcal{B}\left(q\mathrm{,}\bar{q}\right)\right)\rangle \\ +\langle p\mathrm{,}\mathcal{B}\left(\bar{q}\mathrm{,}{\left(2\text{i}{\omega }_0{I}_n-A\right)}^{-1}\mathcal{B}\left(q\mathrm{,}q\right)\right)\rangle ]\\ =\frac{4}{5}>0.\end{array}$

Therefore, the Hopf bifurcation is subcritical.

5. Numerical Simulations

In this section we present some numerical simulations of system (1) to verify the above theoretical analysis. For showing the stability of the equilibrium point S₀ of system (1) under the parameter conditions satisfying Theorem 1, we give the time histories of the solutions by Fourth-order Runge-Kutta methods. The parameters are chosen to be $\left(a\mathrm{,}b\mathrm{,}c\mathrm{,}d\right)=\left(\mathrm{2,}-\mathrm{8,1,}-2\right)$ and $\left(a\mathrm{,}b\mathrm{,}c\mathrm{,}d\right)=\left(\mathrm{2,0.5,2,}-3\right)$ , satisfying the first and second cases of Theorem 1, respectively. According to Figure 2, it can be seen that the equilibrium solution is asymptotically stable when the conditions of Theorem 1 are satisfied.

The bifurcation of the system is analyzed in Section 4. According to the results of the theoretical analysis, choosing $a=-1$ , $b=5$ , $c=1$ , $d=-2$ , at this point, the parameter conditions for the system to generate Hopf bifurcation are satisfied. We give the time histories, Lyapunov exponential spectrum, projections of four-dimensional phase portrait diagram of the model at this time.

Figure 3(a) shows that in this case the system is in periodic motion and also in a steady state. From the Lyapunov exponents spectrum in Figure 3(b), we can see the Lyapunov characteristic index satisfy $\left(\mathrm{0,}-\mathrm{,}-\mathrm{,}-\right)$ which means stable closed orbit will occur. Figure 4 and Figure 5 show the two-dimensional and three-dimensional projections of the system when the Hopf bifurcation condition is satisfied, respectively.

Figure 2. Time histories of system (1) for different parameters. (a) $a=2,b=-8,c=1,d=-2$ ; (b) $a=2,b=0.5,c=2,d=-3$ .

Figure 3. Bifurcation numerical simulations. (a) $a=-1,b=5,c=1,d=-2$ ; (b) $a=-1,b=5,c=1,d=-2$ .

Figure 4. The two-dimensional projection of system (1). (a) The two-dimensional projection of x-y plane; (b) The two-dimensional projection of z-w plane.

Figure 5. The three-dimensional projection of system (1). (a) The three-dimensional projection of x-y-z space; (b) The three-dimensional projection of x-y-w space.

6. Conclusion

The local dynamics of a new four-dimensional Lorenz-like system is investigated with both analytical and numerical methods in this paper. Starting with a discussion of the number of equilibrium points, the stability is analyzed under different parameters by applying the Routh-Hurwitz criterion. The parameter conditions corresponding to the occurrence of pitchfork bifurcation and Hopf bifurcation, respectively, are given, and the related properties of the bifurcations are analyzed. Finally, numerical simulations verify the correctness of the theoretical parts and show the change in the number of equilibrium points as the pitchfork bifurcation parameter crosses the critical value.

Acknowledgements

The project was supported by National Natural Science Foundation of China (11772148, 12172166) and China Postdoctoral Science Foundation (2013T60531).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Perevoznikov, E. and Mikhailova, O. (2015) Neutrality Criteria for Stability Analysis of Dynamical Systems through Lorentz and Rossler Model Problems. Journal of Applied Mathematics and Physics, 3, 569-576. https://doi.org/10.4236/jamp.2015.35070
[2]	Liu, Y.J. and Yang, Q.G. (2009) Dynamics of a New Lorenz-Like Chaotic System. Nonlinear Analysis: Real World Applications, 11, 2563-2572. https://doi.org/10.1016/j.nonrwa.2009.09.001
[3]	Wu, R.C. and Fang, T.B. (2015) Stability and Hopf Bifurcation of a Lorenz-Like System. Applied Mathematics and Computation, 262, 335-343. https://doi.org/10.1016/j.amc.2015.04.072
[4]	Hu, S.G. and Zhou, L.L. (2022) Local Dynamics of a New Four-Dimensional Quadratic Autonomous System. Journal of Applied Mathematics and Physics, 10, 2632-2648. https://doi.org/10.4236/jamp.2022.109177
[5]	Aziz, M.M. and Merie, D.M. (2021) A Four-Dimensional Chaotic System with Hidden Attractor and Its New Proposed Electronic Circuit. Open Access Library Journal, 8, e8101. https://doi.org/10.4236/oalib.1108101
[6]	Fu, H.R., Dong, Y.G. and Zhang, J.G. (2018) Chaotic Secure Communication System of Complex Net Works Based on Three-Dimensional Autonomous System. Journal of Jilin University (Science Edition), 56, 420-425.
[7]	Wang, W.T., Sun, J.Y., Zhang, H. and Zhang, J. (2023) Quantum Cryptosystem and Circuit Design for Color Image Based on Novel 3D Julia-Fractal Chaos System. Quantum Information Processing, 22, Article No. 64. https://doi.org/10.1007/s11128-022-03823-z
[8]	Zhang, X.S., Sun, Z., Zhao, L., Yan, X.S., Zhao, J.J. and Shi, Z.G. (2021) Analysis of Supercritical Pitchfork Bifurcation in Active Magnetic Bearing-Rotor System with Current Saturation. Nonlinear Dynamics, 104, 103-123. https://doi.org/10.1007/s11071-021-06220-w
[9]	Wang, H., Wang, S., Gu, Y. and Yu, Y. (2023) Hopf Bifurcation Analysis of a Two-Dimensional Simplified Hodgkin-Huxley Model. Mathematics, 11, Article No. 717. https://doi.org/10.3390/math11030717
[10]	Wen, S.T., Wang, X.L. and Zhang, G.H. (2023) Hopf Bifurcation in a Reaction-Diffusion-Advection Model with Two Nonlocal Delayed Density-Dependent Feedback Terms. Communications in Nonlinear Science and Numerical Simulation, 119, Article ID: 107080. https://doi.org/10.1016/j.cnsns.2022.107080
[11]	Li, Y.Y., Xiao, L. and Wei, Z.C. (2020) Zero-Hopf Bifurcation Analysis in an Inertial Two-Neural System with Delayed Crespi Function. The European Physical Journal Special Topics, 229, 953-962. https://doi.org/10.1140/epjst/e2020-900159-8
[12]	Wang, X.Y. and Wang, M.J. (2008) A Hyperchaos Generated from Lorenz System. Physica A: Statistical Mechanics and Its Applications, 387, 3751-3758. https://doi.org/10.1016/j.physa.2008.02.020
[13]	Xie, G.B. and Gao, Z.X. (2019) Plaintext Related Color Image Encryption Algorithm Based on Multiple Chaotic Maps. Computer Engineering and Design, 40, 920-930.
[14]	Liao, X.X. (2019) Theory Methods and Application of Stability. 2nd Edition, Huazhong University of Science and Technology Press, Wuhan, 14.
[15]	Lu, Q.S., Peng, L.P. and Yang, Z.Q. (2010) Ordinary Differential Equations and Dynamical Systems. Beijing University of Aeronautics and Astronautics Press, Beijing.
[16]	Wu, W.J., Chen, Z.Q. and Yuan, Z.Z. (2008) Local Bifurcation Analysis of a Four-Dimensional Hyperchaotic System. Chinese Physics B, 17, 2420. https://doi.org/10.1088/1674-1056/17/7/015
[17]	Wiggins, S. (2003) Introduction to Applied Nonlinear Dynamical System and Chaos. Springer, New York.
[18]	Kuznetsov, Y.A. (1998) Elements of Applied Bifurcation Theory. Springer, New York.