Fractional Clique Collocation Technique for Numerical Simulations of Fractional-Order Brusselator Chemical Model

: The primary focus of this research study is in the development of an effective hybrid matrix method to solve a class of nonlinear systems of equations of fractional order arising in the modeling of autocatalytic chemical reaction problems. The fractional operator is considered in the sense of Liouville–Caputo. The proposed approach relies on the combination of the quasi-linearization technique and the spectral collocation strategy based on generalized clique bases. The main feature of the hybrid approach is that it converts the governing nonlinear fractional-order systems into a linear algebraic system of equations, which is solved in each iteration. In a weighted L 2 norm, we prove the error and convergence analysis of the proposed algorithm. By using various model parameters in the numerical examples, we show the computational efﬁcacy as well as the accuracy of our approach. Comparisons with existing available schemes show the high accuracy and robustness of the designed hybrid matrix collocation technique.


Introduction
The Brusselator is a theoretical model for a type of autocatalytic reaction. In fact, this model is a common nonlinear reaction in which a reactant species interacts with other species to increase its production rate. The Brusselator model was proposed by Prigogine and Lefever [1] in 1968. It is also known that the Belousov-Zhabotinsky model and the chemical reactions of the Brusselator are the same [2][3][4]. By U, V, D, A, B, and E, we denote the chemical components in the chemical reaction. Generally, the reaction process can be described by the following four steps: We now assume that the species A and B are sufficiently available and can thus be modeled at a constant concentration. Further, note that the final products E and D are removed once they are produced from the reaction process. Under scaling the rate constant to unity, the rate equations become as follows Fractional integrals and derivatives have attracted considerable attention over the last decades. Due to a wide range of applications from theory to practice, they have gained increasing popularity in the modeling of various natural phenomena in engineering, physics, chemistry, economics, etc. It is found that the non-integer derivatives and integrals are more appropriate for describing the properties of several real processes and materials, see cf. [5,6]. However, the solutions to most fractional differential equations do not exist in terms of elementary functions. Therefore, it is essential to develop computational and approximate procedures for the numerical evaluation of fractional differential equations.
Our main goal is to study the fractional counterpart of the Brusselator model (1). To be precise, this research paper presents a power series solution based upon the (fractional) version of clique functions implemented in matrix formulation for the following nonlinear fractional-order Brusselator system of two equations where θ and η are two positive real numbers. Moreover, LC D λ τ presents the Liouville-Caputo fractional derivative of order λ ∈ (0, 1]. The following initial conditions will accompany the above system, given as If we set λ = 1, the classical system of the Brusselator system (1) will be obtained. The integer-order model of the Brusselator has been solved by three numerical approaches, including the Implicit Runge-Kutta method, the Adams method, and the Backward differential formula in [7].

Literature Review and Related Works
The fractional-order system (2) has been considered in the literature by many research scholars from different points of dynamic systems and numerical behaviors. The stability of the fractional Brusselator system was addressed in [8][9][10]. In [11], the existence of a limit cycle was proven numerically by the Adams-Bashforth-Moulton approach. The authors of [12,13] developed some nonstandard finite difference (NSFD) methods to solve (2) numerically. As a semi-analytical approach, the variational iteration method is devised in [14]. The polynomial least square technique was investigated in [15]. The operational matrix methods based on Bernstein and Legendre wavelet functions were studied in [16,17], respectively. The authors of [18] further developed three explicit and implicit techniques based on product integration, NFSD, and multi-step procedures. In all mentioned works above, the underlying fractional operator was taken as the Liouville-Caputo fractional derivative. However, let us mention that the Brusselator model with fractional derivatives in the sense of Liouville-Caputo, Caputo-Fabrizio, and Atangana-Baleanu was considered in [19] recently. In this paper, the dynamic characteristics of the model under three fractional derivatives have been investigated, and a three-stage iterative approach was also developed for the model under consideration. In addition, in [20], the fractal-fractional differential operators related to the power law, exponential decay, and the generalized Mittag-Leffler kernels were investigated. In the latter research work, the proposed numerical procedure is based on the Lagrange interpolating polynomial together with the theory of fractional calculus. Finally, let us mention that the PDE counterpart of the Brusselator model (2) has been investigated in the literature. Among others, we refer to the published papers [21][22][23][24]

Outline of This Paper
The primary purpose of the current research paper is to propose an effective hybrid technique. Our novel method is based on a combination of the quasi-linearization approach and matrix collocation method for an approximate treatment of the fractional Brusselator equations. The idea of a quasi-linearization method (QLM) is used to convert the nonlinear model into a family of linearized equations. Afterward, the spectral approach based on the (novel) generalized clique functions (GCFs) is employed to solve the quasi-linear equations in an iterative manner. Let us emphasize that the coefficients of all clique polynomials are all positive and integer-compared to the classical set of polynomials, such as Legendre, Chebyshev, Hermit, Laguerre, etc. Consequently, working with positive numbers yields more stable results during the computations. This would be the main motivation to employ the family of clique polynomials in the collocation matrix procedure over others. Another major advantage of the presented hybrid method, namely QLM-GCFs, is that it is not only effective in terms of required CPU time, but it provides high-order accuracy and better resolution characteristics compared to the existing numerical models in the literature. The accurateness and robustness of the spectral collocation strategies have been justified successfully by applying various model equations. Among others, let us mention the works [25][26][27][28][29][30].
The content of this research paper is organized as follows. Some basic facts on fractional calculus are reviewed in Section 2. Section 3 is devoted to the definition of clique basis functions. Moreover, a generalization of these functions is given. Then, the convergence analysis of this basis function is established in a weighted L 2 norm. A detailed description of the present QLM-GCFs technique is provided in Section 4. The results of the performed numerical simulations and experiments are given in Section 5. The concluding summary is given in Section 6.

Fractional Calculus: Basic Facts
Let us give some important facts about fractional calculus that will be used in the subsequent sections. For more detail, we refer the readers to the standard text [6] or some recent expository papers [31,32].
Let us first recall that the Riemann-Liouville fractional integral operator of order λ > 0 is given by where Γ(·) is the Gamma function and we assumed k(τ) ∈ C ξ , ξ > −1. We note that a real function k(τ), τ > 0 belongs to the space C ξ , ξ ∈ R if there exists a number µ ∈ R and a function l(τ) ∈ C ∞ ([0, ∞) such that k(τ) = τ µ l(τ). We also call that k(τ) ∈ C n ξ if and only if k (n) (τ) ∈ C ξ for a n ∈ N.
We are now ready to define the fractional Liouville-Caputo derivative next.

Definition 1.
Assume that k ∈ C n −1 and n − 1 < λ < n, n ∈ N. The Liouville-Caputo fractional derivative of k(τ) of order λ is defined by One should emphasize that the fractional operator LC D λ τ is a linear operator. If C is a constant number, we have LC D λ τ C = 0.
Our next goal is to compute the Liouville-Caputo fractional derivative of the function k(τ) = τ p , where p is a constant. This can be performed through the following relations τ p−λ , for p ∈ N 0 and p ≥ λ or p / ∈ N and p > λ . (5) Note that N 0 := N ∪ {0} and also we have utilized · and · as the ceil and floor functions respectively.

The Fractional-Order Clique Polynomial and Its Convergence Analysis
Here, we first consider the clique polynomials C r (t) related to the cliques in a complete graph. We then introduce the fractional version of these polynomials. Hence, we establish the convergence analysis of these polynomials.

The Clique Functions: The Generalized Form
The clique polynomials were first introduced in [33] and associated with graph theory. However, they have been recently considered for numerical approximations of ordinary and fractional differential equations, see cf. [34][35][36][37]. Below, we first describe the main aspects of them.
By using recursion (7), we derive a few terms of CFs as One can easily check that C r (0) = 1 and C r (1) = 2 r for all values of r ≥ 0. It is not a difficult task to check that these CFs satisfy a second-order differential equation in the form In what follows, we intend to use the CFs on an arbitrary interval D a,b := [a, b]. We are also interested in using the generalized version of these polynomials of fractional order 0 < α ≤ 1.

Definition 3.
Generalized CFs (GCFs) of degree r on D a,b are represented by C α r (τ) and defined by With the help of this transformation, the explicit form in (6) will be given as follows:

L 2 -Convergent of GCFs
Let us investigate the convergence analysis of the GCFs in a weighted L 2 norm. In other words, we will investigate the behavior of the expansion series of a given function with respect to GCFs, especially when we increase the number of bases. We associate the following space to domain D a,b as [25] Here, the related induced norm and inner product are given by Practically, a finite-dimensional subset, say S R , of the space L 2,w (D a,b ) is selected as One observes that dim(S R ) = R + 1 and is also a closed subspace of L 2,w (D a,b ). It follows that S R is a complete subspace of L 2,w (D a,b ). Thus, any given function ∈ L 2,w (D a,b ) has a unique best (finest) approximation ∈ S R in the following sense that Generally, a given function (τ) ∈ L 2,w (D a,b ) can be expressed as a linear combination of GCFs. Thus, we have Here, the coefficients of κ r are unknown for r = 0, 1, . . .. By truncating the former series up to (R + 1) terms, we may approximate (τ) in practice as where we have represented the involved finite series in a compact way by defining Note that the first one is the vector of GCFs while K K K R is the vector of unknowns. To derive an upper bound for the E R (τ) = (τ) − R (τ), we need the following result. A proof of the next theorem can be found in [38].
). If l R (τ) represents the interpolating function of at R Chebyshev points on D a,b , then we have Following [39], we establish the following error bound for the GCFs expansion series.
presents the best (finest) approximation of (τ) out of S R,α , then an error bound is given by Proof. Let us first define the new function z(t) : and for any α > 0. Applying Theorem 1 to function z(t) with R Chebyshev nodes leads to the following error estimate We next substitute t = τ α in the preceding inequality. It follows that According to the theorem's assumption, we know that the approximate solution R (τ) is the finest approximation belonging to the space S R,α . Thus, it holds that The former inequality is valid, particularly for h = l R (τ α ) ∈ S R,α . Employing this fact, as well as (15), we conclude that Using the fact that b a w(τ)dτ = 1, we only require the application of the square roots to the foregoing inequality.

The Methodology of the QLM-GCFs Scheme
Instead of applying the direct collocation procedure to the underlying nonlinear model (2), our main aim is first to employ the quasi-linearization method (QLM) for (2) with initial conditions (3). Then, the generalized CFs (GCFs) collocation matrix technique is applied to the resultant family of linear subequations in an iterative manner.

The Basic Concept of QLM
By employing QLM, we can overcome the nonlinearity of a given model equation. The applicability of the QLM strategy has already been checked through tremendous research studies in the literature. For recent applications, we refer readers to [40][41][42][43].
By rewriting first the nonlinear coupled system (2) in a matrix representation form, we get Here, we have utilized the following notations Suppose that Z Z Z 0 (τ) is a rough first approximation of Z Z Z(τ). Thus, the QLM for (17) is written for p = 0, 1, . . . as Here, by F F F Z Z Z = d dZ Z Z F F F, we denote the corresponding Jacobian matrix. Let us note that the same initial conditions as (3) will be given to the last sequence of equations. After performing some straightforward calculations, we receive the following family of a linear system of equations as the result of QLM from model (17). Thus, we have where Systemically, we can present the initial conditions (3) as To solve the quasi-linear systems (18)-(19) accurately, we will design a matrix collocation procedure based on the GCFs to receive an approximate solution.

The QLM-GCFs Technique
Supposedly, the unknown solutions of quasi-linear model (18) can be expanded as a combination of the cut series form (13) with (R + 1)-terms. Further, assume that for a fixed α ∈ (0, 1] the approximate solutions U R,α (τ) to u p (τ) and v p (τ) in the iteration p for p = 0, 1, . . . are known. For p = 0, we utilize the initial guess Z Z Z 0 (τ) as the starting point. In the next iteration, p + 1, we seek the approximate solutions in the forms (20) for τ ∈ D a,b . Below, our primary job is to find the unknown coefficients {κ (p) r,j } R r=0 for j = 1, 2 and p = 1, 2, . . . by using a spectral matrix collocation approach relying on the GCFs. To this end, we first construct the matrix forms of the approximate solutions given in 20. Similar to (13), the finite series solutions in ( (20)) for j = 1, 2 can be expressed as where the unknown vectors K K K R,j and the vector of GCFs are given by In the next Lemma, we further write the vector of basis functions in terms of monomials multiplied by a constant matrix. Lemma 1. The representation of the vector of GCFs is given by where Π Π Π α R (τ) = 1 τ α τ 2α . . . τ Rα and M M M R of size (R + 1) × (R + 1) is an upper triangular matrix whose components are obtained via (10). It reads Proof. By virtue of relation (10), it is sufficient to multiply matrix M M M R by Π Π Π α R (τ) from the left.
Obviously, the diagonal elements of matrix M M M R are all non-zero. Thus, this matrix is non-singular. In fact, we have det If one combines two former relations (21) and (22), the approximate solutions are written as U Next, we need the λ-derivative of the approximate solutions. To do so, we apply the operator LC D λ τ to both sides of the relation (23). Thus, we get Consequently, we must only compute the λ-derivatives of the vector Π Π Π α R (τ). In this respect, we consider two properties (4) and (5) to calculate the fractional derivatives of Π Π Π α R (τ). As an example, we set R = 3 and λ = 3 4 . Now, using α = 1 and α = 1 2 we get, respectively, Practically, however, we can use the modified version of Algorithm 4.1 in [44] or [45] or Algorithm 3.1 in [46] with linear complexity O(R + 1) to compute the λ-derivative of Π Π Π α R (τ). To continue, let us define the fractional derivative of the vector as We can now place this relation into (24) to arrive at We come back to the matrix differential equation (18). Vector Z Z Z p+1 (τ) and its derivative LC D λ τ Z Z Z p+1 (τ) can be approximated as Lemma 2. In the matrix formats, the approximated solution Z Z Z (27) can be stated as follows: where the following notations are used: Proof. To conclude the results, we need to substitute two relations (23) and (26) into the corresponding vector forms in (27).
We are then looking for a partitioning D a,b that will be used as a set of collocation points. To do so, we utilize (R + 1) equidistant points from interval [a, b]. Let us set Now, adding the aforementioned set of collocation points into the sequence of linear matrix differential equations (18) to get for p = 0, 1, . . .. We next introduce the following matrix and vector notations In the vector representation, we are able to show relation (30) in a compact formulation as Our next aim is to derive the matrix expressions of Σ Σ Σ p and Σ Σ Σ λ p . By collocating two relations (28) at the collocations points, we render Lemma 3. The two relations (28) at the collocation point (29) can be written as follows: where the two matrices Π Π Π and Π Π Π λ are given by Here, the three matrices M M M, Π Π Π, and Π Π Π λ , as well as the vector K K K (p) , are defined in (28) previously.
Finally, we form the so-called fundamental matrix equation at each iteration p by placing the preceding relations (32) into (31). It follows that where It should be noted that the last matrix equation (33) is a linear system with 2(R + 1) unknowns κ (p) r,j for r = 0, 1, . . . , R and j = 1, 2 to be specified as the coefficients of GCFs in the series solutions (20). However, the supplemented initial conditions (3) are not yet implemented and entered into the system (33). First, we consider the matrix representation forms (28) for the approximate solution Z Z Z (p+1) R (τ). We then let τ → 0 arrive at Here, the two constants u 0 and v 0 , are available from (3). The replacement of two rows of the matrix [A A A p ; S S S p ] in (33) will be carried out next by the row matrix A A A 0,p ; S S S 0 . The modified fundamental matrix equation will be shown by To get the unknown coefficients of GCFs, it is sufficient to solve the modified algebraic linear system (34) in each iteration. Now, one requires a linear solver to be used to receive the solution of this system. After finding vector K K K (p) , all unknowns κ (p) r,j , for j = 1, 2, and r = 0, 1, . . . , R as the coefficients in the expansion series (20) will be found in iteration p. Thus, we get an approximate solution of model (2).
Algorithmically, we summarize all of the steps of the proposed QLM-GCs technique in Algorithm 1. Here, by p max we denote the maximum number of iterations required to achieve the desired accuracy in the QLM method. It should be remarked that we have utilized the MATLAB notation A A A[i : j, s : k] to denote the submatrix of A A A formed by all entries in the intersection of rows i, . . . , j and columns s, . . . , k.

Numerical Results and Graphical Representations
In this part, a set of computational examples is provided to describe and support the theoretical findings. In this respect, we apply the QLM-GCFs to the fractional-order Brusselator of Equation (2) by solving the quasi-linear model Equation (17). For performing computational simulations, we use Matlab software version 2021a on a computer with 16 GB of RAM and a CPU with 2.2 GHz Intel® Core™ i7-10870H processor.
In the computational results, we utilize the QLM with parameter p = 5. Furthermore, in the QLM-GCFs, the initial approximation Z Z Z 0 (τ) is selected as the zero function, or we take it as the initial condition (3). As previously mentioned, the exact solutions of this system are not available, especially when the order of the derivative is described in the fractional order. Therefore, we define the residual error functions (REFs) associated with the Brusselator model to measure the accuracy of the proposed spectral QLM-GCFs collocation technique. That is, in iteration p = 1, 2, . . ., we define the error terms as We also compute the L ∞ error norms (for a fixed p) via the relations We further utilize the following relations to compute the obtained numerical order of convergence (Noc) related to the numerical technique applied to both solutions of the coupled system (2) given by Note that these formulae are utilized to check the order of accuracy of our proposed technique in the L ∞ norm for both solutions.

Example 1.
As the first test case, let us consider the fractional Brusselator system by taking two parameters θ = 0 and η = 1 to get The given initial conditions are u(0) = 1, v(0) = 1. This example was considered in [14,16,17] previously.
The second method is the Legendre wavelet operational matrix method (LWOMM) [17], the solutions of which are reported as In Figure 1, we show the above approximate solutions obtained by our method (black lines) and two other existing ones, i.e., the PLSM and LWOMM procedures. From this visualization, we conclude that the alignment between the results of QLM-GCFs and PLSM is more than the outcomes of our method and LWOMM. On the other hand, note that our solutions are obtained by using R = 5, which gives us the approximate polynomial solutions of five degrees compared to the three-degree polynomials reported by the LWOMM and PLSM. However, our proposed procedure can produce more accurate results just by increasing R. To be more precise, we plot the achieved REFs obtained via (35) when R = 5, 10, and 15. These experiments are shown in Figure 2.      Finally, for the integer order λ = 1, we report the numerical results evaluated at some points τ ∈ [0, 1]. For this purpose, we use R = 10 and show the outcomes of the proposed QLM-GCFs for both solutions in Table 1. Table 2 presents the maximum REF values achieved by using R = 2 i , i = 1, 2, 3, 4, 5. The corresponding Noc are also reported in this table related to both solutions u(τ) and v(τ). Higher order accuracy of the proposed method is visible from the results presented in Table 2. The required CPU times to solve modified system A A A p ; S S S p measured in seconds are shown in Table 2. The CPU's spent time clearly behaves linearly as the number of basis functions becomes two-fold.
Let us turn next to the fractional cases and set λ = 0.75. By employing the QLM-GCFs with R = 10, we obtain two approximate solutions for 0 ≤ τ ≤ 1 as given below. The obtained results for α = 1 are given by 10,1 (τ) = 2.2725108 τ 10 − 13.632576 τ 9 + 36.11913 τ 8 − 55.727909 τ 7 + 55.641481 τ 6 − 37.881744 τ 5 + 18.156988 τ 4 − 6.2694104 τ 3 + 1.4527369 τ 2 + 0.079421181 τ + 1.0.   The former approximations for each solution of u(τ) and v(τ) related to two different values of α = 1 and α = 0.75 are depicted in Figures 3 and 4. In addition to the approximate solutions, we also visualize the associated REFs for each solution on the left plots. By looking at these figures, we infer that the approximate solutions related to both α = 1, 0.75 are very close together. However, the achieved REFs for α = 0.75 equal to the fractional order λ = 0.75 are smaller in magnitude than those obtained using α = 1. Therefore, in the next experiments, we only consider the results obtained related to α = λ.
Next, we consider λ = 0.5 and R = 10. Using α = 0.5, the approximate solutions evaluated at some point τ ∈ [0, 1] are reported in Table 3. The corresponding absolute errors defined via relations (35) are also tabulated in Table 3.  When λ = 0.98, different numerical methods, such as the variational iteration method (VIM) [14], the PLSM [15], the method based on an operational matrix of Bernstein polynomials [16], and LWOMM [17], reported the approximate solutions for this value. In all of these approaches, the solutions obtained are three-degree polynomials. However, here, we first consider the case α = 1 and obtain the following approximate solutions Additionally, the maximum absolute values of REFs using α = 0.98 are shown in Table 4 for various values of R as a power of 2. The related Nocs are also tabulated in this table. Moreover, we present the numerical results when λ, α = 0.75 in Table 4. One can obviously observe a high order of accuracy for the proposed QLM-GCFs. Finally, for the first test case, we use various values of λ = 0.25, 0.5, 0.75, 1. Utilizing R = 10 and α = λ, we plot the numerical solutions in Figure 5. Table 3. Numerical results and REFs for u(τ), v(τ) obtained via the QLM-GCFs procedure using R = 10, p = 5 in Example 1 with λ, α = 0.5. τ U (5) 10,0.5 (τ) Res (5) u,10,0.5 (τ) V (5) 10,0.5 (τ) Res (5)     Example 2. In the second model problem, we take θ = 0.5 and η = 0.1. In this case, we consider the nonlinear coupled system Here, we use initial conditions u(0) = 0.4, v(0) = 1.5. This example was considered in [14,16,17] previously.
Let us first take λ, α = 1. By using R = 5, the results of the approximations are as follows The REFs related to the above approximate solutions are shown in Figure 6. To show that the achieved REFs are decreasing as a function of R, we also plot the REFs related to R = 10, 15 in Figure 6. Clearly, the desired level of accuracy is achievable by increasing the number of basis functions.
We next tabulate the numerical results obtained by using R = 10 in Table 5. Here, we have used the midpoints 0.05, 0.15, . . . , 0.95 on the interval [0, 1]. Note that these midpoints are different from the points used in Table 1, which are exactly the same as the collocation points (29) when R = 10. In fact, the smallest magnitude of errors is achieved at the collocation points. Finally, for this test case and for λ = 1, we display the maximum absolute REFs achieved by utilizing various R numbers in Table 6. The associated Nocs are also visible in this table. The results show the exponential behavior in terms of the accuracy of the presented QLM-GCFs.

Conclusions
Generalized (fractional-order) clique basis functions (GCFs) have been used to devise not only an effective but also an accurate spectral matrix collocation approach for finding approximate solutions of the nonlinear Brusselator system of equations of fractional order arising in chemical modeling. The fractional derivative is described in the Liouville-Caputo sense. To overcome the underlying nonlinearity of the model, the method of quasi-linearization (QLM) is first employed to receive a family of linearized equations. Afterward, the spectral clique collocation procedure is used to solve this sequence of equations iteratively. The convergence analysis of the proposed combined QLM-GCFs is established. To support the theoretical findings and in order to show the applicability of the QLM-GCFs, a set of numerical test examples is carried out. The results presented in the tables and figures indicate the accuracy of the proposed approach over the existing numerical models and the gain in computational efficiency in terms of CPU time. The presented technique is straightforward, easy to implement, and computationally less demanding.