EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 4, Article Number 7064 ISSN 1307-5543 – ejpam.com Published by New York Business Global Using the Multi-Domain Spectral Relaxation Method through a Numerical Simulation for the Brusselator System M. Adel1,∗, M. M. Khader2,3, Hijaz Ahmad4, Osama Oqilat5, W. M. Abdelfattah6,7 1 Department of Mathematics, Faculty of Science, Islamic University of Madinah, Madinah, 42351, Saudi Arabia 2 Department of Mathematics and Statistics, College of Science, Imam Mohammad Ibn Saud Islamic University (IMSIU), Riyadh, Saudi Arabia 3 Department of Mathematics, Faculty of Science, Benha University, Benha, Egypt 4 Operational Research Center in Healthcare, Near East University, Nicosia/TRNC, 99138 Mersin 10, Turkey 5 Department of Basic Sciences, Faculty of Arts and Science, Hourani Center for Applied Scientific Research, Al-Ahliyya Amman University, Amman, Jordan 6 College of Engineering, University of Business and Technology, Jeddah 23435, KSA 7 Department of Engineering Mathematics and Physics, Faculty of Engineering, Zagazig University, P.O. 44519, Egypt Abstract. We introduce the numerical solutions of the Brusselator system by applying the multi- domain spectral relaxation method (MSRM). The proposed method is a combination of the Cheby- shev pseudo-spectral method and the Gauss-Seidel relaxation approach. This method breaks down the main interval into several small subintervals, finding solutions within each interval. This ap- proach also transforms the model into a set of algebraic equations. We present the error analysis and the order of convergence for the proposed scheme. We validate the effectiveness and preci- sion of the given procedure by utilizing the fourth-order Runge-Kutta method (RK4M), and the variational iteration method. The considered results are tabularly and graphically demonstrated with different values of the model’s parameters. From a numerical viewpoint, the given simulations and results indicate that the proposed algorithm is a straightforward and appropriate tool with computational efficiency for such models. 2020 Mathematics Subject Classifications: 41A30, 49M27, 65N20 Key Words and Phrases: Brusselator system, MSRM, Chebyshev pseudo-spectral method, Gauss-Seidel relaxation approach, convergence analysis, RK4M ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i4.7064 Email addresses: adel@sci.cu.edu.eg (M. Adel), mmkhader@imamu.edu.sa (M. M. Khader), hijaz.ahmad@neu.edu.tr (H. Ahmad), o.oqilat@ammanu.edu.jo (O. Oqilat), w.abdelfattah@ubt.edu.sa (W. M. Abdelfattah) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 2 of 11 1. Introduction While spectral methods are recognized for their precision, their accuracy diminishes for non-smooth solutions and extensive domain issues, even when more grid points are added ([1], [2]). Domain decomposition has mainly been used with semi-analytical tech- niques to address the numerical solutions of some systems. Such techniques encompass methods like the multistage Adomian decomposition method [3], multistage variational iteration method [4], and multistage homotopy perturbation method [5]. However, a chal- lenge faced by multi-domain techniques based on analytical approximations is the long and complex process of analytical integration within each sub-domain. Motsa and col- leagues ([6], [7]) introduced a spectral method-based multi-domain numerical approach, setting it apart from earlier multi-domain strategies by being entirely numerical. Other numerical strategies that have employed domain decomposition for chaotic systems are the time discretization [8], piecewise successive linearization [9], piecewise spectral homotopy analysis [10], and multi-domain compact finite difference relaxation ([11], [12]). Also, there are numerical algorithms can be applied in many real life problems as in the operation research, data analysis and the financial mathematics ([13], [14]), also in the simulations of some biological models as in ([15],[16]). In this work, we describe the Brusselator mathematical model, which predicts the oscil- lations in chemical reactions. We present an argument for why modern thermodynamics is important in analyzing a system that constantly interacts with its environment, and which operates far from thermodynamic equilibrium. The Brusselator model has been used in a wide range of applications in chemical research, where it has been used to study the formation of spatial patterns in chemical systems, such as Turing patterns and spiral waves. It is considered a theoretical model for a type of self-catalytic reaction [17]. It is also considered an example of an autocatalytic, oscillatory chemical reaction, that is, a reaction in which a species increases the rate of its productive reaction [18]. In this paper, we aim to determine and provide the numerical solutions to the Brussela- tor system using the MSRM. This technique incorporates the Chebyshev spectral method ([19], [20]) and the Gauss-Seidel method process to numerically treat the system across several sub-domains that constitute the entire problem domain. We can confidently assess the effectiveness of the given numerical approach for the model under study. The nu- merical solutions validate that the proposed technique may effectively tackle the specified model, aligning well with existing solutions. By adding terms from the solution series, we can further minimize relative errors. We highlight the efficiency and vast potential of our proposed numerical method by contrasting the estimated solutions with those ob- tained from the RK4M. We also illustrate that the multi-domain strategy’s strength lies in its reduced error accumulation across sub-domains compared to considering a single domain. The proposed method is a powerful technique that uses an approach of lineariza- tion of nonlinear equations for obtaining a better series solution, which provides efficient algorithms for approximate solutions, giving effective precision and convergence. Besides, the solutions converged in a relatively narrow region, and they had sluggish convergence. From the novelty points, the proposed method is considered for the first time for solving M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 3 of 11 the Brusselator system. 2. Basic concepts 2.1. Mathematical formulation Here, we focus on the Brusselator system, which is formulated as follows [21]: θ̇1(t) = δ − (1 + γ) θ1(t) + (θ1(t)) 2 θ2(t), (1) θ̇2(t) = γ θ1(t)− (θ1(t)) 2 θ2(t), (2) with the following initial conditions: θ1(0) = θ01, θ2(0) = θ02, (3) where θ1 and θ2 represent the concentration of the reaction and diffusion processes, re- spectively, the parameters δ > 0, γ > 0, and θ01, θ 0 2 are constants [22]. The system can exhibit sustained oscillations as δ increases beyond this critical value. The parameter γ controls the balance between the reaction and diffusion processes. The reaction dominates at low values of γ, and the system can exhibit sustained oscillations or spiral waves. As γ increases, diffusion becomes more important, and the system can exhibit Turing patterns. 2.2. The Chebyshev-Lagrange interpolation The Chebyshev polynomials Tk(y) of order k are generated from the following formula: Tk(y) = cos (k arccos(y)) , k ∈ N. (4) The Chebyshev-Lagrange interpolation of the function ψ(y) at y = yi, i = 0, 1, ...,m is denoted by ψm(y) and defined as follows: ψm(y) = m∑ i=0 ψiLi(y), (5) where Li(y), i = 0, 1, ...,m are the Lagrange’s polynomials of degree i. The extreme yi of Tm(y), which are the Chebyshev-Gauss-Lobatto points (CGLPs), are chosen to be the collocation points, and defined by: {yi}mi=0 = { cos ( πi m )}m i=0 . (6) Remark 1. This previous choice was made based on our knowledge that the error in Lagrange interpolation is minimized if the interpolation points are zeros of the polynomials. M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 4 of 11 3. Numerical implementation of the MSRM The suggested approach will be used to obtain the numerical solutions of the investi- gated model with the help of the following steps (Algorithm): (i) Without losing the system’s generality, we describe the system (1)-(2) as follows [23]: ξ̇ + C ξ +N(ξ) = 0, (7) where ξ(t) = [θ1(t), θ2(t)] T ; C is an 2 × 2 matrix with entries ci,j , i, j = 1, 2 and is defined by: C = ( 1 + γ 0 −γ 0 ) , here the elements c1,2 and c2,2 are equal zero, because the equations (1)-(2) do not contain on linear terms of θ2(t). Also N(ξ) is a vector of nonlinear terms of Eq.(7) and is defined by: N(ξ) = ( −δ − (θ1(t)) 2 θ2(t) (θ1(t)) 2 θ2(t) ) . (ii) We define the following set of Lagrange polynomials {Li(y)}mi=0 of order i, which are based on the CGLPs as follows: Li(y) = (−1)i+1 ( 1− y2 ) Tm(y) ℏim2 (y − yi) , i = 0, 1, . . . ,m, (8) where ℏ0 = ℏm = 2, ℏi = 1 for i = 1, 2, . . . ,m− 1. Also, the Chebyshev polynomial Tm(y), and the CGLPs yi are defined in (5) and (6), respectively. (iii) We use the Chebyshev-Lagrange interpolation of the function v(y) at y = yi, i = 0, 1, ...,m, which is denoted by vm(y) and defined as follows: vm(y) = m∑ i=0 viLi(y), (9) where the coefficients vi = v(yi). (iv) We calculate the first derivative of the approximate solution (9) at the collocation points yi as follows: dvm(y) dy = m∑ k=0 v (yk) dLk (yi) dy = m∑ k=0 Di k v (yk) = DVi, i = 0, 1, . . . ,m, (10) the elements of the first-order Chebyshev differentiation matrix Dik can be evaluated for i, k = 0, 1, . . . ,m at the collocation points yi (6) as follows [24]: Dik = dLk (yi) dy =  ℏi(−1)k+i ℏk(yk−yi) , i ̸= k, − yi (1−y2i ) , (i = k) ̸= 0,m, 2m2+1 6 , i, k = 0, −2m2+1 6 , i, k = m. (11) M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 5 of 11 (v) We divide the interval I = [0, T ] into n subdomains with a uniform length T n = tj − tj−1, where Ij = [tj−1, tj ] has the characteristic that: n⋃ j=1 [tj−1, tj ] = [0, T ]. (12) Now, by using the following linear transformation, we can convert each subdomain of the form [tj−1, tj ] to the domain of the CGLPs [−1, 1], which is specified in Eq.(6): t̂ = 0.5(tj − tj−1)y + 0.5(tj−1 + tj) = 0.5(T/n) y + 0.5(tj−1 + tj), j = 1, 2, . . . , n. (vi) We estimate the derivatives of the first term in the system of ODEs (7) with the help of the Chebyshev differentiation matrix, to reduce it to a system of algebraic equations. We use Θj s = [ θs ( yj0 ) , θs ( yj1 ) , · · · , θs ( yjm )]T , and yjk, k = 0, 1, ...,m to represent the approximate solution and the collocation points in each subdomain, respectively. These considerations lead us to the following matrix form system:( 2n T D+ cs,s I ) Θj s + 2∑ i=1 i̸=s cs,iΘ j i +Nj s(Θ) = 0, s = 1, 2, (13) where cs,i are the elements of the matrix C defined in (7), and: Θ = [Θ1,Θ2] T , Ns(Θ) = [Ns (Θ (y0)) , Ns (Θ (y1)) , . . . , Ns (Θ (ym))]T . (vii) By using the Gauss-Seidel method to solve the algebraic equations (13) in each subdomain Ij , we can reduce it to the following iteration formula:( 2n T D+ cs,s I ) Θj s,p+1 = − s−1∑ i=1 cs,iΘ j i,p − 2∑ i=s+1 cs,iΘ j i,p −Nj s,p, s = 1, 2, (14) where Θs,p+1 is the numerical estimation of θs at the (p+1)th iteration, and Nj s,p = Ns ( Θj 1,p, Θ j 2,p, . . . ,Θ j s−1,p, Θ j s,p ) are used to evaluate the nonlinear terms. (viii) In light of the above, the numerical solution can be obtained as follows: Θj s,p+1 = A−1 s Bj s , (15) where As = 2n T D+ cs,s I, Bj s = − s−1∑ i=1 cs,iΘ j i,p − 2∑ i=s+1 cs,iΘ j i,p −Nj s,p. We will use the initial conditions in the first iteration at p = 0, of the iteration matrix form (15). The solution of the system (7) is then given by [25]: θs = n⋃ j=1 Θj s (yj) . (16) M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 6 of 11 4. Convergence and error analysis We describe how domain decomposition affects MSRM convergence and error esti- mates. Eq.(15) can be expressed as a block matrix system as follows: AΘ = B, (17) where A = ( 2n T D+ c1,1 I c1,2 I c2,1 I 2n T D+ c2,2 I ) , Θj = ( Θj 1 Θj 2 ) , and Bj = ( −Nj 1 −Nj 2 ) . Theorem 1. The sequence of solutions { Θ(k) }∞ k=0 , which is generated from (15) by Θ = A−1B with any choice of Θ(0) ∈ Rn̄, converges to the unique solution Θ̄ if the matrix A is strictly diagonally dominant. Proof. The details of the proof of this theorem can be found in [26]. Theorem 2. [27] Let Θs(t) be a polynomial of degree ≤ m that interpolates the function θs(t) ∈ Cm+1[0, T ] at m+ 1 various points y0, y1, . . . , ym ∈ [tj , tj+1], with the property (12). According to (8), if the nodes yi are used as the CGLPs, the error term for the polynomial interpolation in [0, T ] using the nodes yi in each sub-domain is estimated as follows: E(t) = |θs(t)−Θs(t)| ≤ λTm+1 (4n)m(m+ 1)! , λ ̸= 0. (18) Remark 2. It can be shown that the error in the multi-domain situation is significantly smaller than the single domain example, n = 1, due to the factor ( 1 n )m << 1 for large n. 5. Numerical simulation By addressing the system (1)-(2) in the range [0, 3] with various values of m, δ, γ with initial conditions θ01 = θ02 = 1, we will show a numerical simulation on a test example to demonstrate the accuracy of the provided numerical scheme. Figure 1 (a,b) shows the behavior of the obtained approximate solutions using different values of the order of the used polynomials, m = 3, m = 7, at δ = 0.25, γ = 1. We compared the results to those obtained using the RK4M (as a reference solution) with δ = 0.0, and γ = 1 in Figure 2 (a,b). This comparison demonstrates that the proposed method is appropriate for solving the suggested model. The behavior of the approximate solution via distinct values of γ = 0.5, 1.0, 1.5 with m = 5, δ = 0 is presented in Figure (a,b). From this figure, we can confirm that the parameter γ controls the balance between the reaction and diffusion processes. Figure (a,b) also displays the behavior of the numerical solutions with various M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 7 of 11 Fig. 1: The approximate solution θ1(t), θ2(t) via various values of m. Fig. 2: The solution θ1(t), θ2(t) by MSRM and RK4M. values of δ = 0.0, 0.5, 1.0, at m = 5, γ = 0.75. From this figure, we can confirm that the system can exhibit sustained oscillations as δ increases beyond this critical value. Through the above results (Figures 1-4), we can see how the numerical solution behaves when the proposed approach is applied to the values of m, γ, δ. From Figure 3, we can confirm that as the value of parameter γ increases, the value of θ1(t) decreases, and the value of θ2(t) increases, and this is consistent with the nature of the effect of this parameter on chemical reactions of this type. From Figure 4, we find that as the value of δ increases, the value of θ1(t) increases, and the value of θ2(t) decreases, and this is consistent with the nature of the effect of this parameter on chemical reactions of this type. In addition, to validate our numerical solutions, we present a comparison of the relative error (RE) in Table 1 with those solutions obtained by the variational iteration method (VIM) [22] with integer derivative α = 1. We considered the system (1)-(2) in the range [0, 2] with initial conditions θ01 = θ02 = 1, δ = 0.25, γ = 1, and m = 5. Here, to evaluate the RE, we took the numerical solutions obtained by the RK4 method as reference solutions. This is due to the lack of knowledge of the exact solutions of the proposed model. This comparison shows the thoroughness of the proposed method and has the largest convergence interval compared to the VIM. Finally, we can confirm that the proposed method in this article addressed and overcame the shortcomings of the VIM method. M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 8 of 11 Fig. 3: The approximate solution θ1(t), θ2(t) via various values of γ. Fig. 4: The approximate solution θ1(t), θ2(t) via various values of δ. Table 1. The relative error for numerical solutions by using MSRM and VIM [22] at m = 5. Relative error-θ1(t) at m = 5 Relative error-θ2(t) at m = 5 t MSRM VIM MSRM VIM 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 5.15975E-06 2.75961E-06 0.18523E-05 7.13289E-06 6.28410E-06 2.63985E-06 8.02514E-05 7.01478E-05 6.52041E-05 2.15671E-04 6.32587E-04 3.25874E-04 2.32104E-03 4.01573E-03 7.25890E-02 0.21546E-02 0.21456E-01 0.25643E-00 8.00145E-06 3.25874E-06 6.65402E-05 0.01478E-06 2.01478E-07 3.25874E-06 0.25641E-05 6.02587E-06 0.21478E-05 2.31456E-05 8.02145E-04 2.01597E-04 0.23145E-03 1.25489E-03 2.32014E-02 0.02546E-02 0.56479E-01 0.36547E-00 6. Conclusions Through this work, the multi-domain approach, known as the MSRM, is modified to solve the given Brusselator system. The numerical simulation of the model under study M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 9 of 11 was investigated by calculating the relative error with various values of the approximation orderm. By including more terms from the approximation solution series, or by raisingm, we may also regulate the precision of the error and lower it. By assessing the relative error, we can ensure that the proposed technique is accurate and efficient. The proposed method shows its advantage in that calculations are relatively easy to follow and understand. The results obtained by the MSRM are accurate, valid for a longer period, and highly compatible with those obtained by the RK4M and the VIM (but in small intervals). Hence, this technique can be used to solve many other nonlinear problems that appear in the applied sciences and engineering applications. The outcomes additionally demonstrate the proposed MSRM’s accuracy, computational effectiveness, and dependability as a solution for solving this model. All the numerical results were obtained with the help of the Mathematica software program. As future work and an attempt to generalize the results obtained in this paper, we will try to treat one or all of the following suggestions: (i) We deal with the same problem, but using a different type of polynomial. (ii) We consider the model in its fractional form and solve it numerically with the help of a suitable numerical/approximate method. (iii) We try to modify the proposed method in the current work to treat some weaknesses of the method, such as the low accuracy of the resulting solutions compared to other methods, and then use it for solving numerically some other models. Competing interests: There is no conflict of interest. Authors’ contributions: This study was written in collaboration the authors. Acknowledgements The researchers wish to extend their sincere gratitude to the Deanship of Scientific Research at the Islamic University of Madinah for the support provided to the Post- Publishing Program. References [1] M. M. Khader. On the numerical solutions for the fractional diffusion equation. Com- munications in Nonlinear Science and Numerical Simulations, 16:2535–2542, 2011. [2] M. M. Khader and M. Adel. Chebyshev wavelet procedure for solving FLDEs. Acta Applicandae Mathematicae, 158(1):1–10, 2018. [3] O. Abdulaziz, N. F. M. Noor, I. Hashim, and M. S. M. Noorani. Further accuracy tests on the Adomian decomposition method for chaotic systems. Chaos Solitons and Fractals, 36:1405–1411, 2008. [4] B. Batiha, M. S. M. Noorani, I. Hashim, and E. S. Ismail. The multistage variational iteration method for a class of nonlinear systems of ODEs. Physica Scripta, 76:388– 392, 2007. M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 10 of 11 [5] M. S. H. Chowdhury, I. Hashim, and S. Momani. The multistage homotopy- perturbation method: a powerful scheme for handling the Lorenz system. Chaos Solitons and Fractals, 40:1929–1937, 2009. [6] S. S. Motsa, P. Dlamini, and M. Khumalo. A new multistage spectral relaxation method for solving chaotic initial value systems. Nonlinear Dynamics, 72:265–283, 2013. [7] S. S. Motsa, P. G. Dlamini, and M. Khumalo. Solving hyperchaotic systems using the spectral relaxation method. Abstract and Applied Analysis, pages 1–18, 2012. [8] M. S. Khan and M. I. Khan. A novel numerical algorithm based on Galerkin-Petrov time-discretization method for solving chaotic nonlinear dynamical systems. Nonlin- ear Dynamics, 91(3):1555–1569, 2018. [9] S. S. Motsa. A new piecewise-quasilinearization method for solving chaotic systems of initial value problems. Central European Journal of Physics, 10(4):936–946, 2012. [10] M. Karimi and H. S. Nik. A piecewise spectral method for solving the chaotic control problems of the hyperchaotic finance system. International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, 31(3):1–14, 2018. [11] D. Mathale, P. G. Dlamini, and M. Khumalo. Compact finite difference relaxation method for chaotic and hyperchaotic initial value systems. Computational and Applied Mathematics, 37(4):5187–5202, 2018. [12] M. M. Khader. Numerical solutions for the problem of the boundary layer flow of a Powell-Eyring fluid over an exponential sheet using the spectral relaxation method. Indian Journal of Physics, 94(9):1369–1374, 2020. [13] T. Omara. Penalized estimators for modified Log-Bilal regression: simulations and applications. Statistics, Optimization & Information Computing, 14(3):1566–1583, 2025. [14] T. Omara. Robust Lasso Estimator for the Liu-Type regression model and its appli- cations. Statistics, Optimization & Information Computing, 13(5):1819–1831, 2025. [15] N. H. Sweilem, M. M. Khader, and M. Adel. Numerical simulation of fractional Cable equation of spiny neuronal dendrites. Journal of Advanced Research, 5(2):253–259, 2014. [16] M. M. Khader, N. H. Sweilam, and A. M. S. Mahdy. Two computational algorithms for the numerical solution for systems of fractional differential equations. Arab Journal of Mathematical Sciences, 21(1):39–52, 2015. [17] K. S. Nisar, R. Jagatheeshwari, C. Ravichandran, and P. Veeresha. An effective ana- lytical method for the fractional Brusselator reaction-diffusion system. Mathematical Methods in the Applied Sciences, 46(18):18749–18758, 2023. [18] M. M. Khader, J. E. M. Dı́az, K. M. Saad, and W. M. Hamanah. Vieta-Lucas poly- nomials for the Brusselator system with the Rabotnov fractional-exponential kernel fractional derivative. Symmetry, 15(9):1–10, 2023. [19] M. M. Khader and K. M. Saad. A numerical study using Chebyshev collocation method for a problem of biological invasion: fractional Fisher equation. International Journal of Biomathematics, 11(8):1–15, 2018. [20] K. M. Saad, M. M. Khader, J. F. Gomez-Aguilar, and D. Baleanu. Numerical solutions M. Adel et al. / Eur. J. Pure Appl. Math, 18 (4) (2025), 7064 11 of 11 of the fractional Fisher’s type equations with Atangana-Baleanu fractional derivative by using spectral collocation methods. Chaos, 29:1–5, 2019. [21] V. Gafiychuk and B. Datsko. Stability analysis and limit cycle in a fractional system with Brusselator nonlinearities. Physics Letters A, 372(29):4902–4904, 2008. [22] H. Jafari, A. Kadem, and D. Baleanu. Variational iteration method for a fractional- order Brusselator system. Abstract and Applied Analysis, page 496323, 2014. [23] P. Reiterer, C. Lainscsek, Sch, and J. Maquet. A nine-dimensional Lorenz system to study high-dimensional chaos. Journal of Physics A: Math. Gen., 31:7121–7139, 1998. [24] L. N. Trefethen. Spectral Methods in MATLAB. SIAM, Philadelphia, Pa, USA, 2000. [25] J. N. Kouagou, P. G. Dlamini, and S. M. Simelane. On the multi-domain compact finite difference relaxation method for high dimensional chaos: The nine-dimensional Lorenz system. Alexandria Engineering Journal, 59(4):2617–2625, 2020. [26] M. Adel, M. M. Khader, T. A. Assiri, and W. Kallel. Simulating COVID-19 model research using a multidomain spectral relaxation technique. Symmetry, 15:931, 2023. [27] P. Dlamini and S. Simelane. An efficient spectral method-based algorithm for solv- ing a high-dimensional chaotic Lorenz system. Journal of Applied Computational Mechanics, 7(1):225–234, 2021.