Heteroclinic Loop and Homoclinic Loop in a Controlled Chen System

Suqi Ma^1*, Meisheng Li², Bohan Ma², Jie Sheng¹

¹Department of Mathematics, China Agricultural University, Beijing, China.
²Department of Mathematics, Beihang University, Beijing, China.

DOI: 10.4236/ijmnta.2025.144004                

Abstract

The simulation of the Heteroclinic Loop and Homoclinic Loop in a controlled Chen system is finished. The controlled Chen system is $Z_2$ symmetric, and the limit cycle loop or the attractor is observed as a varying free parameter. As the loop becomes the boundary of the unstable manifold of the equilibrium solution, the heteroclinic orbit from the unstable equilibrium solution to the loop is formed. Usually, the twins’ unstable manifold appears due to $Z_2$ symmetry. The Generalized Hopf point brings forth the limit point cycle bifurcation, and nearby, the homoclinic bifurcation is observed. The homoclinic bifurcation arises since the equilibrium solution undergoes a Bogdanov-Takens bifurcation of codimension two. A novel phenomena of multi-loop coexistence are observed near the intersection point of the homoclinic bifurcation curve and Hopf line. Later, homoclinic curve tangents to the Bautin bifurcation line, the stable limit cycle expands into a homoclinic solution, which is called a limit point homoclinic solution.

Keywords

Homoclinic Solution, Chen System, Bogdanov-Takens Bifurcation

1. Introduction

In 1963, Lorenz observed the first chaos attractor and developed the chaos system named the Lorenz system [1]. Even from then on, people devote enormous work to finding the chaos, and the research work has more focus on the new singularity of system complex dynamics to understand it [2]-[5]. In 1999, Chen developed the new chaos system, which is the Chen system, to have new symmetry. Chen system is similar to the Lorenz system; however, they aren’t topologically equivalent [6] [7]. The development of the Chen system control method has referred to the papers [8] [9] due to its broad application in secure communication, noise control fields, etc.

Chen system has $Z_2$ symmetry, its unstable manifold, which is two-dimensional, is mirror symmetrical. In general, we have developed the controlled Chen system as follows:

$\begin{array}{l}{x}^{\prime }=a\left(y-x\right)\hfill \\ y^{\prime }=\left(c-a\right)x-xz+cy+e_1{x}^2+e_{2}yz\hfill \\ z^{\prime }=xy-bz\hfill \end{array}$ (1)

With $e_1=0$ , the symmetry unstable manifold is separated by two twin attractors from each other after parameter perturbation, which inspires us further new enthusiasm with interests in its heteroclinic orbits and homoclinic orbits simulation by numerical computation. We apply the Matcont software to compute some novel bifurcation points, and in this way, the homoclinic orbits are computed. The continuation of homoclinic orbits as the free parameters are continuously varied is done, and the homoclinic bifurcation line is tracked.

With $e_1=0$ , we assume system (1) has zero $O$ and the other two equilibrium points $P$ and $Q$ . Suppose $J={\left(-1,-1,1\right)}^{\text{T}}$ , it is seen that $F\left(J\left(\begin{array}{l}x\hfill \\ y\hfill \\ z\hfill \end{array}\right)\right)=JF\left(\begin{array}{l}x\hfill \\ y\hfill \\ z\hfill \end{array}\right)$ , hence after

$JP=Q,JQ=P$ (2)

And we have the following proposition,

If the equilibrium $P,Q$ is unstable, that is, the characteristic equation has eigenvalues

$\Re \left({\lambda }_1\right)<0<\Re \left({\lambda }_{2,3}\right)$ (3)

(A1) Suppose there exists a heteroclinic orbit $he{t}_1$ from $P$ to limit cycle loop Γ, that is ${\text{limit}}_{t\to -\infty }he{t}_1=P$ , ${\text{limit}}_{t\to +\infty }he{t}_1=\Gamma$ , the new heteroclinic orbits $he{t}_2$ from $Q$ to loop Γ is proved to satisfy the condition ${\text{limit}}_{t\to -\infty }he{t}_2=Q$ , ${\text{limit}}_{t\to +\infty }he{t}_2=\Gamma$ .

Proof: It is easily seen that $he{t}_2=Jhe{t}_1$ ,$he{t}_1=Jhe{t}_2$ . In addition, the two unstable manifolds originated from $P$ and $Q$ is coincide with each other. We also easily get that ${\text{limit}}_{t\to +\infty }het=\Gamma$ since ${\text{limit}}_{t\to -\infty }het=O$ .

By Poincaré section computation, the bifurcation diagram of the periodicity of limit cycles is drawn versus the varying free parameter. The routes to chaos are verified by the Poincaré section scheme. With $e_1>0$ , the symmetry boundary Γ destroys itself and separates into two loops ${\Gamma }_1\left(x,y,z,e_1\right)$ and ${\Gamma }_2\left(-x,-y,z,-e_1\right)$ , hence we derive the conclusion.

(A2) For $e_1>0$ , Suppose there exists a limit cycle loop ${\Gamma }_1\left(x,y,z,e_1\right)$ , with heteroclinic $he{t}_1$ given that ${\lim }_{t\to +\infty }he{t}_1={\Gamma }_1$ and ${\lim }_{t\to -\infty }he{t}_1=P$ , then the new limit cycle loops ${\Gamma }_2\left(-x,-y,z,-e_1\right)$ appears, which lie on the boundary of the unstable manifold of the equilibrium $Q$ and the heteroclinic loop $he{t}_2$ satisfying

${\lim }_{t\to +\infty }he{t}_2={\Gamma }_2,{\lim }_{t\to -\infty }he{t}_2=Q$ (4)

Proof: It is easy to get that $Q\left(-e_1\right)=JP\left(e_1\right)$ , and suppose that $he{t}_2=Jhe{t}_1$ , the conclusion is given.

For $e_1=0$ , the heteroclinic orbit is symmetric, however, the unstable manifold of $P,Q$ are coincide with saddle $O$ , that is, the same cycle loop lying at the boundary of its unstable manifold. Sometimes, the limit cycle loop is evolving into two homoclinic cycles across the zero point, that is, ${\lim }_{t\to \pm \infty }\Gamma =O$ . We often observe the homoclinic bifurcation phenomena as varying two free parameters by continuing the limit cycle loop. For example, if the equilibrium happens with Bogdanov-Takens (BT) bifurcation, which is of singularity 2 [10]-[12], the homoclinic bifurcation may occur nearby. The generalized Hopf bifurcation happens with the first Lyapunov exponent of normal from equal to zero. Since the hysteresis phenomena near the generalized Hopf (GH) point [10], there are two coexisting limit cycles nearby. BT bifurcation brings forth the homoclinic bifurcation phenomena; therefore, if a homoclinic loop exists, which becomes the boundary of a stable limit cycle? In fact, a novel phenomenon of two loops with one homoclinic orbit coexisting occurs at parameters near the intersection of the homoclinic line and Hopf curve. With the perturbation of $e_1$ , the continuation of the homoclinic bifurcation is completed further. Later on, the homoclinic line emerges into a limit point cycle bifurcation line, which expands the related stable limit cycle into a limit point homoclinic solution.

The whole paper is organized into three sections. In section 2, the heteroclinic orbit links the limit cycle and the two unstable equilibrium solutions, respectively, through the unstable manifold. By observing the periodic windows of the Poincaré bifurcation diagram, two separated limit cycles that lie in the boundary of the twin unstable manifold with $\pm e_1$ are simulated. In section 3, the generalized Hopf bifurcation is discussed, and its near dynamics are classified geometrically. In section 4, the homoclinic bifurcation is continued by using Matcont software [13]-[15]. The continuation of a pair of symmetrical homoclinic orbits is tracked along the homoclinic line bifurcating from the BT point. The discussion is given finally.

2. The Heteroclinic Loop Simulation

In this section, we suppose the equilibrium solution $P$ and $Q$ are unstable saddle-focus. The general algorithm for solving the ODE equation is carried out to compute the heteroclinic solution, which is plotted in Figure 1. It is directly drawn with $e_1=0$ , $e_2=2.2$ and

$W^{u}O=he{t}_1=W^{s}P,W^u\left(O\right)=he{t}_2=W^s\left(Q\right),JP=Q,he{t}_2=Jhe{t}_1$

in Figure 1(a). The limit cycle Γ is observed with $e_1=0,e_2=22$ and

$W^{u}P=he{t}_1,W^u\left(Q\right)=he{t}_2,JP=Q,he{t}_2=Jhe{t}_1$

which also satisfy

${\lim }_{t\to +\infty }he{t}_1=\Gamma ,{\lim }_{t\to +\infty }he{t}_2=\Gamma ,{\lim }_{t\to -\infty }he{t}_1=O,{\lim }_{t\to -\infty }he{t}_2=O$

This verifies that the unstable manifold of $P,Q$ and $O$ are in coincidence. Further, we set $e_1=0$ , $e_2=0.6$ , then the chaos is observed, which is the unstable manifold of $P$ and $Q$ , as shown in Figure 2. Other parameters are fixed as $a=35$ , $c=28$ , $b=3$ .

As the parameter $e_1$ varies, the symmetry limit cycle, which lies at the boundary of the unstable manifold, is destroyed. As shown in Figures 3(a)-(f) and Figures 4(a)-(c). The periodicity of the limit cycles is analyzed by Poincaré section computation,

Figure 1. The heteroclinic orbit is observed. (a) As $e_2=22$ , the heteroclinic orbit starts from the zero solution to the nontrivial equilibrium; (b) As $e_2=2.2$ , the heteroclinic orbit starts from an unstable equilibrium to a limit cycle.

Figure 2. The chaos is simulated with $e_2=0.6$ . (a) Sight in $X-Y-Z$ view; (b) View in $X-Y$ plane.

Figure 3. The perturbation of the parameter $e_1$ brings forth the periodical solutions ${\text{Γ}}_{1,2}$ , which have a relationship ${\Gamma }_2\left(-x,-y,z,-e_1\right)=J{\Gamma }_1\left(x,y,z,e_1\right)$ . Fixed parameters are chosen as $a=10$ , $b=5$ , $c=8$ , $e_2=0.6$ . (a) With $e_1=\pm 0.12$ ; (b) With $e_1=\pm 0.14$ ; (c) With $e_1=\pm 0.16$ ; The period doubling bifurcation happens with varying parameter $e_1$ further; (d) With $e_1=\pm 0.29$ ; (e) With $e_1=\pm 0.32$ ; (f) With $e_1=\pm 0.336$ . The observed attractors are plotted with $a=15$ , $b=0.6$ , $c=12$ , $e_2=0.6$ ; (g) With $e_1=0.28$ ; (h) With $e_1=0.3$ ; (i) With $e_1=0.316$ .

Figure 4. The near dynamics of the generalized Hopf point. The bifurcation diagram has three different regimes, which represent the topological and geometrical classification of the near dynamical behavior.

which is often the route to chaos via period-doubling bifurcation of periodic solutions. Via the routes to chaos, the periodical window of solutions is observed as a varying parameter $e_1$ . For example, if we choose parameters with $a=10$ , $b=5$ , $c=8$ , $e_2=0.6$ , we draw Figure 3 with varying parameter $e_1$ . The periodic solutions of its return map with different periodicities are observed. Choose $e_1=0.12,0.14,0.16$ and $e_1=0.29,0.32,0.336$ , the periodical oscillations are observed with asymmetrical periodical solutions, which have periodicity with routes to quasi-periodical solutions shown in Figures 2(a)-(c), and the period-doubling bifurcation shown in Figures 2(d)-(f).

With $a=15$ , $b=0.6$ , $c=12$ , $e_2=0.6$ , the similar numerical simulation is carried out to compute the periodical window as varying $e_1$ , the asymmetry solutions are drawn with opposite value $\pm e_1$ . With the above conclusion, we can express the periodical solutions by ${\Gamma }_2\left(-x,-y,z,-e_1\right)=J{\Gamma }_1\left(x,y,z,e_1\right)$ . As shown in Figures 3(a)-(c), the attractors with multi-period cycles are simulated.

3. Bifurcation of GH Analysis

The generalized Hopf bifurcation (GH) is a Hopf bifurcation that occurs when the first Lyapunov exponent equals zero. For system (1), as the imaginary roots cross the imaginary axis, Hopf bifurcation and the stability property are lost in its instability state. In general, Sub-critical Hopf or super-critical Hopf bifurcation occurs, which bifurcates periodic solutions. The stability of the bifurcating periodical solutions is determined by normal form coefficients. In general, the periodic solution is stable if the first Lyapunov exponent is negative, and is unstable otherwise. Therefore, Bautin bifurcation occurs when the first Lyapunov exponent is zero, which is also called the GH bifurcation [13] [14]. Usually, two limit cycles collide at the GH bifurcation line, and the near dynamics of the GH point are analyzed by the classifying method. Herein, with the assumption $e_1=0$ , we do system stability analysis by the linearizing method.

Suppose system (1) has an equilibrium solution $P\left(x_{\text{*}},y_{\text{*}},z_{\text{*}}\right)$ and $Q\left(-x_{\text{*}},-y_{\text{*}},z_{\text{*}}\right)$ . Do axis transformation and the linearized system at $P$ is written as

$\begin{array}{l}{x}^{\prime }=a\left(y-z\right)\hfill \\ y^{\prime }=\left(-a+c-z_{\text{*}}\right)x+\left(e_2{z}_{\text{*}}+c\right)y+\left(e_2{y}_{\text{*}}-x_{\text{*}}\right)z\hfill \\ z^{\prime }=y_{\text{*}}x+x_{\text{*}}y-bz\hfill \end{array}$ (5)

Therefore, the characteristic equation is written as

$\begin{array}{c}\text{Δ}\left(\lambda \right)=-{\lambda }^3+\left(e_2{z}_{\text{*}}-a-b+c\right){\lambda }^2+\left(a{e}_2{z}_{\text{*}}+b{e}_2{z}_{\text{*}}+e_2{x}_{\text{*}}{y}_{\text{*}}\right)\lambda \\ +\left(-a^2-ab+2ac-a{z}_{\text{*}}+bc-x_{\text{*}}^2\right)\lambda \\ +\left(ab{e}_2{z}_{\text{*}}+a{e}_2{x}_{\text{*}}{y}_{\text{*}}+a{e}_2{y}_{\text{*}}^2-a^{2}b+2abc-ab{z}_{\text{*}}-a{x}_{\text{*}}^2-a{x}_{\text{*}}{y}_{\text{*}}\right)\\ =0\end{array}$ (6)

Set $\lambda =i\omega$ , one gets that

$\begin{array}{l}{\omega }^3-\left(a{e}_2{z}_{\text{*}}+b{e}_2{z}_{\text{*}}+e_2{x}_{\text{*}}{y}_{\text{*}}-a^2-ab+2ac-a{z}_{\text{*}}+bc-x_{\text{*}}^2\right)\omega =0\\ \left(e_2{z}_{\text{*}}-a-b+c\right){\omega }^2=\left(ab{e}_2{z}_{\text{*}}+a{e}_2{x}_{\text{*}}{y}_{\text{*}}+a{e}_2{y}_{\text{*}}^2-a^{2}b+2abc-ab{z}_{\text{*}}-a{x}_{\text{*}}^2-a{x}_{\text{*}}{y}_{\text{*}}\right)\end{array}$ (7)

The BT bifurcation occurs at the intersection point of the fold line and Hopf line, so we get

$\begin{array}{l}a{e}_2{z}_{\text{*}}+b{e}_2{z}_{\text{*}}+e_2{x}_{\text{*}}{y}_{\text{*}}-a^2-ab+2ac-a{z}_{\text{*}}+bc-x_{\text{*}}^2=0\\ ab{e}_2{z}_{\text{*}}+a{e}_2{x}_{\text{*}}{y}_{\text{*}}+a{e}_2{y}_{\text{*}}^2-a^{2}b+2abc-ab{z}_{\text{*}}-a{x}_{\text{*}}^2-a{x}_{\text{*}}{y}_{\text{*}}=0\end{array}$ (8)

Suppose a Hopf bifurcation occurs at the equilibrium solution $P$ , with $\lambda =\pm i\omega \left(\omega >0\right)$ , and the corresponding right eigenvector $q$ and left eigenvector $p$ are written as

$\begin{array}{l}q=\left(\begin{array}{c}-a\\ -a-i\omega \\ -\frac{y_{\text{*}}\left(i\omega +2a\right)}{b+i\omega }\end{array}\right),\\ p=m\left(\begin{array}{c}\frac{\left(2e_1{y}_{\text{*}}-a+c\right)b^2+\left(-\left(2i{y}_{\text{*}}{e}_1-ia+ic\right)\omega +y_{\text{*}}^2\left(e_2-2\right)\right)b+i\omega y_{\text{*}}^2}{b\left(a-i\omega \right)}\\ b-i\omega \\ e_2{y}_{\text{*}}-x_{\text{*}}\end{array}\right)\end{array}$ (9)

After axis transformation, the multilinear form of system (1) is truncated to its third-order form in the style of a Taylor expansion.

$X^{\prime }=LX+B\left(X,X\right)+C\left(X,X,X\right)+o\left({\Vert X\Vert }^3\right)$ (10)

Herein, we derive that

$Lq=i\omega q,L\overline{q}=-i\omega \overline{q},pL=-i\omega p,\overline{p}L=i\omega \overline{p}$

and $\langle p,q\rangle =1$ , $\langle p,\overline{q}\rangle =0$ . Set $X=qz+\overline{q}\overline{z}+W$ , then Equation (10) is written as

$\begin{array}{l}{z}^{\prime }=i\omega z+\langle p,B\left(q,q\right)\rangle z^2+2\langle p,B\left(q,\overline{q}\right)\rangle z\overline{z}+\langle p,B\left(\overline{q},\overline{q}\right)\rangle {\overline{z}}^2\\ +2B\left(qz,W\right)+2B\left(\overline{q}\overline{z},W\right)\\ W^{\prime }=LW+\left(B\left(q,q\right)-q\langle p,B\left(q,q\right)\rangle -\overline{q}\overline{\langle p,B\left(q,q\right)\rangle }\right)z^2\\ +\left(B\left(q,\overline{q}\right)-2q\langle p,B\left(q,\overline{q}\right)\rangle -2\overline{q}\overline{\langle p,B\left(q,\overline{q}\right)\rangle }\right)z\overline{z}\\ +\left(B\left(\overline{q},\overline{q}\right)-q\langle p,B\left(\overline{q},\overline{q}\right)\rangle -\overline{q}\overline{\langle p,B\left(\overline{q},\overline{q}\right)\rangle }\right){\overline{z}}^2+o\left({\Vert z\Vert }^2\right)\\ =LW+H_{20}{z}^2+H_{11}z\overline{z}+H_{02}{\overline{z}}^2+o\left({\Vert z\Vert }^2\right)\end{array}$ (11)

with

$G_{20}=\langle p,B\left(q,q\right)\rangle ,$

$G_{11}=2\langle p,B\left(q,\overline{q}\right)\rangle$

$G_{02}=\langle p,B\left(\overline{q},\overline{q}\right)\rangle$

$G_{21}=3\langle p,C\left(q,q,\overline{q}\right)\rangle$

$H_{20}=\left(B\left(q,q\right)-q\langle p,B\left(q,q\right)\rangle -\overline{q}\overline{\langle p,B\left(q,q\right)\rangle }\right)$

$H_{11}=\left(B\left(q,\overline{q}\right)-2q\langle p,B\left(q,\overline{q}\right)\rangle -2\overline{q}\overline{\langle p,B\left(q,\overline{q}\right)\rangle }\right)$

$H_{02}=\left(B\left(\overline{q},\overline{q}\right)-q\langle p,B\left(\overline{q},\overline{q}\right)\rangle -\overline{q}\overline{\langle p,B\left(\overline{q},\overline{q}\right)\rangle }\right)$

With the assumption $W=W_{20}{z}^2+W_{11}z\overline{z}+W_{02}{\overline{z}}^2+\cdots$ , differentiate both sides of the equality with respect to time $t$ to get

$\begin{array}{l}L{W}_{20}=2i\omega W_{20}-H_{20}\hfill \\ L{W}_{11}=-H_{11}\hfill \\ L{W}_{02}=-2i\omega W_{02}-H_{02}\hfill \end{array}$ (12)

Do inverse computation to get

$\begin{array}{l}{W}_{20}={\left(2i\omega -L\right)}^{\left(-1\right)}{H}_{20}\hfill \\ W_{11}=-L^{\left(-1\right)}{H}_{11}\hfill \\ W_{02}={\left(-2i\omega -L\right)}^{\left(-1\right)}{H}_{02}\hfill \end{array}$ (13)

By near identity transformation, the normal form is written as

$z^{\prime }=i\omega z+C\left(0\right)z^2\overline{z}+o\left({\Vert z\Vert }^3\right)$

with

$C\left(0\right)=2B\left(q,W_{11}\right)+2B\left(\overline{q},W_{20}\right)+\frac{1}{i\omega }\left(G_{20}{G}_{11}-{\left|G_{11}\right|}^2-\frac{2}{3}{\left|G_{02}\right|}^2\right)$ (14)

With Matcont software, the generalized Hopf bifurcation point is calculated along the Hopf line, which satisfies $C\left(0\right)=0$ , the exact condition for the first Lyapunov exponent. The generalized Hopf bifurcation curve is also calculated on the parameter $\left(a,e_2\right)$ plane as varying two free parameters, as shown in Figure 5(a). As often, people intend to understand the near dynamics of the GH point by the classification method. We also draw it in Figure 5(a), by parameter regimes with numbers (1)-(3) to indicate that the near dynamics of the GH point can be classified into different dynamical behaviors topologically. In regime (1), two symmetry stable limit cycles exist, with the equilibrium solutions being unstable. In regime (2), two symmetry pairs of limit cycles arise. In coexisting limit cycle pairs, the red pair denotes the unstable periodical solutions, whilst the blue pairs represent the stable periodical solutions. It is named the hysteresis phenomenon of the limit cycle, which is induced by sub-critical Hopf bifurcation along the right line of the GH point. In regime (3), the equilibrium solutions are stable. It is seen that a pair of limit point cycles is observed along the GH bifurcation curve with green color. The continuation job is finished with varying free parameters, as shown in Figure 5(b), which illustrates the sub-critical Hopf bifurcation phenomena near the GH point. The listed GH point is calculated with parameters chosen as $a=3.8503588$ , $b=2$ ; $c=4.3414$ ,$e_2=2.2321482$ .

Figure 5. The homoclinic bifurcation phenomena happen on the $\left(a-e_2\right)$ parameter plane. (a) The homoclinic bifurcation line with the novelty to cross over GH point, with fixed parameter $b=3$ , $c=3.478447389$ ; (b) The intersection phenomena of homoclinic bifurcation line and Hopf line, and the tangency of homoclinic line with the GH bifurcation line, with fixed parameter $b=2$ , $c=4.341405792$ ; (c) The bifurcating homoclinic solutions of (a); (d) The bifurcating homoclinic solutions of (b).

4. Homoclinic Bifurcation Phenomena

Homoclinic bifurcation phenomena are induced by the separation distance of the orbit itself. One general method to get a homoclinic orbit is to compute the homoclinic orbit starting from the Bogdanov-Takens (BT) bifurcation point. For the BT bifurcation point, the near dynamics of homoclinic bifurcation phenomena are very interesting both in theory and numerical computation. However, it is the novel job of dynamical bifurcation software, and the simulation work of the homoclinic loops is exploited by steps, from the BT point to periodical solutions and the homoclinic orbit. We are familiar with the dynamic hand tools, such as Matcont software, which runs in Windows environments and has a quick speed.

Due to the symmetry property of system (1), we list the following conclusion.

Suppose system (1) has two stable equilibrium solutions $P$ and $Q$ :

(A1) For $e_1=0$ , since $Q=JP$ , a pair of homoclinic bifurcation orbits, ${\Gamma }_1$ and ${\Gamma }_2$ are obtained near the BT bifurcation point, which satisfy

${\lim }_{t\to \pm \infty }{\Gamma }_1=O,{\Gamma }_2=J{\Gamma }_1.$

(A2) Varying free parameter $e_1$ , the homoclinic bifurcation orbit is calculated starting from $e_1=0$ , then continued. $Z-2$ symmetry is destroyed, however, $Q\left(-x_{\text{*}},-y_{\text{*}},z_{\text{*}},-e_1\right)=JP\left(x_{\text{*}},y_{\text{*}},z_{\text{*}},e_1\right)$ . Suppose BT bifurcation happens at $P\left(x_{\text{*}},y_{\text{*}},z_{\text{*}},e_1\right)$ with $e_1=0$ , and a homoclinic orbit ${\Gamma }_1\left(x,y,z,e_1,e_2\right)$ is obtained near the BT point via free parameters perturbation, then at the BT bifurcation point $Q$ , and the homoclinic orbit ${\Gamma }_2\left(-x,-y,z,-e_1,e_2\right)$ is derived by the same computation method, which relates to ${\Gamma }_1$ by ${\Gamma }_2=J{\Gamma }_1$ .

For $e_1=0$ , we carry out the bifurcation computation work by seeking the BT bifurcation point. By computing the roots ${\lambda }_{1,2}=0$ of the characteristic equation in section 2, the BT point is found by tracking the Hopf bifurcation line or fold line on $\left(a,e_2\right)$ -parameter plane, as shown in Figure 6. With chosen parameters $b=3$ ,$c=3.478447389$ , one BT point is calculated at $a=6.957$ , $e_2=-1$ , and the GH point is found at $a=3.813$ , $e_2=0.9023$ . The simulation job of Matcont manifests that the homoclinic bifurcation line through the GH point, and we emphasize that the homoclinic line coincides with the GH bifurcation curve (wherein two limit cycles collide), which is denoted by the pink line in Figure 6(a). As varying two free parameters on the $\left(a,e_2\right)$ plane, the continuation of the homoclinic bifurcation orbits is plotted in Figure 6(c), which is in coincidence with the occurrence of limit cycles collision phenomena. It is called the SNhom orbit, which is significantly noticed as the homoclinic line crossing the GH point.

For $e_1=0$ , with chosen parameters $b=2$ , $c=4.341405792$ , the bifurcation diagram is also drawn in Figure 6(b). The pink color line is the bifurcation line of homoclinic bifurcation, whilst the green color line denotes bifurcation of the limit point cycles. The lines are drawn using the continuation method in Matcont by varying $a$ and $e_2$ continuously. The continuation job of the homoclinic orbit is done, and the red homoclinic orbits are drawn near the BT point, and the blue homoclinic orbits are varying parameters quickly before $P_2$ point, then beyond

Figure 6. The amplification of the homoclinic bifurcation line in the $\left(a-e_2\right)$ parameter plane with $b=2$ , $c=4.341405792$ .

it. The overlap phenomena of the homoclinic orbits manifest the possible complex multi-loop dynamics in system (1). As shown in Figure 7, the intersection point of the homoclinic bifurcation line and Hopf line is denoted as $P_1$ , and it is noticed that the nontrivial equilibrium herein is a center. However, the homoclinic bifurcation line is tangent to the limit point cycle line again at the point $P_2$ . The multi-loop dynamics phenomena appear between points $P_1$ and $P_2$ . In addition to the homoclinic orbit, which usually lies on the red line, one limit point cycle is evolving into one stable and one unstable limit cycle on this segment. Further, after passing into the tangent point $P_2$ , the stable and unstable limit cycles collide into the limit point cycle, which looks like the homoclinic orbit. The continuation of homoclinic orbits and the overlap phenomena are manifest in Figure 6(d).

Perturbation with the free parameter $e_1$ , the homoclinic bifurcation is also simulated with overlapping phenomena, as shown in Figure 7. With chosen parameters $b=3$ , $c=3.478447389$ , the continuation of the homoclinic orbit is done on the $\left(e_1,e_2\right)$ parameter plane. Hence, no symmetry homoclinic orbit is observed. The bifurcation line in the diagram, as shown in Figure 7(a) and Figure 7(b), manifests

Figure 7. The homoclinic bifurcation phenomena happen on the $\left(a-e_2\right)$ parameter plane with $e_1$ parameter perturbation. (a) The homoclinic bifurcation lines simulated with $a=4.03564$ , $b=3$ , $c=3.47844$ ; (b) The homoclinic bifurcation lines simulated with $a=6.2426188$ , $b=2$ , $c=4.341405792$ ; (c) The continuation of homoclinic solutions of (a), with $e_1>0$ , the red line is the homoclinic solution starting at $e_1=0$ , $e_2=1.25$ ; (d) The continuation of homoclinic solutions of (b), with $e_1>0$ , the red line is the homoclinic solution starting at $e_1=0$ , $e_2=0.8$ ; (e) The homoclinic solution overlaps with the periodic solution of doubly period; (f) The phase portraits of the periodical limit cycles of doubly period near $P\left(-3.007,-2.874\right)$ .

the homoclinic bifurcation about the parameter $e_1$ . Suppose ${\Gamma }_1$ is the homoclinic orbit with $e_1>0$ , then another homoclinic orbit ${\Gamma }_2$ is simulated that satisfies ${\Gamma }_2\left(-x.-y,z,-e_1\right)=J{\Gamma }_1\left(x,y,z,e_1\right)$ . As shown in Figure 7(c) and Figure 7(d), the continuation job of the homoclinic orbit is simulated. The red orbit denotes the homoclinic solution at $e_1=0$ . Noticed in Figure 7(d), a novel overlap phenomenon of a homoclinic loop to the phase of a doubly period solution is observed at $P$ with $e_1=-3.007$ , $e_2=-2.874$ . We suppose that the homoclinic bifurcation line is tangent to the period doubling bifurcation line near the point $P$ . By varying the parameter $e_1$ , four $P-2$ periodical limit cycles with double period are simulated, as shown in Figure 7(f).

5. Discussion

The controlled Chen system is dynamically interesting since it involves the simulation of heteroclinic and homoclinic loops. The controlled Chen system is $Z_2$ symmetric, and we observed the symmetrical attractor as a varying free parameter. The system manifests the heteroclinic orbit from the unstable equilibrium solution to the loop, and the boundary of the unstable manifold of the equilibrium solution was formed. Near the generalized Hopf bifurcation point, the limit point cycle was observed. The near dynamics of the Bautin point were classified topologically and geometrically, and the Bautin bifurcation line was drawn. The homoclinic bifurcation line, which starts from the BT point, was computed by Matcont software. It was seen that the homoclinic line either coincides with the Bautin bifurcation line and crosses the GH point, or is tangent to the Bautin line with homoclinic overlap phenomena observed. The perturbation of the free parameter destroys system symmetry. Hence, with $e_1$ either positive or negative, the new homoclinic bifurcation phenomena were discussed. A novel homoclinic loop near the doubly period solutions was simulated due to overlap phenomena as a continuation of homoclinic loops.

Conflicts of Interest

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

References

[1]	Masoller, C., Schifino, A.C.S. and Romanelli, L. (1995) Characterization of Strange Attractors of Lorenz Model of General Circulation of the Atmosphere. Chaos, Solitons & Fractals, 6, 357-366.[CrossRef]
[2]	Lorenz, E.N. (2004) Deterministic Nonperiodic Flow. In: Hunt, B.R., Li, T.Y., Ken-nedy, J.A. and Nusse, H.E., Eds., The Theory of Chaotic Attractors, Springer, 25-36.[CrossRef]
[3]	Sarathy, R. and Sachdev, P.L. (1994) On the Gluing and Ungluing of Strange Attractors in the Study of the Lorenz System. Physics Letters A, 191, 238-244.[CrossRef]
[4]	Al-Sawalha, M.M. and Noorani, M.S.M. (2009) Application of the Differential Transformation Method for the Solution of the Hyperchaotic Rössler System. Communications in Nonlinear Science and Numerical Simulation, 14, 1509-1514.[CrossRef]
[5]	Mackey, M.C. and Glass, L. (1997) Oscillation and Chaos in Physiological Control Systems. Science, 197, 287-289.[CrossRef] [PubMed]
[6]	Chen, G. and Ueta, T. (1999) Yet Another Chaotic Attractor. International Journal of Bifurcation and Chaos, 9, 1465-1466.[CrossRef]
[7]	Ueta, T. and Chen, G. (2000) Bifurcation Analysis of Chen’s Equation. International Journal of Bifurcation and Chaos, 10, 1917-1931.[CrossRef]
[8]	Chen, G. and Dong, X. (1998) From Chaos to Order. World Scientific.[CrossRef]
[9]	Lü, J.H., Chen, G.R. and Zhang, S.C. (2003) A Unified Chaotic System and Its Research. Journal of University of Chinese Academy of Science, 20, 123-129.
[10]	Kuznetsov, Y.A. (1998) Elements of Applied Bifurcation Theory. 2nd Edition, Springer-Verlag.
[11]	Lü, J., Chen, G., Cheng, D. and Celikovsky, S. (2002) Bridge the Gap Between the Lorenz System and the Chen System. International Journal of Bifurcation and Chaos, 12, 2917-2926.[CrossRef]
[12]	Kuznetsov, Y.A. (2011) Practical Computation of Normal Forms on Center Manifolds at Degenerate Bogdanov-Takens Bifurcations. International Journal of Bifurcation & Chaos, 15, 3535-3546.
[13]	Dhooge, A., Govaerts, W. and Kuznetsov, Y.A. (2003) MATCONT: A MATLAB Package for Numerical Bifurcation Analysis of ODEs. ACM Transactions on Mathematical Software, 29, 141-164.[CrossRef]
[14]	Doedel, E.J., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Y.A., Sandstede, B. and Wang X.J. (2000) AUTO97-AUTO2000: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont), User’s Guide. Concordia University, Montreal. http://indy.cs.concordia.ca
[15]	Govaerts, W.J.F. (2000) Numerical Methods for Bifurcations of Dynamical Equilibria. Society for Industrial and Applied Mathematics. Society for Industrial and Applied Mathematics.[CrossRef]