Vieta–Lucas Polynomials for the Brusselator System with the Rabotnov Fractional-Exponential Kernel Fractional Derivative

Abstract

In this study, we provide an efficient simulation to investigate the behavior of the solution to the Brusselator system (a biodynamic system) with the Rabotnov fractional-exponential (RFE) kernel fractional derivative. A system of fractional differential equations can be used to represent this model. The fractional-order derivative of a polynomial function $t^p$ is approximated in terms of the RFE kernel. In this work, we employ shifted Vieta–Lucas polynomials in the spectral collocation technique. This process transforms the mathematical model into a set of algebraic equations. By assessing the residual error function, we can confirm that the provided approach is accurate and efficient. The outcomes demonstrate the effectiveness and simplicity of the technique for accurately simulating such models.

1. Introduction

For the past 30 years, numerous scientists have been interested in fractional calculus [1]. As a result, scientists were able to offer new definitions of the fractional derivative with non-singular kernels, which were required to satisfy the demand for mathematical modeling of numerous real-life problems in various spheres of our lives, including physics, biology, engineering, viscoelasticity, and fluid mechanics. It was essential to develop and use many of the approximation approaches because it is well known that the majority of fractional differential equations (FDEs) are challenging to be solved exactly [2,3]. The reader is invited to check Reference [4] for more information about the definitions and properties of these fractional derivatives. These days, a wide variety of fractional operators have been developed as generalizations of classical derivatives. The Caputo and Riemann–Liouville fractional derivatives possess a power kernel that extends the classical derivatives. However, we obtain a new generalized class of fractional derivatives if we replace that kernel with exponential or Mittag–Leffler kernels. It is worth pointing out that the derivative with a Mittag–Leffler kernel is known as the Atangana–Baleanu derivative, whereas the derivative with an exponential kernel is known in groundwater flow [5], medical sciences [6], chaos theory [7], and other areas [8].

Spectral methods are among the most significant and efficient tools for simulating differential equations of various types [9]. As an application of these methods, we can see that the spectral Tau method is implemented to solve a general class of FDEs [10], the spectral collocation method (SCM) with shifted third-kind Chebyshev polynomials is used to numerically treatment a type of nonlinear fractional pantograph differential equations [11], and the coefficients of differentiated expansions of double and triple Jacobi polynomials are considered in [12]. The capacity of these techniques to provide us with reliable findings with very small errors is one of their most significant characteristics. For instance, the Vieta–Lucas polynomials’ (VLPs’) orthogonality property is used to approximate periodic functions in some bounded and closed interval $[a,b]$ [13,14]. Spectral methods heavily rely on polynomials. In this work, our goal is to use the RFE kernel to evaluate the fractional derivative. The fractional-order derivative of a polynomial function $t^p$ is approximated in terms of the RFE kernel. We provide the numerical solution to the proposed system using this approximation formula and the characteristics of VLPs. By taking a fractional derivative of polynomials, we demonstrate the validity of the newly derived formula. Additionally, we investigate the Brusselator system using the fractional derivative of the RFE kernel. This allows us to predict the viability of our numerical approach for this model ([15,16]). The SCM and the shifted VLPs are combined to take advantage of the properties of each of them in approximation, on the one hand, and in increasing the applicability and accuracy of the proposed technique, on the other hand, which in turn gives good and close solutions to the real solutions to the problem under study.

It is important to recall that several researchers have recently examined the fractional Brusselator system [17,18]. Gafiychuk and Datsko investigated the stability of this system [18]. The existence of a limit cycle in the solutions of the fractional Brusselator system was demonstrated by Wang and Li using a numerical approach [18]. To approximate the proposed RFE kernel problem, we applied the current approximation technique. We are able to show that this method can be used to solve the given model successfully by using the numerical approach and that there is excellent agreement with the solutions that are currently available in the literature. To control and decrease the relative errors, we can include additional terms from the solution series. The use and potential of the proposed numerical method are demonstrated by comparing the exact and approximate solutions. The comparison with previously published work using a different numerical approach and a different fractional derivative allows us to conclude that the operator without singularity is more suitable for numerical simulations for the model under consideration in this research. It is worth pointing out that the numerical simulations were carried out using the computer software Mathematica.

This study is organized as follows. In Section 2, we consider some preliminaries and notations concerning definitions of the fractional derivatives and shifting Vieta–Lucas polynomials. In Section 3, we present the numerical implementation of the proposed method. The conclusion is given in Section 4.

2. Preliminaries and Notations

2.1. Definitions of Fractional Derivatives

We provide a fresh start by recalling some definitions for fractional derivatives available in the literature, along with some of their properties.

Definition 1.

For the function $p\left(\eta \right)\in H_1(0,b)$, the fractional derivative of order $0<\nu \le 1$ in the Caputo sense is given by

$${}^C{D}^{\nu }p\left(\eta \right)=\frac{1}{\Gamma (1-\nu )}{\int }_0^{\eta }\frac{p^{\prime }\left(\tau \right)}{{(\eta -\tau )}^{\nu }}d\tau ,\eta >0.$$

Definition 2.

For a function $\Theta \left(\eta \right)$, the left-sided Caputo fractional derivative of order κ on the interval $[0,1]$ is defined by

$${}^{RFE}{D}^{\kappa }\Theta \left(\eta \right)={\int }_0^{\eta }{\Theta }^{\left(n\right)}\left(\xi \right){\mathbb{R}}_{\kappa }[-\Omega {(\eta -\xi )}^{\kappa }]d\xi ,n-1<\kappa \le n.$$

(1)

Here, the Rabotnov fractional exponential function for the function $\Omega \in R^{+}$ is defined as follows:

$${\mathbb{R}}_{\kappa }[-\Omega {\left(\eta \right)}^{\kappa }]=\sum _{j=0}^{\infty }\frac{{(-\Omega )}^j{\eta }^{(j+1)(\kappa +1)-1}}{\Gamma \left[\right(j+1\left)\right(\kappa +1\left)\right]}.$$

See [19,20] for additional information regarding the RFE-operator derivative. In the next step, we recall the approximation formula for the fractional derivative concerning the RFE kernel using a readily available numerical integration scheme.

Theorem 1.

The function $g\left(\eta \right)={\eta }^{\theta }$ with $\theta \ge n$ ($n=\lceil \kappa \rceil$) has an RFE derivative of order $n-1<\kappa <n$, which is given by [20]

$$\begin{array}{cc}\hfill {}^{RFE}{D}^{\kappa }{\eta }^{\theta }=& \frac{\Gamma (\theta +1)}{\Gamma (\theta +1-\lceil \kappa \rceil )}\times \frac{h}{3}[G_{\kappa ,\theta }(\eta ,{\gamma }_0)+G_{\kappa ,\theta }(\eta ,{\gamma }_m)\hfill \\ \hfill & +4\sum _{j=1,j-odd}^{m-1}{G}_{\kappa ,\theta }(\eta ,{\gamma }_j)+2\sum _{j=2,j-even}^{m-2}{G}_{\kappa ,\theta }(\eta ,{\gamma }_j)].\hfill \end{array}$$

(2)

Each segment of length h in the domain $[0,1]$ is separated into m equal segments so that, for each $j=0,1,\dots ,m$, we have $h=\frac{1}{m}$, ${\gamma }_j=\frac{j}{m}$, and

$$G_{\kappa ,\theta }(\eta ,\gamma )={\gamma }^{\theta -\lceil \kappa \rceil }{\mathbb{R}}_{\kappa }[-\Omega {(\eta -\gamma )}^{\kappa }].$$

Remark 1.

Since the integral in the formula (1) is difficult to calculate in exact form, we can estimate it by using the Simpson $\frac{1}{3}$ method or any other numerical methodology. It is well known that the formula in (1) is used to generate the formula in (2).

2.2. Shifting Vieta–Lucas Polynomials

In this subsection, we define the shifted Vieta–Lucas polynomials and the properties which are necessary in the present study. We focus our investigation on a class of orthogonal polynomials. It is possible to create a new family of orthogonal polynomials using the recurrence relations and analytical equations of the VLPs as follows.

The VLP, ${\mathrm{VL}}_i\left(t\right)$, of degree $i\in {\mathbb{N}}_0$ is defined by [21]

$${\mathrm{VL}}_i\left(t\right)=2cos\left(i\psi \right),$$

where $\psi =arccos\left(0.5t\right)$, $\psi \in [0,\pi ]$, and $-2\le t\le 2$. It is easy to demonstrate that the following recurrence relation holds for the VLPs:

$${\mathrm{VL}}_i\left(t\right)=t{\mathrm{VL}}_{i-1}\left(t\right)-{\mathrm{VL}}_{i-2}\left(t\right),\forall i=2,3,\dots$$

Here, ${\mathrm{VL}}_0\left(t\right)=2$ and ${\mathrm{VL}}_1\left(t\right)=t$. A new class of orthogonal polynomials on the interval $[0,1]$ is created from VLPs using $t=4\eta -2$. This new class of polynomials is denoted by the symbol ${\mathrm{SVL}}_i\left(\eta \right)$, and the elements are given by the relation

$${\mathrm{SVL}}_i\left(\eta \right)={\mathrm{VL}}_i(4\eta -2).$$

As a consequence, it is easy to check that the following recurrence relation holds:

$${\mathrm{SVL}}_{i+1}\left(\eta \right)=(4\eta -2){\mathrm{SVL}}_{i-1}\left(\eta \right)-{\mathrm{SVL}}_{i-2}\left(\eta \right),\forall i=2,3,\dots ,$$

where ${\mathrm{SVL}}_0\left(\eta \right)=2$ and ${\mathrm{SVL}}_1\left(\eta \right)=4\eta -2.$ Also, observe that ${\mathrm{SVL}}_i\left(0\right)=2{(-1)}^i$ and ${\mathrm{SVL}}_i\left(1\right)=2$ are satisfied for all $i=0,1,2,...$.

The analytical formula for the polynomials ${\mathrm{SVL}}_i\left(\eta \right)$ is given by

$${\mathrm{SVL}}_i\left(\eta \right)=2i\sum _{j=0}^i{(-1)}^j\frac{4^{i-j}\Gamma (2i-j)}{\Gamma (j+1)\Gamma (2i-2j+1)}{\eta }^{i-j},\forall i=2,3,\dots$$

Notice that the polynomials ${\mathrm{SVL}}_i\left(\eta \right)$ are orthogonal polynomials on $[0,1]$ concerning the weight function $\frac{1}{\sqrt{\eta -{\eta }^2}}$. As a consequence, we have that

$$〈{\mathrm{SVL}}_i\left(\eta \right),{\mathrm{SVL}}_j\left(\eta \right)〉={\int }_0^1\frac{{\mathrm{SVL}}_i\left(\eta \right){\mathrm{SVL}}_j\left(\eta \right)}{\sqrt{\eta -{\eta }^2}}d\eta =\left\{\begin{array}{cc}0,\hfill & i\ne j\ne 0,\hfill \\ 4\pi ,\hfill & i=j=0,\hfill \\ 2\pi ,\hfill & i=j\ne 0.\hfill \end{array}\right.$$

Let $\varphi \left(\eta \right)\in L^2[0,1]$. Using the polynomial ${\mathrm{SVL}}_j\left(\eta \right)$, we readily obtain that

$$\varphi \left(\eta \right)=\sum _{j=0}^{\infty }{c}_j{\mathrm{SVL}}_j\left(\eta \right),$$

(3)

where $c_j$ must be evaluated to transform $\varphi \left(\eta \right)$ into the terms of ${\mathrm{SVL}}_i\left(\eta \right)$. By taking into account only the first $m+1$ terms in (3), we define

$${\varphi }_m\left(\eta \right)=\sum _{j=0}^m{c}_j{\mathrm{SVL}}_j\left(\eta \right).$$

(4)

Consequently, it is possible to calculate $c_j$ for each $j=0,1,2,\cdots$ as follows:

$$c_j=\frac{1}{{\delta }_j}{\int }_0^1\frac{{\varphi }_m\left(\eta \right){\mathrm{SVL}}_j\left(\eta \right)}{\sqrt{\eta -{\eta }^2}}d\eta ,{\delta }_j=\left\{\begin{array}{c}4\pi ,j=0,\hfill \\ 2\pi ,j=1,2,\dots ,m.\hfill \end{array}\right.$$

(5)

Lemma 1.

Let us assume that $\varphi \left(\eta \right)\in L_{\tilde{w}}^2[0,1]$, regarding the weight function $\tilde{w}\left(\eta \right)=\frac{1}{\sqrt{\eta -{\eta }^2}}$, with the assumption $|{\varphi }^{″}\left(\eta \right)|\le \epsilon$, for some constant ε. Then, the series (4) uniformly converges to the function $\varphi \left(\eta \right)$ as $m\to \infty$. Additionally, the following estimates are satisfied:

• In Equation (4), the coefficients’ series are bounded, that is,

  $$\left|c_j\right|\le \frac{\epsilon }{4j\left(j^2-1\right)},\forall j>2.$$
• The following inequality applies to the error estimate norm $(L_{\tilde{w}}^2[0,1]$-norm):

  $${∥\varphi \left(\eta \right)-{\varphi }_m\left(\eta \right)∥}_{\tilde{w}}<\frac{\epsilon }{12\sqrt{m^3}}.$$
• The following absolute error bound applies if ${\varphi }^{\left(m\right)}\left(\eta \right)\in C[0,1]$:

  $$∥\varphi \left(\eta \right)-{\varphi }_m\left(\eta \right)∥\le \frac{\Delta {\Pi }^{m+1}}{(m+1)!}\sqrt{\pi }.$$
  Here, $\Delta ={max}_{\eta \in [0,1]}{\varphi }^{(m+1)}\left(\eta \right)$ and $\Pi =max\left\{1-{\eta }_0,{\eta }_0\right\}$.

See the reference [22] for additional information on these polynomials and the convergence analysis for the approximation in (4).

Theorem 2.

The κ-order of the RFE fractional derivative for the function ${\varphi }_i\left(\eta \right)$ which is defined in Equation (4) can be evaluated as follows [20]:

$$\begin{array}{cc}\hfill {}^{RFE}{D}^{\kappa }{\varphi }_i\left(\eta \right)& =\sum _{j=\lceil \kappa \rceil }^i{\chi }_{i,j,\kappa }[G_{\kappa ,p}(\eta ,{\gamma }_0)+G_{\kappa ,p}(\eta ,{\gamma }_m)\hfill \\ \hfill & +4\sum _{\ell =1,\ell -odd}^{m-1}{G}_{\kappa ,p}(\eta ,{\gamma }_{\ell })+2\sum _{\ell =2,\ell -even}^{m-2}{G}_{\kappa ,p}(\eta ,{\gamma }_{\ell })],\hfill \end{array}$$

(6)

where

$${\chi }_{i,j,\kappa }=\frac{h\Gamma (i-j+1)}{3\Gamma (i-j+1-\lceil \kappa \rceil )}\times \frac{{(-1)}^{j}2i{4}^{i-j}\Gamma (2i-j)}{\Gamma (j+1)\Gamma (2i-2j+1)}$$

and

$$G_{\kappa ,p}(\eta ,\gamma )={\gamma }^{p-\lceil \kappa \rceil }{\mathbb{R}}_{\kappa }{[-\Omega {(\eta -\gamma )}^{\kappa }]}_{p=i-j}.$$

Proof. 

By using Theorem 1, we have

$$\begin{array}{cc}\hfill {}^{RFE}{D}^{\kappa }{t}^{i-j}=& \frac{\Gamma (i-j+1)}{\Gamma (i-j-\lceil \kappa \rceil +1)}\times \frac{h}{3}[G_{\kappa ,p}(\eta ,{\gamma }_0)+G_{\kappa ,p}(\eta ,{\gamma }_m)\hfill \\ \hfill & +4\sum _{\ell =1,\ell -odd}^{m-1}{G}_{\kappa ,p}(\eta ,{\gamma }_{\ell })+2\sum _{\ell =2,\ell -even}^{m-2}{G}_{\kappa ,p}(\eta ,{\gamma }_{\ell })].\hfill \end{array}$$

(7)

Each segment in the domain $[0,1]$ has a length of h and is divided into m equal segments. This fact yields that $h=\frac{1}{m}$ and ${\gamma }_{\ell }=\frac{\ell }{m}$ for each $\ell =0,1,2,...,m$, and

$$G_{\kappa ,p}(\eta ,\gamma )={\gamma }^{p-\lceil \kappa \rceil }{\mathbb{R}}_{\kappa }{[-\Omega {(t-\gamma )}^{\kappa }]}_{p=i-j}.$$

Now using Equations (4) and (7), we can evaluate the RFE derivative of ${\varphi }_i\left(\eta \right)$ as follows:

$$\begin{array}{cc}\hfill {}^{RFE}{D}^{\kappa }{\varphi }_i\left(\eta \right)& =\sum _{j=0}^i\frac{{(-1)}^{j}2i{4}^{i-j}\Gamma (2i-j)}{\Gamma (j+1)\Gamma (2i-2j+1)}{}^{RFE}{D}^{\kappa }{\eta }^{i-j}\hfill \\ & =\sum _{j=\lceil \kappa \rceil }^i\frac{\Gamma (i-j+1)}{\Gamma (i-j+1-\lceil \kappa \rceil )}\times \frac{{(-1)}^{j}2i{4}^{i-j}\Gamma (2i-j)}{\Gamma (j+1)\Gamma (2i-2j+1)}\hfill \\ & \times \frac{h}{3}[G_{\kappa ,p}(\eta ,{\gamma }_0)+G_{\kappa ,p}(\eta ,{\gamma }_m)+4\sum _{\ell =1,\ell -odd}^{m-1}{G}_{\kappa ,p}(\eta ,{\gamma }_{\ell })+2\sum _{\ell =2,\ell -even}^{m-2}{G}_{\kappa ,p}(\eta ,{\gamma }_{\ell })].\hfill \end{array}$$

(8)

From this result, we can easily obtain the required Formula (6), as desired. □

3. Numerical Implementation

In this stage of our study, the approach designed in the previous section is used to numerically obtain a semi-analytical solution for the proposed fractional model. In this way, we drastically simplify the problem to a nonlinear system of algebraic equations. To obtain the unknown coefficients of the series solution, this system is converted into a restricted optimization problem. More precisely, we investigate the Brusselator system as formulated in an RFE-fractional form as follows ([17,18]):

$$\begin{array}{cc}\hfill {}^{RFE}{D}^{\varrho }{\psi }_1\left(\eta \right)& =\rho -(\lambda +1){\psi }_1\left(\eta \right)+{\psi }_1^2\left(\eta \right){\psi }_2\left(\eta \right),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {}^{RFE}{D}^{\varrho }{\psi }_2\left(\eta \right)& =\lambda {\psi }_1\left(\eta \right)-{\psi }_1^2\left(\eta \right){\psi }_2\left(\eta \right),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\psi }_1\left(0\right)& ={\psi }_{1,0},{\psi }_2\left(0\right)={\psi }_{2,0}.\hfill \end{array}$$

(11)

In this model, the parameters $\rho >0$, $\lambda >0$, and ${\psi }_{1,0},{\psi }_{2,0}$ are real constants [23]. These equations describe how the concentrations of the two species change over time due to the reaction and diffusion processes. The behavior of the Brusselator model (BM) depends critically on the values of the parameters $\rho$ and $\lambda$. These parameters control the strength of the reaction and diffusion processes and, hence, determine the system’s overall behavior. In particular, the parameter $\rho$ controls the distance from the Hopf bifurcation, a critical point at which the system transitions from a stable to an unstable steady state. The system can exhibit sustained oscillations as $\rho$ increases beyond this critical value. Parameter $\lambda$ controls the balance between the reaction and diffusion processes. The reaction dominates at low values of $\lambda$, and the system can exhibit sustained oscillations or spiral waves. As $\lambda$ increases, diffusion becomes more important, and the system can exhibit Turing patterns.

As we pointed out, we use the approach proposed in this work to numerically solve the system (9)–(11). In the first stage, we approximate ${\psi }_1\left(\eta \right)$ and ${\psi }_2\left(\eta \right)$ by ${\psi }_1^N\left(\eta \right),$ and ${\psi }_2^N\left(\eta \right)$, respectively, in the following formulas:

$${\psi }_1^N\left(\eta \right)=\sum _{i=0}^N{\alpha }_i{\mathrm{SVL}}_i\left(\eta \right),{\psi }_2^N\left(\eta \right)=\sum _{i=0}^N{\beta }_i{\mathrm{SVL}}_i\left(\eta \right).$$

(12)

Now substituting (4) and (12) into the system (9)–(10) and evaluating them at the N points ${\eta }_s$ which are the roots of ${\mathrm{SVL}}_N\left(\eta \right)$, we obtain the identities

$$\begin{array}{cc}\hfill & \sum _{j=\lceil \varrho \rceil }^N{\alpha }_j{\chi }_{N,j,\varrho }[G_{\varrho ,p}({\eta }_s,{\gamma }_0)+G_{\varrho ,p}({\eta }_s,{\gamma }_m)\hfill \\ \hfill & +4\sum _{k=1,k-odd}^{m-1}{G}_{\varrho ,p}({\eta }_s,{\gamma }_k)+2\sum _{k=2,k-even}^{m-2}{G}_{\varrho ,p}({\eta }_s,{\gamma }_k)]\hfill \\ \hfill & =\rho -(1+\lambda )\left(\sum _{i=0}^N{\alpha }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)+{\left(\sum _{i=0}^N{\alpha }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)}^2\left(\sum _{i=0}^N{\beta }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right),\hfill \end{array}$$

(13)

and

$$\begin{array}{cc}\hfill & \sum _{j=\lceil \varrho \rceil }^N{\beta }_j{\chi }_{N,j,\varrho }[G_{\varrho ,p}({\eta }_s,{\gamma }_0)+G_{\varrho ,p}({\eta }_s,{\gamma }_m)\hfill \\ \hfill & +4\sum _{k=1,k-odd}^{m-1}{G}_{\varrho ,p}({\eta }_s,{\gamma }_k)+2\sum _{k=2,k-even}^{m-2}{G}_{\varrho ,p}({\eta }_s,{\gamma }_k)]\hfill \\ \hfill & =\lambda \left(\sum _{i=0}^N{\alpha }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)+{\left(\sum _{i=0}^N{\alpha }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)}^2\left(\sum _{i=0}^N{\beta }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right).\hfill \end{array}$$

(14)

Substituting Equation (12) into (11), the initial conditions (11) are transformed into the following system of algebraic equations:

$$\sum _{i=0}^{N}2{(-1)}^i{\alpha }_i={\psi }_{1,0},\sum _{i=0}^{N}2{(-1)}^i{\beta }_i={\psi }_{2,0}.$$

(15)

As a consequence, the following cost functions (CFs) can be used to express the preceding system as a restricted optimization problem:

$$\begin{array}{cc}\hfill CF1=& \sum _{r=0}^N\left|\sum _{j=\lceil \varrho \rceil }^N{\alpha }_j{\chi }_{N,j,\varrho }\right[G_{\varrho ,p}({\eta }_s,{\gamma }_0)+G_{\varrho ,p}({\eta }_s,{\gamma }_m)\hfill \\ \hfill & +4\sum _{k=1,k-odd}^{m-1}{G}_{\varrho ,p}({\eta }_s,{\gamma }_k)+2\sum _{k=2,k-even}^{m-2}{G}_{\varrho ,p}({\eta }_s,{\gamma }_k)]-\rho \hfill \\ \hfill & +(1+\lambda )\left(\sum _{i=0}^N{\alpha }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)-{\left(\sum _{i=0}^N{\alpha }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)}^2\left(\sum _{i=0}^N{\beta }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)|,\hfill \end{array}$$

(16)

and

$$\begin{array}{cc}\hfill CF2=& \sum _{r=0}^N\left|\sum _{j=\lceil \varrho \rceil }^N{\beta }_j{\chi }_{N,j,\varrho }\right[G_{\varrho ,p}({\eta }_s,{\gamma }_0)+G_{\varrho ,p}({\eta }_s,{\gamma }_m)\hfill \\ \hfill & +4\sum _{k=1,k-odd}^{m-1}{G}_{\varrho ,p}({\eta }_s,{\gamma }_k)+2\sum _{k=2,k-even}^{m-2}{G}_{\varrho ,p}({\eta }_s,{\gamma }_k)]\hfill \\ \hfill & -\lambda \left(\sum _{i=0}^N{\alpha }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)+{\left(\sum _{i=0}^N{\alpha }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)}^2\left(\sum _{i=0}^N{\beta }_i{\mathrm{SVL}}_i\left({\eta }_s\right)\right)|.\hfill \end{array}$$

(17)

In this case, the constraint has the form

$$\mathrm{Cons}=|\sum _{i=0}^{N}2{(-1)}^i{\alpha }_i-{\psi }_{1,0}|+|\sum _{i=0}^{N}2{(-1)}^i{\beta }_i-{\psi }_{2,0}|.$$

(18)

We use the Penalty Leap Frog procedure [24] for solving the constrained optimization problem (16)–(18) for the unknowns ${\alpha }_i$ and ${\beta }_i$, for each $i=0,1,2,...,N$. This, in turn, leads us to formulate the approximate solution by substitution into (12).

The remainder of this section is devoted to presenting numerical simulations on a test example associated with system (9)–(11) in the range $[0,3]$, with varying values of N, $\rho$, and $\lambda$, and initial conditions ${\psi }_{1,0}=1$ and ${\psi }_{2,0}=1$. As a result of our experiments, we confirm the accuracy and quality of the numerical scheme. Figure 1 illustrates the behavior of the approximation using several values of $\varrho$ (namely, $\varrho =1.0,0.9,0.8,0.7$), $\lambda =1$, $\rho =0.25$, and $N=8$. Figure 2 presents a comparison of our approach against the well-known fourth-order Runge–Kutta method (RK4), using $\varrho =1$, $\rho =0.25$, $\lambda =1$, and $N=8$. The behavior of the approximate solution employing the different values of $\lambda =0.5,1.0,1.5,$ along with $\varrho =0.95$, $N=8$, and $\rho =0.25$, is presented in Figure 3. Also, the behavior of the approximate solution using the values $\rho =0.3,0.6,0.9,$ and $\varrho =0.95$, $N=8$ and $\lambda =1.0$ is shown in Figure 4. From Figure 3 and Figure 4, we can confirm that the behavior of the BM depends critically on the values of the parameters $\lambda$ and $\rho$. These parameters control the strength of the reaction and diffusion processes and, hence, determine the system’s overall behavior. In turn, Figure 5 shows the residual error function (REF) of the approximate solution using the values $N=5,9,$ with $\varrho =0.95$, $\rho =0.1$ and $\lambda =1$. Finally, Figure 6 presents the behavior of the REF for the approximate solution using $\rho =0.2,0.4,0.6,$ along with $\varrho =0.95$, $N=8$, and $\lambda =0.5$. These results show the behavior of the numerical solution resulting from the application of the proposed method concerning the parameters $\varrho$, N, $\lambda$, and $\rho$. We observe that the proposed method is adequate for solving the proposed model in its fractional form with the RFE operator.

4. Conclusions

The Brusselator model is a simple yet powerful mathematical model of chemical reaction dynamics that has significantly impacted the field of chemical research. By understanding the behavior of this model, researchers can gain insights into the underlying mechanisms of chemical reactions and develop new strategies for controlling and manipulating these reactions. This model is an important example of a nonlinear dynamical system and has applications in various fields beyond chemistry. In this manuscript, the numerical solutions for the fractional Brusselator model were calculated using the Rabotnov fractional-exponential kernel fractional derivative. Various fractional orders were considered, and we employed the residual error function in our investigation. The results confirmed that the proposed method is suitable for approximating the fractional mathematical model considered in this work. Additionally, by including additional terms from the approximation solution series, we were able to control the precision of the error and decrease it conveniently. Moreover, we concluded that the numerical simulations of the model under consideration in the current study are better suited for the RFE operator without singularity. By assessing the residual error function, we confirmed the effectiveness and accuracy of the proposed approach. As a future generalization of this work, we will attempt to deal with the same problem using a different kind of fractional derivative and different type of polynomials. Finally, we must point out that the Mathematica software program was used to carry out the numerical simulation work.