Feedback Control for a Diffusive and Delayed Brusselator Model: Semi-Analytical Solutions

: This paper describes the stability and Hopf bifurcation analysis of the Brusselator system with delayed feedback control in the single domain of a reaction–diffusion cell. The Galerkin analytical technique is used to present a system equation composed of ordinary differential equations. The condition able to determine the Hopf bifurcation point is found. Full maps of the Hopf bifurcation regions for the interacting chemical species are shown and discussed, indicating that the time delay, feedback control, and diffusion parameters can play a signiﬁcant and important role in the stability dynamics of the two concentration reactants in the system. As a result, these parameters can be changed to destabilize the model. The results show that the Hopf bifurcation points for chemical control increase as the feedback parameters increase, whereas the Hopf bifurcation points decrease when the diffusion parameters increase. Bifurcation diagrams with examples of periodic oscillation and phase-plane maps are provided to conﬁrm all the outcomes calculated in the model. The beneﬁts and accuracy of this work show that there is excellent agreement between the analytical results and numerical simulation scheme for all the ﬁgures and examples that are illustrated.


Introduction
Over the last fifty years, many nonlinear chemical applications with diffusion reactions have been explained and studied analytically and numerically via the use of partial differential equations (PDEs), for instance, the Belousov-Zhabotinsky (BZ) reaction [1], the Brusselator model [2], the pellet model [3], the reversible Selkov [4], and the Gray and Scott cubic autocatalytic systems [5,6]. Furthermore, the nonlinear (and quasilinear) diffusive models are of great interest for a large and increasing number of practical applications in many applied mathematics fields, for instance see [7,8] and [9]. These models discuss various chemical oscillatory phenomena in our daily lives using a continuous-flow stirred-tank reactor (CSTR) [10]. The CSTR drives a perfect result in both theoretical and experimental examinations of oscillatory chemical models [2,6,11,12].
One of the most important chemical reaction models was discovered in 1968 [13] and is called the Brusselator model. The nonlinear diffusive Brusselator system consists of two species of BZ model. This model has been discussed by many researchers. For example, ref. [2] examined the analytical results of the Brusselator model in both 1D and 2D domains. Galerkin's technique was applied to determine the ordinary differential equations (ODEs) drawn from the PDEs in the model. As a consequence, perfect solutions were explored to show the stability of the Hopf bifurcation analysis. Many researchers have referred to references concerning the Brusselator model with reaction-diffusion, see [14][15][16][17][18][19] and the references therein.
The existence of feedback with a delay time term in PDE models is extremely important as it can strongly affect the stability of the model. It can also change the results of the Hopf bifurcation analysis [12,20,21]. The authors of [22] considered delayed feedback control for the Brusselator system with reaction-diffusion. They determined that the delayed feedback control caused the formation of extreme events [22,23], where there was a large parameter value, and that the delayed feedback control caused instability. The authors illustrated numerical simulation examples to prove their theoretical outcomes.
A dynamic analytical solution to a nonlinear PDE system can be obtained in many ways. However, some of the models are difficult to examine. Here, I provide an analytical method that is a reliable technique used widely to drive nonlinear PDEs. The Galerkin analytical technique [24] has been used to provide an analytical result that presents a comparison between the numerical simulation PDE model and the theoretical analysis of the ODE system. This method provides a comparison with the analytical solution to demonstrate the accuracy of analytical results. Some examples that use this method include the Selkov-Schnakenberg reaction-diffusion system [25], the BZ model [1], Nicholson's blowflies system [26], neural network systems [27], logistic equations with delays [12,[28][29][30], the Gray and Scott cubic autocatalytic system [6], viral infection models [11], and the foodlimited equation [31]. In general, all the researchers that have applied this method reported an extremely agreeable comparison between the theoretical ODE results and the numerical scheme of the PDE model.
The Brusselator system with delayed feedback control in 1D cells was studied in this work. Research studies are lacking on the effect of a delayed feedback control term in the diffusive Brusselator model. Furthermore, there is not enough detail explaining the effects of delayed feedback control on the overall stability and Hopf bifurcation analysis. Hence, additional research is needed on the effect of these parameters on the PDE model. Thus, our primary objective was to derive the theoretical ODE results by using the Galerkin method. This work offers an explanation of how to find the Hopf bifurcation points. In addition, to construct full map diagrams of the Hopf bifurcation points and the stability analysis (stable and unstable regions), I use examples of periods in the limit cycle maps and the 2D phase-plane maps.
This paper is arranged as follows. Section 2 details the nonlinear PDE model with delayed feedback control under the Dirichlet boundary condition. In Section 3, Galerkin's method is described, which helps to find a system for the delay ODE equations. Section 4 presents the methodology and theoretical framework to present how the Hopf bifurcation maps devolve. Section 5 develops the stability of the concentration cells in the Hopf bifurcation maps and provides the bifurcation diagrams and limited cycle of stable and unstable results. This is in addition to the 2D phase-plane graphs for both the theoretical results of the ODEs and the numerical simulation of the PDE solutions.

Mathematical Formulations Model
The diffusive Brusselator model with delayed feedback control is a nonlinear system of PDEs as follows: where X and Y refer to the activator and inhibitor, respectively, in the dimensionless concentrations of two reactants (interacting chemical species) at time t; A and B indicate two control parameter concentrations during the reaction process; D denotes the diffusion coefficient of two reactants in the system; H refers to the delayed feedback parameter, which is the strength of the feedback control; and τ describes the time delay in the feedback term. Note that the approximation that uses the delay value is useful in relation to the limit where the strength of the delay has both small and large values. In this term, the chemical feedback concentrations are small enough for their development to be explored via a single term of delay. This system is open, and the result has a symmetrical pattern at the center of the domain x = 0. Furthermore, X s > 0 and Y s > 0 are the positive initial concentrations when the time ranges between (−τ, 0). Note that X s = Y s = 0.5 are used in all numerical simulations exemplified in this work. A Crank-Nicolson finite-difference scheme [30,32] is used to obtain the numerical solutions of the PDE model (1). Furthermore, I consider the Runge-Kutta method [3,33] to show the numerical simulation results of the ODE semi-analytical system (3). The numerical examples used for spatial and temporal discretization are ∆x = 0.01 and ∆t = 10 −3 . Note that these values are appropriate to obtain numerical convergence (this implicit scheme is unconditionally stable).

The Galerkin Technique Methodology
Galerkin's technique [24] applies orthogonality to the base function collections to convert PDE into ODE systems. This technique is important and useful. It considers a temporal-spatial separation and explores a spatial form of profile density [3,6]. Galerkin's technique is an analytical method that employs the orthogonality of the rudimentary roles set to find the ODE system from the PDE model. This method has been widely applied to solve nonlinear PDEs in many models and applications [26,27,29,30]. The following equations are trial functions: where ψ 1 = πx 2 and ψ 2 = 3πx 2 . The trial equation function is built as X i = X 1 + X 2 , and variable Y i = Y 1 + Y 2 denotes the middle of the profile concentrations x = 0. The trial equation expansions (2) meet the boundary conditions in the PDE system in (1). The free variables of the values in the model are then explored by computing the values of the PDEs with delay. Next, the PDEs are weighted by two trial function expansions: cos(ψ 1 ) and cos(ψ 2 ). The resulting system of four ODE equations is presented as follows: Note that the equations in (2) have been abbreviated at the two-term expansion point. This is so that the two-term outcome provides sufficient reliability in terms of accuracy without an extreme swell in expression. It also provides good results compared with the numerical scheme of the PDE system. In addition, the system of the one-term equations can be shown via letting X 2 = Y 2 = 0 in the first two equations in the ODE system (3).

Theoretical Framework for the Existence of Hopf Bifurcation
This section analytically describes the methodologies that explore the existence of Hopf bifurcation points. There is a curve that divides the region stability into stable and unstable areas; this curve consists of Hopf bifurcation points. Hopf bifurcation denotes the periodic limit cycle oscillation occurrence in the neighborhood of the steady state that can move the solutions from stable to unstable outcomes. As a result, the conjugated pair of eigenvalues passes into the imaginary axis. For more details, see [20,21,34]. The points of the Hopf bifurcation are determined by discovering the evolution of the Taylor series built for the steady-state parameters. This can be examined as: Therefore, the expressions X i and Y i as in (4) are replaced in the ODE system (3). Then, the steady-state parameters are linearized. Next, the Jacobian matrix of eigenvalues considers a small system perturbation that demonstrates the typical growth value k. This is achieved by placing k = iω in the single equation of the characteristic. This equation can be separated into two different equations: real (Real) and imaginary (Imag). Therefore, the next condition equation can solve the Hopf bifurcation point:

Hopf Bifurcation Areas
This section provides the Hopf bifurcation maps for the stable and unstable regions presented as semi-analytical maps of ODEs (3), comparing the analytical outcomes to the numerical simulation of the PDE results. Furthermore, I studied the effects of all free parameters in the system (τ, H, D) on the regions of the Hopf bifurcation during the reaction process. At the end of this section, I provide a special numerical example to demonstrate the accuracy and advantage of the analytical method over the numerical simulation PDE results. Note that the Hopf bifurcation point can be determined by solving the condition Equation (3); these points compose a curve that divides the stability region into stable and unstable regions. Figure 1a,b show maps of the two areas of the Hopf bifurcation in τ versus concentrations of input reagents A and B, respectively. The positive values that were used are D = 0.01, H = 0.1, and B = 10 in Figure 1a and A = 2 in Figure 1b. In Figure 1a, there are two areas provided: the upper part is the stable zone and the lower part is the unstable zone. It shows that as delay parameter τ increases, the rate of the Hopf bifurcation points for control A reduce slowly. However, in Figure 1b, two different regions are plotted but oppositely to Figure 1a. In contrast, this figure shows that the Hopf bifurcation points for chemical input control B increase steadily as delay parameter τ increases. Therefore, we find good agreement amongst the analytical outcomes of the two terms with the PDE numerical simulation with no more than a 4% error for all of the selected parameters in the domain up to τ = 50.     [24]. Therefore, the values of feedback parameter H and diffusion term D can significantly impact the control parameter concentrations during the reaction process.   Figure 4a, the two areas considered as the Hopf bifurcation points of the diffusion coefficient D increase linearly as feedback parameter H increases. In addition, Figure 4b shows that the diffusion parameters grow slowly in relation to the small range of the D value versus the increase in delay parameter τ. Therefore, the diffusion coefficient values have more of an impact on feedback control parameter H than delay variable τ.    Comparisons were performed for the special parameters of D = 0.3, H = 0.2, A = 3, and τ = 15. In this case, the points of the Hopf bifurcation exist at B c = 1.455 for the oneterm outcome, whereas B c = 1.469 in the analysis of the two-term result. The numerical simulation result of the PDE system is B c = 1.478. All these results predict the incidence of the Hopf bifurcation analysis shown in this example as reliable and in agreement. The error rate between the numerical simulation and the theoretical results is less than 1%.

Bifurcation Diagrams, Periodic Oscillation, and Phase-Plane Map
In this part, I chose control parameter concentration B as the bifurcation parameter. The steady-state results where τ = 0 (no delay case) are provided. I constructed bifurcation diagrams and studied the impact of the other free parameters in the diagrams and limit cycle maps. Furthermore, the maps of the phase plane in the 2D domain are also plotted. The maps of the bifurcation diagrams were examined at the center of the domain and the limit cycles are plotted over a long time period for time t. Figure 6 displays the steady state of the two chemical concentrations of reaction products X and Y versus control parameter B, where H = A = D = 1 and τ = 0. It reveals the analytical outcomes for the two-and one-term systems and the numerical simulation of the presented PDE scheme. There are unique pattern steady-state results for the concentrations of activator X and inhibitor Y. Figure 6a illustrates that as control parameter B increases, the dimensionless concentration activator X decreases significantly before approaching a minimum point at B = 15. The curve of reactant concentration inhibitor Y increases gradually as control parameter B increases in Figure 6b. When the concentration B value is large, the reactant concentration of activator X has a small value (near zero), while reactant concentration Y increases exponentially. There is very good agreement between the numerical model of the PDE simulations and the analytical results of the ODE system, with no more than a 1% error rate in the domain of the concentration of the reactant.  Figure 7 depicts the bifurcation diagram of the dimensionless chemical concentrations of reactants X and Y versus control parameter concentration B. In this illustration, the PDE numerical simulation and the analytical outcomes of both one and two terms are presented and are shown as being closely related. The positive parameters applied in these figures are A = τ = 1, H = 0.1, and D = 0.01. This example shows the importance of the existence of a time delay in the model with a feedback control term, which can make a positive steady state unstable and thus induce instability. In both cases, the periodic solutions have a supercritical Hopf bifurcation from which the reactant concentrations of X and Y are derived. Thus, the Hopf bifurcation points for the analytical outcomes for both one and two terms are B c = 2.560 and 2.428, respectively. The Hopf bifurcation of the numerical simulation result of the PDE is B c = 2.459. It can be seen that as the value of the chemical control parameter B increases, so does the maximum amplitude of the oscillations. The minimum amplitude is still low, and therefore, the predictions concerning the analytical outcomes of the two terms agree with the numerical simulation of the PDE predictions.   1.25 and Y(0, t) 1.46 with an increase in time after a few periodic oscillations. The agreement between the two-term analytical and numerical simulation PDE results is excellent, with an error rate less than 2% in the steady state over a long time t. Figure 10 shows a limited cycle curve for chemical concentrations X and Y versus time t. The parameters in these figures are the same as in Figure 10, where B = 3. This point forms the unstable region in Figure 7. The periodic solution is shown over a long time. The length of the numerical simulation of the period of the limit cycle for the reactant concentrations is 5.96 and the amplitudes of oscillation are 3.81 and 4.17 for the chemical concentrations of the two reaction products X and Y, respectively. The analytical twoterm results in the limit cycle period and amplitude were approached according to the PDE simulation scheme. The period is 5.89 and the amplitudes are 3.87 and 4.04 for the concentrations of the two reactants of X and Y, respectively, with less than a 1.5% error rate between them.   Figure 11a shows the phase plane of the 2D maps for chemical concentration X against Y (activator and inhibitor), where B = 3, D = 0.01, H = 0.1, τ = 1, and A = 1. The scheme of the numerical PDE simulation and one-and two-term analytical output was determined. The numerical simulation of the PDE outcomes approaches the analytical two-term results. Figure 11b presents the phase plane of the 2D maps showing the analytical results of the two terms with four different examples of chemical control concentration B: B = 3, 3.5, 4, and 4.5. It is evident that as control concentration rate B increases, the limit cycle increases regularly. This result indicates that there is no doubling or chaos in the limit cycle in this system [26].   Figure 12a, when H < 0.15, the solution will be stable; in Figure 12b, the result can move from unstable to stable outcomes after point D 0.7.

Conclusions
In this study, I constructed a lower-order analytical system of delayed feedback control for the Brusselator model in the 1D domain. A system of delay ODEs was derived from Galerkin's technique. Hopf bifurcation analysis was determined both theoretically and numerically. The Hopf bifurcation points consisted of a curve that divided each graph into two different areas of stability (stable and unstable). Hopf bifurcation diagrams for the two concentrations of the respective control parameters were plotted. In addition, the effects of the time delay, feedback control, and diffusion parameters were studied in detail. These parameters have a strong impact on the instability and are significant for the two concentration reactants (interacting chemical species) in the system. The solutions in this work were developed as shown by the various examples of stable and unstable limit cycles and the phase plane in the 2D maps. The comparison of the analytical two-term outcome against the PDE scheme for the numerical simulation results facilitated the use of the semi-analytical technique. It is worth using this method as it is a very effective and precise analytical method in relation to a PDE system. In future work, I plan to apply this method to another chemical reaction cell with two different delays in the feedback control.