# Stability, Chaos, and Bifurcation Analysis of a Discrete Chemical System

Authors: Abdul Qadeer Khan; Ibraheem M. Alsulami; Umbreen Sadiq
Journal: Complexity
 (2022)

---

## Abstract

The local dynamics, chaos, and bifurcations of a discrete Brusselator system are investigated. It is shown that a discrete Brusselator system has an interior fixed pointP1,r if r>0. Then, by linear stability theory, local dynamical characteristics are explored at interior fixed point P1,r. Furthermore, for the discrete Brusselator system, the existence of periodic points is investigated. The existence of bifurcations around an interior fixed point is also investigated and proved that the discrete Brusselator model undergoes hopf and flip bifurcations if r,h∈ℋℬ|P1,r=r,h,h=2−r and r,h∈ℱℬ|P1,r=r,h,h=4/2−r−r2−4r, respectively. The next feedback control method is utilized to stabilize the chaos that exists in the discrete Brusselator system. Finally, obtained results are verified numerically.


---

## Body

## 1. Introduction

### 1.1. Mathematical Formulation of a Continuous-Time Brusselator System

The Brusselator is a theoretical system for a specific form of (cubic) autocatalytic chemical reaction found in a variety of chemical and biological systems. At least one product of an autocatalytic process is also a reactant. Glycolysis, clock reactions (product and reactant concentrations change on a regular basis), ozone formation from atomic oxygen (via triple oxygen atom collisions and enzymatic reactions), and the more esoteric but widely quoted BZ-reaction generate autocatalysis from a mixture of potassium bromate, malonic acid, and manganese sulphate in a solution of heated sulfuric acid. In chemistry and biology, the full Brusselator system is a reaction-diffusion system. Reaction-diffusion models show how the concentration of one or more substances distributed in space (e.g., in a cellular space) changes as a result of local chemical reactions in which the substances are transformed into one another as well as diffusion of the substances spread out over the surface. The reaction components of the Brusselator system consist of a pair of intermediate state variablesX and Y with time varying concentrations x and y reacting with product chemicals A, B, and C whose constraints a, b, and c are kept constants. Therefore, reaction models are established by straightforward application of mass action to the following Brusselator system reactions:(1)A⟶k1X,B+X⟶k2Y+C,2X+Y⟶k33X,X⟶k4D.Moreover, the full Brusselator system is obtained by adding diffusion terms∇2x and ∇2y with μ a parameter proportional to the diffusion ratio of the two species x and y. Therefore, the full continuous-time Brusselator system representing systems of differential equations takes the form as follows:(2)dxdt=a−b+1x+x2y+∇2x,dydt=bx−x2y+μ∇2y.Furthermore, without loss of generality, one considered the rate constantski to 1, and moreover continuous-time Brusselator system (2) takes the following form under the assumption of negligible diffusion, i.e., μ≈0, and ∇2x,μ∇2y [1]:(3)dxdt=ab+1x+x2y,dydt=bx−x2y.Finally, ifa=1, then continuous-time Brusselator system becomes the following form (see page 230 in [2]):(4)dxdt=1−r+1xt+xt2yt,dydt=rxt−xt2yt,with b=r.


### 1.2. Review of Literature and Statement of the Problem

In recent years, many papers have appeared on investigating the dynamical characteristics and numerical solution of continuous-time chemical models. For example, the dynamical characteristics of the following chemical system were studied by Zafar et al. [3]:(5)dxdt=x+α2−α1+xy,dydt=1α3x−α2y−xy,where x is the substrate concentration and y is the intermediate complex. The parameters α1, α2, and α3 are all dimensionless. Edeki et al. [4] used a hybrid approach to investigate the numerical solution of the following biochemical system:(6)α3dxdt=y−α1x−xy,dydt=−y+α1−α2,where xandy represent the substrate concentration at a given time t and α1,α2,α3 represent dimensionless parameters. Brown and Davidson [5] have studied the global bifurcation of a Brusselator system. Guo et al. [6] have investigated the hopf bifurcation of a general Brusselator model with diffusion. Li and Wang [7] have explored the hopf bifurcation and diffusion-driven instability of a Brusselator system. Li and Zhang [8] have studied stochastic bifurcation and stability of a continuous Brusselator system with multiplicative white noise. More precisely, Li and Zhang [8] have transformed the original continuous Brusselator system to an Itô averaging system by polar coordinates transformation and the stochastic averaging method. Furthermore, the local and global stabilities are analyzed by the singular boundary theory and largest Lyapunov exponent, and phenomenological bifurcation is also studies by examining the associated Fokker–Planck equation. Finally, numerical simulations are given to show the effectiveness of the obtained results. Zuo and Wei [9] have studied bifurcations of a diffusive Brusselator system with delayed feedback control. Zhao and Ma [10] have studied global and local bifurcations of a general Brusselator model. Ji and Shen [11] have studied the turing instability of a Brusselator in the reaction-diffusion network. Luo and Guo [12] have explored period-1 evolutions to chaos in a Brusselator. Ma [13] has studied the stochastic hopf bifurcation of a Brusselator system. Fu et al. [14] have investigated bifurcation in a diffusive Brusselator model. On the other hand, due to their efficient computing results and rich dynamical behavior, discrete models governed by difference equations are more appropriate than continuous models [15]. For instance, Kang and Pesin [16] explored the local dynamics of a discrete Brusselator model. Khan [17] recently investigated bifurcations of a discrete chemical model:(7)xn+1=xn+hα1−α2xn−xnyn2,yn+1=yn+hα2xn+xnyn2−yn,which is the discrete version of the following continuous glycolytic oscillator model:(8)dxdt=α1−α2x−xy2,dydt=α2x+xy2−y,by Euler’s forward formula, whereα1,α2 are positive constants, x and y, respectively, represent fructose-6-phosphate and adenosine diphosphate. Naik et al. [18] have explored the bifurcations of the following discrete chemical model:(9)xn+1ch=−xnchynch2−βxnch+α+xnch,yn+1ch=xnchynch2+βxnch,where xnch and ynch are the substrate concentrations at time t, and dimensionless parameters are denoted by α,β. Motivated by the aforementioned studies, the purpose of this paper is to explore the chaos and bifurcations of the following discrete Brusselator system:(10)xt+1=h+1−h1+rxt+hxt2yt,yt+1=yt+rhxt−hxt2yt,which is the discrete form of continuous-time Brusselator system (4) by Euler’s forward formula.


### 1.3. Main Contribution

The purpose of this paper is to investigate local dynamics, chaos, and bifurcations of a discrete Brusselator system (10). More precisely, our main contribution in this paper includes(1)
Local dynamics at a fixed pointP1,r of the discrete Brusselator system (10)(2)
Existence of periodic point of discrete Brusselator system (10)(3)
Existence of possible bifurcations aroundP1,r(4)
Controlling the chaos(5)
Verification of theoretical results numerically


### 1.4. Paper Layout

We study fixed points with the linearized form of the discrete Brusselator system (10) in Section 2. Topological classifications at P1,r are explored in Section 3. In Section 4, periodic points of the Brusselator system are investigated. The bifurcation of a Brusselator system (10) around a fixed point P1,r is explored in Section 5. Section 6 is about the controlling of chaos. The main findings are numerically validated in Section 7. The closing remarks are given in Section 8.


## 1.1. Mathematical Formulation of a Continuous-Time Brusselator System

The Brusselator is a theoretical system for a specific form of (cubic) autocatalytic chemical reaction found in a variety of chemical and biological systems. At least one product of an autocatalytic process is also a reactant. Glycolysis, clock reactions (product and reactant concentrations change on a regular basis), ozone formation from atomic oxygen (via triple oxygen atom collisions and enzymatic reactions), and the more esoteric but widely quoted BZ-reaction generate autocatalysis from a mixture of potassium bromate, malonic acid, and manganese sulphate in a solution of heated sulfuric acid. In chemistry and biology, the full Brusselator system is a reaction-diffusion system. Reaction-diffusion models show how the concentration of one or more substances distributed in space (e.g., in a cellular space) changes as a result of local chemical reactions in which the substances are transformed into one another as well as diffusion of the substances spread out over the surface. The reaction components of the Brusselator system consist of a pair of intermediate state variablesX and Y with time varying concentrations x and y reacting with product chemicals A, B, and C whose constraints a, b, and c are kept constants. Therefore, reaction models are established by straightforward application of mass action to the following Brusselator system reactions:(1)A⟶k1X,B+X⟶k2Y+C,2X+Y⟶k33X,X⟶k4D.Moreover, the full Brusselator system is obtained by adding diffusion terms∇2x and ∇2y with μ a parameter proportional to the diffusion ratio of the two species x and y. Therefore, the full continuous-time Brusselator system representing systems of differential equations takes the form as follows:(2)dxdt=a−b+1x+x2y+∇2x,dydt=bx−x2y+μ∇2y.Furthermore, without loss of generality, one considered the rate constantski to 1, and moreover continuous-time Brusselator system (2) takes the following form under the assumption of negligible diffusion, i.e., μ≈0, and ∇2x,μ∇2y [1]:(3)dxdt=ab+1x+x2y,dydt=bx−x2y.Finally, ifa=1, then continuous-time Brusselator system becomes the following form (see page 230 in [2]):(4)dxdt=1−r+1xt+xt2yt,dydt=rxt−xt2yt,with b=r.


## 1.2. Review of Literature and Statement of the Problem

In recent years, many papers have appeared on investigating the dynamical characteristics and numerical solution of continuous-time chemical models. For example, the dynamical characteristics of the following chemical system were studied by Zafar et al. [3]:(5)dxdt=x+α2−α1+xy,dydt=1α3x−α2y−xy,where x is the substrate concentration and y is the intermediate complex. The parameters α1, α2, and α3 are all dimensionless. Edeki et al. [4] used a hybrid approach to investigate the numerical solution of the following biochemical system:(6)α3dxdt=y−α1x−xy,dydt=−y+α1−α2,where xandy represent the substrate concentration at a given time t and α1,α2,α3 represent dimensionless parameters. Brown and Davidson [5] have studied the global bifurcation of a Brusselator system. Guo et al. [6] have investigated the hopf bifurcation of a general Brusselator model with diffusion. Li and Wang [7] have explored the hopf bifurcation and diffusion-driven instability of a Brusselator system. Li and Zhang [8] have studied stochastic bifurcation and stability of a continuous Brusselator system with multiplicative white noise. More precisely, Li and Zhang [8] have transformed the original continuous Brusselator system to an Itô averaging system by polar coordinates transformation and the stochastic averaging method. Furthermore, the local and global stabilities are analyzed by the singular boundary theory and largest Lyapunov exponent, and phenomenological bifurcation is also studies by examining the associated Fokker–Planck equation. Finally, numerical simulations are given to show the effectiveness of the obtained results. Zuo and Wei [9] have studied bifurcations of a diffusive Brusselator system with delayed feedback control. Zhao and Ma [10] have studied global and local bifurcations of a general Brusselator model. Ji and Shen [11] have studied the turing instability of a Brusselator in the reaction-diffusion network. Luo and Guo [12] have explored period-1 evolutions to chaos in a Brusselator. Ma [13] has studied the stochastic hopf bifurcation of a Brusselator system. Fu et al. [14] have investigated bifurcation in a diffusive Brusselator model. On the other hand, due to their efficient computing results and rich dynamical behavior, discrete models governed by difference equations are more appropriate than continuous models [15]. For instance, Kang and Pesin [16] explored the local dynamics of a discrete Brusselator model. Khan [17] recently investigated bifurcations of a discrete chemical model:(7)xn+1=xn+hα1−α2xn−xnyn2,yn+1=yn+hα2xn+xnyn2−yn,which is the discrete version of the following continuous glycolytic oscillator model:(8)dxdt=α1−α2x−xy2,dydt=α2x+xy2−y,by Euler’s forward formula, whereα1,α2 are positive constants, x and y, respectively, represent fructose-6-phosphate and adenosine diphosphate. Naik et al. [18] have explored the bifurcations of the following discrete chemical model:(9)xn+1ch=−xnchynch2−βxnch+α+xnch,yn+1ch=xnchynch2+βxnch,where xnch and ynch are the substrate concentrations at time t, and dimensionless parameters are denoted by α,β. Motivated by the aforementioned studies, the purpose of this paper is to explore the chaos and bifurcations of the following discrete Brusselator system:(10)xt+1=h+1−h1+rxt+hxt2yt,yt+1=yt+rhxt−hxt2yt,which is the discrete form of continuous-time Brusselator system (4) by Euler’s forward formula.


## 1.3. Main Contribution

The purpose of this paper is to investigate local dynamics, chaos, and bifurcations of a discrete Brusselator system (10). More precisely, our main contribution in this paper includes(1)
Local dynamics at a fixed pointP1,r of the discrete Brusselator system (10)(2)
Existence of periodic point of discrete Brusselator system (10)(3)
Existence of possible bifurcations aroundP1,r(4)
Controlling the chaos(5)
Verification of theoretical results numerically


## 1.4. Paper Layout

We study fixed points with the linearized form of the discrete Brusselator system (10) in Section 2. Topological classifications at P1,r are explored in Section 3. In Section 4, periodic points of the Brusselator system are investigated. The bifurcation of a Brusselator system (10) around a fixed point P1,r is explored in Section 5. Section 6 is about the controlling of chaos. The main findings are numerically validated in Section 7. The closing remarks are given in Section 8.


## 2. Fixed Point with Linearized Form

In this section, we study the existence of equilibria with the linearized form of the discrete Brusselator system (10). We first summarize the result regarding the existence of a fixed point of the discrete Brusselator system (10) in the following Lemma.Lemma 1.
For allr>0, Brusselator system (10) has an interior fixed pointP1,r.Proof.
IfPx,y is an equilibrium of the discrete Brusselator system (10), then(11)x=h+1−h1+rx+hx2y,y=y+rhx−hx2y.
From (11), the computation yields(12)1=xr+x−x2y,y=rx.
Using2nd equation of system (12) into 1st equation, one gets x=1. Now, if x=1, then from 2nd equation of system (12), one gets y=r. Therefore, for all r>0, P1,r is interior fixed point of (10).
Now, one has the map in order to write the linearized form of (10) at Px,y:(13)f1,g1↦xn+1,yn+1,where(14)f1=h+1−h1+rxt+hxt2yt,g1=yt+rhxt−hxt2yt.
Finally, variational matrixJ|Px,y at Px,y of the discrete Brusselator system (10) under the map (13) is as follows:(15)J|Px,y=1−r+1h+2hxtythxt2rh−2hxtyt1−hxt2.


## 3. Topological Classifications atP1,r

This section investigates the topological classifications of system (10) at P1,r. J|P1,r at P1,r is(16)J|P1,r=1+r−1hh−rh1−h.In addition, the corresponding characteristic equation is(17)λ2−p1,rλ+q1,r=0,where(18)p1,r=2−2h+rh,q1,r=1+rh−2h+h2.Finally, roots of (17) are(19)λ1,2=p1,r±Δ2,where(20)Δ=p1,r2−4q1,r=r2h2−4rh2.Now, topological classifications of the Brusselator system aboutP1,r are given based on the sign of Δ<0resp.Δ≥0.Lemma 2.
IfΔ=r2h2−4rh2<0, then the following topological classifications hold aroundP1,rof the Brusselator system(10):(i)
P1,r is stable focus if(21)h<2−rwhere0<r<2;(ii)
P1,r is unstable focus if(22)h>2−r;(iii)
P1,r is nonhyperbolic if(23)h=2−r.Lemma 3.
IfΔ=r2h2−4rh2>0, then the following topological classifications hold aroundP1,rof the Brusselator system(10):(i)
P1,r is stable node if(24)0<h<min42−r−r2−4r,42−r+r2−4r;(ii)
P1,r is unstable node if(25)h>max42−r−r2−4r,42−r+r2−4r;(iii)
P1,r is nonhyperbolic if(26)h=42−r−r2−4r;or(27)h=42−r+r2−4r.


## 4. Periodic Points of the Brusselator System

Theorem 1.
P1,rof the Brusselator system(10) is a periodic point with prime period-1.Proof.
From (10), define(28)Fx,y:=f1,g1,where f1 and f2 are represented in (14). Now, the computation yields(29)F|P1,r=P1,r.
From (29), we can deduce that P1,r of the Brusselator system (10) is a periodic point with prime period-1.Theorem 2.
P1,rof Brusselator system (10) is a periodic point of period-2,3,…,l.Proof.
From (28), one has the required statement(30)F2=h+1−h1+rf1+hf12g1,g1+rhf1−hf12g1⇒F2|P1,r=P1,r,F3=h+1−h1+rf12+hf122g12,g12+rhf12−hf122g12⇒F3|P1,r=P1,r,Fl=h+1−h1+rf1l−1+hf1l−12g1l−1,g1l−1+rhf1l−1−hf1l−12g1l−1⇒Fl|P1,r=P1,r.


## 5. Bifurcation

In the present section, we explore bifurcation analysis of the Brusselator system (10) around a fixed point P1,r.


### 5.1. Hopf Bifurcation atP1,r

If (23) holds, then from (19), one gets λ1,223=1. Hence, from this, one can conclude that at P1,r, the Brusselator system (10) may undergoes hopf bifurcation if respective parameters r,h located in the set as follows:(31)ℋℬ|P1,r:=r,h,h=2−r.However, the following theorem guarantees the fact that for the Bruesslator system (10), there must exist a hopf bifurcation at P1,r.Theorem 3.
AtP1,r, the Bruesslator system (10) undergoes the hopf bifurcation ifr,h∈ℋℬ|P1,r.Proof.
It is noted that ifh is considered as a bifurcation parameter, then the Bruesslator system (10) becomes(32)xt+1=h∗+ε+1−h∗+ε1+rxt+h∗+εxt2yt,yt+1=yt+rh∗+εxt−h∗+εxt2yt,where ε≪1. Further, for ε-dependence the Bruesslator system (32), complex eigenvalues are as follows:(33)λ1,2=pε1,r±ι4qε1,r−pε21,r2,where(34)pε1,r=2−2h∗+ε+rh∗+ε,qε1,r=1+rh∗+ε−2h∗+ε+h∗+ε2.
From (33) and (34), one has(35)λ1,2=2−2h∗+ε+rh∗+ε2±ιh∗+ε24r−r2.
From (35), the computation yields λ1,2=qε1,r and dλ1,2/dε|ε=0=2−r/2≠0. Now, for the occurrence of hopf bifurcation, it is required that λ1,2l≠1,l=1,2,3,4 which is equivalent to p01,r≠−2,0,1,2. However, q01,r=1 if (23) holds, and thus p01,r≠−2,2. Therefore, its only requires that p01,r≠0,1, i.e., r≠2,3. Now, the following transformations are utilized in order to shift the fixed point P1,r into O0,0:(36)ut=xt−1,vt=yt−r.
In view of (32) and (36), one gets(37)ut+1=h∗+ε+1−h∗+εr+1ut+1+h∗+εut+12vt+r−1,vt+1=vt+rh∗+εut+1−h∗+εut+12vt+r.
Ifε=0, normal form of (37) is explored. By Taylor series expansion about O0,0, (37) gives(38)ut+1=1+rh−hut+hvt+rhut2+2hutvt+hut2vt,vt+1=−rhut+1−hvt−rhut2−2hutvt−hut2vt.
Now, the following invertible transformation is easily constructed, which converts the linear part of (38) into canonical form:(39)utvt=h0−rh2−h24r−r2ΦtΨt.
From (38) and (39), one gets(40)Φt+1=2−2h+rh2Φt−h24r−r2Ψt+P¯,Ψt+1=h24r−r2Φt+2−2h+rh2Ψt+Q¯,where(41)P¯=−h24r−r2ΦtΨt−rh32Φt3−h324r−r2Φt2Ψt,Q¯=h2r−2ΦtΨt−rh24r−r2h+1Φt3+h32r−2Φt2Ψt.
From (41), one has(42)∂2P¯∂Φt2|O=0,∂2P¯∂Φt∂ΨtO=−h24r−r2,∂3P¯∂Φt3|O=−3rh3,∂3P¯∂Φt2∂Ψt|O=−h34r−r2,∂2P¯∂Ψt2|O=∂3P¯∂Φt∂Ψt2|O=∂3P¯∂Ψt3|O=0,∂2Q¯∂Φt∂Ψt|O=h2r−2,∂2Q¯∂Φt2|O=0,∂3Q¯∂Φt3|O=−6rh24r−r21+h,∂3Q¯∂Φt2∂Ψt|O=h3r−2,∂2Q¯∂Ψt2|O=∂3Q¯∂Φt∂Ψt2|O=∂3Q¯∂Ψt3|O=0.
Now, the nondegeneracy condition for hopf bifurcation is given by [19–27](43)π=−ℜ1−2λ¯λ¯21−λδ11δ20−12δ112−δ022+ℜλ¯δ21,where(44)δ02=18∂2P¯∂Φt2−∂2P¯∂Ψt2+2∂2Q¯∂Φt∂Ψt+ι∂2Q¯∂ϕt2−∂2Q¯∂Ψt2+2∂2P¯∂Φt∂Ψt|O,δ11=14∂2P¯∂ϕt2+∂2P¯∂Ψt2+ι∂2Q¯∂Φt2+∂2Q¯∂Ψt2|O,δ20=18∂2P¯∂Φt2−∂2P¯∂Ψt2+2∂2Q¯∂Φt∂Ψt+ι∂2Q¯∂Φt2−∂2Q¯∂Ψt2−2∂2P¯∂Φt∂Ψt|O,δ21=116∂3P¯∂Φt3+∂3P¯∂Ψt3+∂3Q¯∂ϕt2∂Ψt+∂3Q¯∂Ψt3+ι∂3Q¯∂Ψt3+∂3Q¯∂Φt∂Ψt2−∂3P¯∂Φt2∂Ψt−∂3P¯∂Ψt3|.
Using values from (42) in (44), one gets(45)δ02=h24r−2−ι4r−r2,δ11=0,δ20=h24r−2+ι4r−r2,δ21=116−2r−2h3+ιh2h4r−r2−6r4r−r21+h.
Finally, utilizing (45) in (43) and after manipulation, one gets(46)π=2−r3324+6r−6r2+r3−2−r44.
Asr,h∈ℋℬ|P1,r, which shows that around P1,r, Bruesslator system (10) undergoes hopf bifurcation. Also, supercritical (resp., subcritical) hopf bifurcation occurs if π<0 (resp. π>0).


### 5.2. Flip Bifurcation aroundP1,r

If (26) holds, then from (19), one gets λ1|26=−1 but λ2|26=−2−r−3r2−4r/2−r−r2−4r≠1or−1, which further implies that P1,r Brusselator system (10) may undergoe flip bifurcation if respective parameters r,h are located in the set:(47)ℱℬ|P1,r=r,h,h=42−r−r2−4r.However, the following result guarantees that aroundP1,r, the Bruesslator system (10) undergo the flip bifurcation.Theorem 4.
AtP1,r, the Bruesslator system (10) undergoes flip bifurcation if r,h∈ℱℬ|P1,r.Proof.
Ifh considered as a bifurcation parameter, then Bruesslator system (10) takes the form that is already presented in (32). Furthermore, Bruesslator system (10) becomes(48)ut+1=1+rh−hut+hvt+rhut2+2hutvt+hut2vt+r−1utε+vtε+rut2ε+2utvtε+ut2vtε,vt+1=−rhut+1−hvt−rhut2−2hutvt−hut2vt−rutε−vtε−rut2ε−2utvtε−ut2vtε,by using (36). Now, by utilizing the transformation,(49)utvt=hh−2+rh−hλ2−1+rh−hXtYt,
equation (48) further takes the form as(50)Xt+1Yt+1=−100λ2XtYt+P^εQ^ε,where(51)P^ε=rhλ2−1+rh−2hh1+λ2ut2+2hλ2−1+rh−2hh1+λ2utvt+hλ2−1+rh−2hh1+λ2ut2vt+r−1λ2−1+rh−2hh1+λ2utε+λ2−1+rh−2hh1+λ2vtε+rλ2−1+rh−2hh1+λ2ut2ε+2λ2−1+rh−2hh1+λ2utvtε+λ2−1+rh−2hh1+λ2ut2vtε,Q^ε=rh2+rh−2hh1+λ2ut2+2h2+rh−2hh1+λ2utvt+h2+rh−2hh1+λ2ut2vt+r−12+rh−2hh1+λ2utε+2+rh−2hh1+λ2vtε+r2+rh−2hh1+λ2ut2ε+22+rh−2hh1+λ2utvtε+2+rh−2hh1+λ2ut2vtε,ut=hXt+htYt,vt=−2+rh−hXt+λ2−1+rh−hYt,ut2=h2Xt2+2XtYt+Yt2,utvt=−h2+rh−hXt2+hλ2−1+rh−h−h2+rh−hXtYt+hλ2−1+rh−hYt2,ut2vt=h2−2+rh−hXt3+λ2−1+rh−h−22+rh−hXt2Yt+2λ2−1+rh−h−2+rh−hXtYt2+λ2−1+rh−hYt3,utε=hXtε+hYtε,vtε=−2+rh−hXtε+λ2−1+rh−hYtε,ut2ε=h2Xt2ε+2XtYtε+Yt2ε,utvtε=−h2+rh−hXt2ε+hλ2−1+rh−h−h2+rh−hXtYtε+hλ2−1+rh−hYt2ε,ut2vtε=h2−2+rh−hXt3ε+λ2−1+rh−h−22+rh−hXt2Ytε+2λ2−1+rh−h−2+rh−hXtYt2ε+λ2−1+rh−hYt3ε.
Now, the center manifold is(52)ZcO=Xt,Yt:Yt=l0ε+l1Xt2+l2Xtε+l3ε3+OXt+ε3.
After calculation, one gets(53)l0=0,l1=h2+rh−h4h−rh−4−rh21−λ22,l2=rh2+rh−h−2+rh−h2−rh2h1−λ22,l3=0.
So, map (50) restricted to ZcO is(54)fxt=−xt+v1xt2+v2xtε+v3xt2ε+v4xtε2+v5xt3+OXt+ε4,where(55)v1=11+λ2rh2λ2−1+rh−h+rh3−2+rh−h2hλ2−1+rh−h+2h2,v2=1h1+λ2r−1λ2−1+rh−h+rhh−λ2−1+rh−h+h2+rh−h,v3=1h1+λ22l2h2rhλ2−1+rh−h−rh2+l2h2hλ2−1+rh−h+2h2×λ2−21+rh−h−1+l1hr−1λ2−1+rh−h+rh+l1λ2−1+rh−hλ2−1+rh−h+1+rh−h+hhrλ2−1+rh−h+rh−2+rh−h2λ2−1+rh−h+2h,v4=l21+λ2hr−1λ2−1+rh−h+rh+λ2−1+rh−h+hλ2−1+rh−h,v5=11+λ22hrhλ2−1+rh−h+rh2+2hλ2−1+rh−h+2h2×λ2−21+rh−h−1l1−h1+1+rh−hhλ2−1+rh−h+h2.
In order to show that the Brusselator system undergoes flip bifurcation, the following discriminantΓ1,Γ2 must be nonzero [19, 20]:(56)Γ1=∂2f∂xt∂ε+12∂f∂ε∂2f∂xt2|O=12r2−4rr2−6r+4−2r2−4r≠0,Γ2=16∂3f∂xt3+12∂2f∂xt22|O=10242−r−r2−4r63r2+3r−2+r−1r2−4r2r2−r−r2−4r2−r+r2−4r+18r+r2−4r2−r−r2−4r2+−r+r2−4r2.
From (56), it is noted that Γ2≠0 if r,h∈ℱℬ|P1,r as P1,r. Therefore, the Bruesslator system (10) undergoes flip bifurcation. Further, if Γ2>0, respectively, Γ2<0, then period-2 points bifurcating from P1,r are stable (respectively, unstable).


## 5.1. Hopf Bifurcation atP1,r

If (23) holds, then from (19), one gets λ1,223=1. Hence, from this, one can conclude that at P1,r, the Brusselator system (10) may undergoes hopf bifurcation if respective parameters r,h located in the set as follows:(31)ℋℬ|P1,r:=r,h,h=2−r.However, the following theorem guarantees the fact that for the Bruesslator system (10), there must exist a hopf bifurcation at P1,r.Theorem 3.
AtP1,r, the Bruesslator system (10) undergoes the hopf bifurcation ifr,h∈ℋℬ|P1,r.Proof.
It is noted that ifh is considered as a bifurcation parameter, then the Bruesslator system (10) becomes(32)xt+1=h∗+ε+1−h∗+ε1+rxt+h∗+εxt2yt,yt+1=yt+rh∗+εxt−h∗+εxt2yt,where ε≪1. Further, for ε-dependence the Bruesslator system (32), complex eigenvalues are as follows:(33)λ1,2=pε1,r±ι4qε1,r−pε21,r2,where(34)pε1,r=2−2h∗+ε+rh∗+ε,qε1,r=1+rh∗+ε−2h∗+ε+h∗+ε2.
From (33) and (34), one has(35)λ1,2=2−2h∗+ε+rh∗+ε2±ιh∗+ε24r−r2.
From (35), the computation yields λ1,2=qε1,r and dλ1,2/dε|ε=0=2−r/2≠0. Now, for the occurrence of hopf bifurcation, it is required that λ1,2l≠1,l=1,2,3,4 which is equivalent to p01,r≠−2,0,1,2. However, q01,r=1 if (23) holds, and thus p01,r≠−2,2. Therefore, its only requires that p01,r≠0,1, i.e., r≠2,3. Now, the following transformations are utilized in order to shift the fixed point P1,r into O0,0:(36)ut=xt−1,vt=yt−r.
In view of (32) and (36), one gets(37)ut+1=h∗+ε+1−h∗+εr+1ut+1+h∗+εut+12vt+r−1,vt+1=vt+rh∗+εut+1−h∗+εut+12vt+r.
Ifε=0, normal form of (37) is explored. By Taylor series expansion about O0,0, (37) gives(38)ut+1=1+rh−hut+hvt+rhut2+2hutvt+hut2vt,vt+1=−rhut+1−hvt−rhut2−2hutvt−hut2vt.
Now, the following invertible transformation is easily constructed, which converts the linear part of (38) into canonical form:(39)utvt=h0−rh2−h24r−r2ΦtΨt.
From (38) and (39), one gets(40)Φt+1=2−2h+rh2Φt−h24r−r2Ψt+P¯,Ψt+1=h24r−r2Φt+2−2h+rh2Ψt+Q¯,where(41)P¯=−h24r−r2ΦtΨt−rh32Φt3−h324r−r2Φt2Ψt,Q¯=h2r−2ΦtΨt−rh24r−r2h+1Φt3+h32r−2Φt2Ψt.
From (41), one has(42)∂2P¯∂Φt2|O=0,∂2P¯∂Φt∂ΨtO=−h24r−r2,∂3P¯∂Φt3|O=−3rh3,∂3P¯∂Φt2∂Ψt|O=−h34r−r2,∂2P¯∂Ψt2|O=∂3P¯∂Φt∂Ψt2|O=∂3P¯∂Ψt3|O=0,∂2Q¯∂Φt∂Ψt|O=h2r−2,∂2Q¯∂Φt2|O=0,∂3Q¯∂Φt3|O=−6rh24r−r21+h,∂3Q¯∂Φt2∂Ψt|O=h3r−2,∂2Q¯∂Ψt2|O=∂3Q¯∂Φt∂Ψt2|O=∂3Q¯∂Ψt3|O=0.
Now, the nondegeneracy condition for hopf bifurcation is given by [19–27](43)π=−ℜ1−2λ¯λ¯21−λδ11δ20−12δ112−δ022+ℜλ¯δ21,where(44)δ02=18∂2P¯∂Φt2−∂2P¯∂Ψt2+2∂2Q¯∂Φt∂Ψt+ι∂2Q¯∂ϕt2−∂2Q¯∂Ψt2+2∂2P¯∂Φt∂Ψt|O,δ11=14∂2P¯∂ϕt2+∂2P¯∂Ψt2+ι∂2Q¯∂Φt2+∂2Q¯∂Ψt2|O,δ20=18∂2P¯∂Φt2−∂2P¯∂Ψt2+2∂2Q¯∂Φt∂Ψt+ι∂2Q¯∂Φt2−∂2Q¯∂Ψt2−2∂2P¯∂Φt∂Ψt|O,δ21=116∂3P¯∂Φt3+∂3P¯∂Ψt3+∂3Q¯∂ϕt2∂Ψt+∂3Q¯∂Ψt3+ι∂3Q¯∂Ψt3+∂3Q¯∂Φt∂Ψt2−∂3P¯∂Φt2∂Ψt−∂3P¯∂Ψt3|.
Using values from (42) in (44), one gets(45)δ02=h24r−2−ι4r−r2,δ11=0,δ20=h24r−2+ι4r−r2,δ21=116−2r−2h3+ιh2h4r−r2−6r4r−r21+h.
Finally, utilizing (45) in (43) and after manipulation, one gets(46)π=2−r3324+6r−6r2+r3−2−r44.
Asr,h∈ℋℬ|P1,r, which shows that around P1,r, Bruesslator system (10) undergoes hopf bifurcation. Also, supercritical (resp., subcritical) hopf bifurcation occurs if π<0 (resp. π>0).


## 5.2. Flip Bifurcation aroundP1,r

If (26) holds, then from (19), one gets λ1|26=−1 but λ2|26=−2−r−3r2−4r/2−r−r2−4r≠1or−1, which further implies that P1,r Brusselator system (10) may undergoe flip bifurcation if respective parameters r,h are located in the set:(47)ℱℬ|P1,r=r,h,h=42−r−r2−4r.However, the following result guarantees that aroundP1,r, the Bruesslator system (10) undergo the flip bifurcation.Theorem 4.
AtP1,r, the Bruesslator system (10) undergoes flip bifurcation if r,h∈ℱℬ|P1,r.Proof.
Ifh considered as a bifurcation parameter, then Bruesslator system (10) takes the form that is already presented in (32). Furthermore, Bruesslator system (10) becomes(48)ut+1=1+rh−hut+hvt+rhut2+2hutvt+hut2vt+r−1utε+vtε+rut2ε+2utvtε+ut2vtε,vt+1=−rhut+1−hvt−rhut2−2hutvt−hut2vt−rutε−vtε−rut2ε−2utvtε−ut2vtε,by using (36). Now, by utilizing the transformation,(49)utvt=hh−2+rh−hλ2−1+rh−hXtYt,
equation (48) further takes the form as(50)Xt+1Yt+1=−100λ2XtYt+P^εQ^ε,where(51)P^ε=rhλ2−1+rh−2hh1+λ2ut2+2hλ2−1+rh−2hh1+λ2utvt+hλ2−1+rh−2hh1+λ2ut2vt+r−1λ2−1+rh−2hh1+λ2utε+λ2−1+rh−2hh1+λ2vtε+rλ2−1+rh−2hh1+λ2ut2ε+2λ2−1+rh−2hh1+λ2utvtε+λ2−1+rh−2hh1+λ2ut2vtε,Q^ε=rh2+rh−2hh1+λ2ut2+2h2+rh−2hh1+λ2utvt+h2+rh−2hh1+λ2ut2vt+r−12+rh−2hh1+λ2utε+2+rh−2hh1+λ2vtε+r2+rh−2hh1+λ2ut2ε+22+rh−2hh1+λ2utvtε+2+rh−2hh1+λ2ut2vtε,ut=hXt+htYt,vt=−2+rh−hXt+λ2−1+rh−hYt,ut2=h2Xt2+2XtYt+Yt2,utvt=−h2+rh−hXt2+hλ2−1+rh−h−h2+rh−hXtYt+hλ2−1+rh−hYt2,ut2vt=h2−2+rh−hXt3+λ2−1+rh−h−22+rh−hXt2Yt+2λ2−1+rh−h−2+rh−hXtYt2+λ2−1+rh−hYt3,utε=hXtε+hYtε,vtε=−2+rh−hXtε+λ2−1+rh−hYtε,ut2ε=h2Xt2ε+2XtYtε+Yt2ε,utvtε=−h2+rh−hXt2ε+hλ2−1+rh−h−h2+rh−hXtYtε+hλ2−1+rh−hYt2ε,ut2vtε=h2−2+rh−hXt3ε+λ2−1+rh−h−22+rh−hXt2Ytε+2λ2−1+rh−h−2+rh−hXtYt2ε+λ2−1+rh−hYt3ε.
Now, the center manifold is(52)ZcO=Xt,Yt:Yt=l0ε+l1Xt2+l2Xtε+l3ε3+OXt+ε3.
After calculation, one gets(53)l0=0,l1=h2+rh−h4h−rh−4−rh21−λ22,l2=rh2+rh−h−2+rh−h2−rh2h1−λ22,l3=0.
So, map (50) restricted to ZcO is(54)fxt=−xt+v1xt2+v2xtε+v3xt2ε+v4xtε2+v5xt3+OXt+ε4,where(55)v1=11+λ2rh2λ2−1+rh−h+rh3−2+rh−h2hλ2−1+rh−h+2h2,v2=1h1+λ2r−1λ2−1+rh−h+rhh−λ2−1+rh−h+h2+rh−h,v3=1h1+λ22l2h2rhλ2−1+rh−h−rh2+l2h2hλ2−1+rh−h+2h2×λ2−21+rh−h−1+l1hr−1λ2−1+rh−h+rh+l1λ2−1+rh−hλ2−1+rh−h+1+rh−h+hhrλ2−1+rh−h+rh−2+rh−h2λ2−1+rh−h+2h,v4=l21+λ2hr−1λ2−1+rh−h+rh+λ2−1+rh−h+hλ2−1+rh−h,v5=11+λ22hrhλ2−1+rh−h+rh2+2hλ2−1+rh−h+2h2×λ2−21+rh−h−1l1−h1+1+rh−hhλ2−1+rh−h+h2.
In order to show that the Brusselator system undergoes flip bifurcation, the following discriminantΓ1,Γ2 must be nonzero [19, 20]:(56)Γ1=∂2f∂xt∂ε+12∂f∂ε∂2f∂xt2|O=12r2−4rr2−6r+4−2r2−4r≠0,Γ2=16∂3f∂xt3+12∂2f∂xt22|O=10242−r−r2−4r63r2+3r−2+r−1r2−4r2r2−r−r2−4r2−r+r2−4r+18r+r2−4r2−r−r2−4r2+−r+r2−4r2.
From (56), it is noted that Γ2≠0 if r,h∈ℱℬ|P1,r as P1,r. Therefore, the Bruesslator system (10) undergoes flip bifurcation. Further, if Γ2>0, respectively, Γ2<0, then period-2 points bifurcating from P1,r are stable (respectively, unstable).


## 6. Chaos Control

Several methods can be used for obtaining chaos control in discrete-time models. A few of these methods are the pole-placement technique, the hybrid control method, and the state feedback method [28, 29]. We give a feedback control method to stabilize chaotic orbits at P1,r of the discrete Bruesslator system (10), with the following controlled form:(57)xt+1=h+1−hr+1xt+hxt2yt+ut,yt+1=yt+rhxt−hxt2yt,where ut=−k1xt−x∗−k2yt−y∗ with k1,k2 denote feedback gains. For controlled system (57), JC|P1,r is(58)JC|P1,r=1+hr−1−k1h−k2−rh1−h.Moreover, auxiliary equation ofJC|P1,r is(59)λ2−2−2h+rh−k1λ+h−1k1−rhk2+1−2h+rh+h2=0.Ifλ1,2 are roots of (59), then(60)λ1+λ2=2−2h+rh−k1,and(61)λ1λ2=h−1k1−rhk2+1−2h+rh+h2.The lines of marginal stability are calculated by the solutionλ1=±1 and λ1λ2=1, and these conditions confirm that λ1,2<1. If λ1λ2=1, then from (61), one has(62)L1:h−1k1−rhk2−2h+rh+h2=0.Ifλ1=1, then from (60) and (61), one has(63)L2:k1−rk2+h=0.Finally, ifλ1=−1, then from (60) and (61), one gets(64)L3:2−hk1+rhk2−4−2hr+4h−h2=0.Hence,L1,L2, and L3 in k1,k2-plane define a triangular region that gives λ1,2<1.


## 7. Numerical Simulations

The following three cases are to be considered for the completeness of this section.Case 1.
If one chooses0<r=1.9<2, then from (23), one gets h=0.1, which is the value of the bifurcation parameter. From Lemma 2, interior fixed point P1,1.9 of the discrete Brusselator model (10) is stable (respectively unstable) focus if 0<h<0.1 (respectively, h>0.1). So, if h=0.001<0.1, then Figure 1(a) shows that P1,1.9 of the Brusselator model (10) is stable focus, and moreover Figures 1(b)–1(h) also show the same qualitative behavior for indicated bifurcation value h. On the other hand, if h=0.11>0.1, then Figure 2(a) shows that P1,1.9 of the Brusselator system (10) changes behavior, and as a consequence, stable curves appear. Now, if h=0.11>0.1, then from (46), one gets π=−6.281250000000146×10−6<0, which shows that model (10) undergoes supercritical hopf bifurcation. , Figures 2(b)–2(h) also show that model (10) also undergoes supercritical hopf bifurcation for indicated bifurcation value h. The hopf bifurcation diagrams are drawn in Figure 3. Moreover, the maximum Lyapunov exponent corresponding to Figure 3 is drawn in Figure 4. Finally, 3D bifurcation diagram is plotted in Figure 5.Figure 1
A stable focus of the Brusselator system (10). (a) h = 0.001 with (0.194, 0.195), (b) h = 0.001 with (0.94, 0.95), (c) h = 0.01 with (0.4, 0.5), (d) h = 0.05 with (0.4, 0.5), (e) h = 0.09 with (0.14, 0.15), (f) h = 0.098 with (0.9, 0.9), (g) h = 0.0987 with (0.9, 0.9), and (h) h = 0.0987 with (0.7345, 0.7345).
(a)(b)(c)(d)(e)(f)(g)(h)Figure 2
An unstable focus of the Brusselator system (10). (a) h = 0.11 with (0.3, 0.3), (b) h = 0.13 with (0.53, 0.53), (c) h = 0.2 with (0.93, 0.93), (d) h = 0.28 with (0.93, 0.93), (e) h = 0.3 with (0.63, 0.63), (f) h = 0.38 with (0.63, 0.63), (g) h = 0.4 with (0.63, 0.63), and (h) h = 0.5 with (0.1, 0.1).
(a)(b)(c)(d)(e)(f)(g)(h)Figure 3
Hopf bifurcation diagrams of the Brusselator system (10) with 1.5,1.5 and h∈0.0003,0.55.
(a)(b)Figure 4
Maximum Lyapunov exponents (M. L. E) corresponding to Figure3 of the Brusselator system (10).Figure 5
The 3D bifurcation diagram of the Brusselator system (10).Case 2.
Hereafter, obtained results in subsection5.2 have been verified numerically if r<0. If one choose r=−0.1<0, then from (26), one gets h=2.7403124237432848, which is the value of the bifurcation parameter. From Lemma 3, P1,−0.1 is a stable node if h<min2.7403124237432848,1.459687576256715. So, if h=0.06<min2.7403124237432848,1.459687576256715, then Figure 6(a) shows that P1,−0.1 of the Brusselator system (10) is a stable node and moreover, Figures 6(b)–6(d) also show the same qualitative behavior for the indicated bifurcation value h. On the other hand, if h>2.7403124237432848, then P1,−0.1of the Brusselator system (10) changes behavior, and as a consequence, flip bifurcation occurs, i.e., if h=−0.1, then from (56), one gets Γ1=2.5998052115323693≠0 and Γ2=158.67618818571708>0. Since Γ2=158.67618818571708>0, period-2 points bifurcate from P1,−0.1. So, the flip bifurcation diagrams are drawn in Figure 7. Moreover, the maximum Lyapunov exponent corresponding to Figure 7 is drawn in Figure 8. Finally, the 3D bifurcation diagram is plotted in Figure 9.Figure 6
A stable node of the Brusselator system (10). (a) h = 0.06 with (0.09, 0.09), (b) h = 0.06 with (1.09, 1.09), (c) h = 0.06 with (1.09, 0.09), and (d) h = 0.076 with (0.9, 0.9).
(a)(b)(c)(d)Figure 7
Flip bifurcation diagrams of the Brusselator system (10) with 0.5,0.05 and h∈1.1,2.75.
(a)(b)Figure 8
M. L. E corresponding to Figure7 of the Brusselator system (10).Figure 9
The 3D bifurcation diagram of the Brusselator system (10).Case 3.
Now, numerical simulations are provided to verify obtained results in Section6. For instance, if h=3.5,r=−0.1, then from (62)–(64), one gets(65)L1:2.5k1+0.35000000000000003k2+4.9=0,(66)L2:k1+0.1k2+3.5=0,and(67)L3:−1.5k1−0.35000000000000003k2−1.5499999999999998.
Hence, the lines that are presented in (65), (66), and (67) determine triangular region that gives λ1,2<1 (see Figure 10).Figure 10
Region of stability whereλ1,2<1.Remark 1.
It is noted here that our simulation revels that no flip bifurcation occurs ifr>0.


## 8. Conclusion

This work includes the study of topological classifications around a fixed point, periodic points, chaos, and bifurcations of a discrete Buesslator system (10). We proved that for all parametric values r,h, the Bruesslator system (10) has an interior fixed point P1,r. Furthermore, we studied the local stability with different topological classifications around a fixed points P1,r and proved that it is a stable focus if 0<h<2−r where 0<r<2, unstable focus if h>2−r, nonhyperbolic if h=2−r; stable node if 0<h<min4/2−r−r2−4r,4/2−r−r2−4r, unstable node if h>max4/2−r−r2−4r,4/2−r+r2−4r, and nonhyperbolic if h=4/2−r−r2−4r or h=4/2−r−r2−4r. Next for the Bruesslator system (10), we also studied the existence of periodic points by existing theory. The existence of possible bifurcations around P1,r of the Bruesslator system (10) is explored. We also investigated that around P1,r both hopf and flip bifurcations if r,h∈ℋℬ|P1,r=r,h,h=2−r and r,h∈ℱℬ|P1,r=r,h,h=4/2−r−r2−4r, respectively. Further, the state feedback control method is utilized in order to stabilize the chaos that exists in the discrete Bruesslator system. Finally, the obtained results are numerically verified.



---

*Source: 6921934-2022-09-25.xml*