EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 2, 2024, 753-771 ISSN 1307-5543 – ejpam.com Published by New York Business Global Computational treatment for the coupled system of viscous Burger’s equations through non-central formula in the method of lines Rawan Alharbi1, Abeer Alshareef2, Huda O. Bakodah1,∗, Aisha Alshaery1 1 Department of Mathematics and statistics, Faculty of Science, University of Jeddah, P.O. Box 80327, Jeddah, Saudi Arabia 2 Department of Basic Sciences, College of Science and Theoretical Studies, Saudi Electronic University, Riyadh 11673, Saudi Arabia Abstract. Viscous Burger’s equation is one of the most celebrated models with immense applica- tions cutting across all aspects of mathematical physics. Thus, the present communication makes use of the non-central formula infused in the method of lines coupled with Runge-Kutta spatial discretization to computationally and efficiently treat the coupled system of viscous Burger’s equa- tions. Further, we numerically tested the derived schemes of the governing model amidst suitable initial and boundary data. In fact, we eventually analyzed the effectiveness of the method via L2 and L∞ norms on some test problems and found it to be robust; indeed, a comparison of the current scheme with some notable approaches in the literature has been established, which are realized to be in perfect agreement. 2020 Mathematics Subject Classifications: 35A25, 35A99, 35C08, 65G99, 65L06, 65M20 Key Words and Phrases: Non-central formula, method of lines, Nonlinear coupled system, viscous Burger’s equations 1. Introduction Viscous Burger’s equation is among the most celebrated models with immense rele- vance in many science and engineering processes, including fluid dynamics, optics, wave scattering, shocks, and shallow water waves to list a few. In particular, a coupled variant of Burger’s equation was devised as the coupled viscous Burger’s equation, which was utilized in modeling the transmission of poly-dispersion in sedimentation processes; read Esipov [8]. More so, in line with the abundant applications of the coupled model, several researchers have in the past and present times deployed dissimilar computational ap- proaches in search of an optimal numerical solution for the governing model. For instance, ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v17i2.5102 Email addresses: rawan mmh@hotmail.com (Rawan Alharbi), as.alshareef@seu.edu.sa (A. Alshareef) hobakodah@uj.edu.sa (H. O. Bakodah), aaal-shaery@uj.edu.sa (A. A. Alshaery) https://www.ejpam.com 753 © 2024 EJPAM All rights reserved. H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 754 Mubaraki et al [1] examined the special class of a generalized reaction-advection-diffusion dynamical model that is called the system of coupled Burger’s equations by employing two promising analytical integration schemes. Also, Rashid and Ismail [2] utilized the Fourier pseudo-spectral analysis method to effectively derive a universal computational scheme for treating the coupled one-dimensional Burger’s equation. El-Shiekh and Al-Nowehy [20] applied the symmetry group method to obtain similarity reductions for generalized symmetric coupled Burgers-type equations. Abazari and Borhanifar [16] used the differ- ential transformation technique on the class of Burger’s and coupled Burger’s equations to acquire some interesting semi-numerical results. Moreover, one can equally find the appli- cation of finite difference numerical method on the acquisition of optimal computational solution for the coupled (2+1)-dimensional Burger’s equation in Srivaslava et. al [22–24]. In addition, other numerical methods that were used on the coupled Burger’s equation include the higher-order trigonometric B-spline [6], semi-Lagrangian method [21], quartic B-spline approach [27], and the quintic B-spline coupled with adaptive time integrator method [26] to mention a few. To see more related work, we refer the reader to see [10, 13, 28]. However, the key aim of the present study is to come up with an effective numerical scheme for the complete treatment of the coupled system of Burger’s equations in mathe- matical physics. In the proposed method, we improve the method of lines (MOL) by using a non-central formula instead of the classical schemes of the first and second derivatives to get a more accurate result. This technique focuses on converting a partial differential equation (PDE) into an ordinary differential equation (ODE), since several methods can used to solve the ODE. In fact, the non-central formula in the method of lines (MOL) with an infusion of time-discretization via the Runge-Kutta numerical method will be deployed for the derivation of the scheme. Besides, the classical MOL has been used by several authors to find efficient numerical solutions for different partial differential equations, in- cluding for instance, the Korteweg-de Vries (KdV) equations [25], the one-dimensional wave propagation problem [9], the KdV-Burger’s equation [17] and the system of coupled elastic beam equations [15], to mention but just a few. Additionally, Sharaf and Bakodah [4] introduced a reliable spatial discretization in MOL for the computational treatment of a class of partial differential equations. Further, other related studies have made use of the same approach, given as in [4], alongside the non-central formula to successfully tackle the regularized-long wave [14] and the equal-width wave [12] equations, respectively; read also the good work of Bakodah [11], which made use of the 7-point formula infused in the spatial-discretization via MOL to solve various forms of Burgers’ equations. More precisely, the current study would devise an appropriate spatial-discretization in the MOL for the coupled system of viscous Burger’s equations. In addition, this approach sets to convert the coupled model endowed with sufficient auxiliary conditions to a coupled system of ordinary differential equations (ODEs). What is more, the unfailing fourth-order Runge-Kutta numerical method will further be sought to numerically treat the resulting ODEs, as against the most widely used finite difference method, which was proven to be superior in the literature; read the computational efficiency of the fourth-order Runge- Kutta numerical method that was reported in [11, 12, 14] and some of the references H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 755 therewith. Lastly, the manuscript is arranged as follows: Section 2 features the central coupled model, Section 3 gives the delineation of the MOL, and Section 4 gives the resulting numerical experiments and discussion; while Section 5 gives some finishing notes. 2. Coupled system for viscous Burger’s equation The coupled nonlinear system of viscous Burger’s equation, which is extensively used in modeling turbulent fluids, theory of shock waves, and stochastic processes among others, is a nonlinear evolution equation, or precisely a nonlinear partial differential equation (PDE) that reads as follows [7, 8] ut = uxx − ηuux − γ(vu)x, vt = vxx − ηvvx − ξ(vu)x, c < x < d, t > 0, (1) where γ, ξ and η are real constants, η is an arbitrary constant depending on the system parameters, more so, the first two constants γ, ξ depend on the parameters of the system, like the Brownian diffusivity, Stokes velocity, and Peclet number, to mention a few. In addition, the celebrating coupled model in Eq.(1) admits broad relevance in various areas of mathematical physics, including, the nonlinear acoustic, fluid mechanics, traffic flow, and gas dynamics to mention but just a few. Furthermore, for a successful deployment and further implementation of the aiming computational approach, it is imperative to prescribe appropriate initial data. Thus, we make consideration of the initial data as follows u(x, 0) = f(x), v(x, 0) = g(x), c ≤ x ≤ d. Once more, we restate here that we would devise an appropriate spatial discretization in the MOL to solve the above initial-value problem for the coupled system of viscous Burgers’ equations. Above and beyond, other recent methods applied to the model can be found in [16] that used the differential transformation method, and [19] that utilized finite-difference scheme among others. 3. Methodology This section delineates the methodology of interest; the MOL coupled with an infusion of a reliable spatial discretization via finite difference method for the computational treat- ment of a class of PDE, the coupled system of viscous Burger’s equations. Specifically, the present section starts off the methodology by converting the coupled model, which is endowed with sufficient auxiliary conditions to a coupled system of ODEs; in addition to the deployment of the consistent fourth-order Runge-Kutta technique for the numerical solution of the resulting ODEs. Accordingly, to begin with, we portray the spatial dis- cretization by considering the spatial variable x, and discretize it into N + 1 uniformly H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 756 spaced grid points as follows xi = xi−1 +∆x, i = 1, 2, 3, ..., N, where ∆x = xi − xi−1 = 1/N is a constant spacing, x0 and xN are the two end-points (specifically, the boundary points), while x′is for i = 2, 3, 4, ..., (N − 1) are the interior points. Next, upon making use of the non-central difference 5-point scheme [4], the con- cerning first- and second-order derivatives are thus approximated as follows dui dx = 1 4!(∆x) [ui−2 − 8ui−1 + 8ui+1 − ui+2] , d2ui dx2 = 1 4!(∆x)2 [−ui−2 + 16ui−1 − 30ui + 16ui+1 − ui+2] . (2) dvi dx = 1 4!(∆x) [vi−2 − 8vi−1 + 8vi+1 − vi+2] , d2vi dx2 = 1 4!(∆x)2 [−vi−2 + 16vi−1 − 30vi + 16vi+1 − vi+2] . (3) Additionally, on using the approximations expressed above in Eq. (1), where uux = 1 2(u 2)x, and vvx = 1 2(v 2)x, the following coupled system of first-order ODEs is thus revealed H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 757 U̇1 = 1 4(∆x)2 [35U0 − 104U1 + 114U2 − 56U3 + 11U4] − η 2(4!)(∆x) [ −25U2 0 + 48U2 1 − 36U2 2 + 16U2 3 − 3U2 4 ] − γ 4!(∆x) (U1 [−25V0 + 48V1 − 36V2 + 16V3 − 3V4]− V1 [−25U0 + 48U1 − 36U2 + 16U3 − 3U4]) , U̇2 = 1 4(∆x)2 [11U1 − 20U2 + 6U3 + 4U4 − 1U5] − η 2(4!)(∆x) [ −3U2 1 − 10U2 2 + 18U2 3 − 6U2 4 + 1U2 5 ] − γ 4!(∆x) (U2 [−3V1 − 10V2 + 18V3 − 6V4 + 1V5]− V2 [−3U1 − 10U2 + 18U3 − 6U4 + 1U5]) , U̇i = 1 4!(∆x)2 [−Ui−2 + 16Ui−1 − 30Ui + 16Ui+1 − Ui+1] − η 2(4!)(∆x) [ U2 i−2 − 8U2 i−1 + 8U2 i+1 − U2 i+2 ] − γ 4!(∆x) (Ui [Vi−2 − 8Vi−1 + 8Vi+1 − Vi+2]− Vi [Ui−2 − 8Ui−1 + 8Ui+1 − Ui+2]) , U0 = U(x0), Ui = U(xi), for i = 3(1)N − 3, U̇N−2 = 1 4!(∆x)2 [−UN−5 + 4UN−4 + 6UN−3 − 20UN−2 + 11UN−1] − η 2(4!)(∆x) [ −U2 N−5 + 6U2 N−4 − 18U2 N−3 + 10U2 N−2 + 3U2 N−1 ] − γ 4!(∆x) (UN−2 [−VN−5 + 6VN−4 − 18VN−3 + 10VN−2 + 3VN−1] − VN−2 [−UN−5 + 6UN−4 − 18UN−3 + 10UN−2 + 3UN−1]), U̇N−1 = 1 4!(∆x)2 [11UN−4 − 56UN−3 + 114UN−2 − 104UN−1 + 35UN ] − η 2(4!)(∆x) [ 3U2 N−4 − 16U2 N−3 + 36U2 N−2 − 48U2 N−1 + 25U2 N ] − γ 4!(∆x) (UN−1 [3VN−4 − 16VN−3 + 36VN−2 − 48VN−1 + 25VN ] − VN−1 [3UN−4 − 16UN−3 + 36UN−2 − 48UN−1 + 25UN ]), (4) H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 758 V̇1 = 1 4(∆x)2 [35V0 − 104V1 + 114V2 − 56V3 + 11V4] − η 2(4!)(∆x) [ −25V 2 0 + 48V 2 1 − 36V 2 2 + 16V 2 3 − 3V 2 4 ] − ξ 4!(∆x) (U1 [−25V0 + 48V1 − 36V2 + 16V3 − 3V4]− V1 [−25U0 + 48U1 − 36U2 + 16U3 − 3U4]) , V̇2 = 1 4(∆x)2 [11V1 − 20V2 + 6V3 + 4V4 − 1V5] − η 2(4!)(∆x) [ −3V 2 1 − 10V 2 2 + 18V 2 3 − 6V 2 4 + 1V 2 5 ] − ξ 4!(∆x) (U2 [−3V1 − 10V2 + 18V3 − 6V4 + 1V5]− V2 [−3U1 − 10U2 + 18U3 − 6U4 + 1U5]) , V̇i = 1 4!(∆x)2 [−Vi−2 + 16Vi−1 − 30Vi + 16Vi+1 − Vi+1] − η 2(4!)(∆x) [ V 2 i−2 − 8V 2 i−1 + 8V 2 i+1 − V 2 i+2 ] − ξ 4!(∆x) (Ui [Vi−2 − 8Vi−1 + 8Vi+1 − Vi+2]− Vi [Ui−2 − 8Ui−1 + 8Ui+1 − Ui+2]) , V0 = V (x0), Vi = V (xi), for i = 3(1)N − 3, V̇N−2 = 1 4!(∆x)2 [−VN−5 + 4VN−4 + 6VN−3 − 20VN−2 + 11VN−1] − η 2(4!)(∆x) [ −V 2 N−5 + 6V 2 N−4 − 18V 2 N−3 + 10V 2 N−2 + 3V 2 N−1 ] − ξ 4!(∆x) (UN−2 [−VN−5 + 6VN−4 − 18VN−3 + 10VN−2 + 3VN−1] − VN−2 [−UN−5 + 6UN−4 − 18UN−3 + 10UN−2 + 3UN−1]), V̇N−1 = 1 4!(∆x)2 [11VN−4 − 56VN−3 + 114VN−2 − 104VN−1 + 35VN ] − η 2(4!)(∆x) [ 3V 2 N−4 − 16V 2 N−3 + 36V 2 N−2 − 48V 2 N−1 + 25V 2 N ] − ξ 4!(∆x) (UN−1 [3VN−4 − 16VN−3 + 36VN−2 − 48VN−1 + 25VN ] − VN−1 [3UN−4 − 16UN−3 + 36UN−2 − 48UN−1 + 25UN ]). (5) Therefore, to numerically solve the coupled system of ODEs outlined in Eqs. (4)-(5), we then deploy the reliable fourth-order Runge-Kutta method. In fact, we re-express these equations in a more compact form as follows U̇i = F (Ui, Vi), i = 1, 2, ..., N − 1, (6) H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 759 V̇i = F (Ui, Vi), i = 1, 2, ..., N − 1. (7) Moreover, the implementation of MOL begins with the replacement of the involving spatial derivatives, which arise from the governing coupled Burger’s model with discretized alge- braic approximations. In light of this, explicit representation of the spatial derivatives is not realistic. Hence, a coupled system of ODEs is revealed and further solved numerically, which then gives the approximation of the novel PDE. 3.1. MOL via the fourth-order Runge-Kutta method Here, the approach starts by deploying the fourth-order Runge-Kutta method to solve the resulting ODEs [19], in Eqs. (4)-(5), which were revealed independently with the help of spatial discretization of the finite difference approach. Certainly, in trying to incorporate the application of MOL, it is pertinent to recall that highly accurate computational results are achieved when deploying multistep approaches, as against the use of one-step methods, where the intermediate values of the solutions and the related derivatives are produced and utilized at each computational stage. At first, the following definition Un = [ Un 1 , U n 2 , ..., U n N−1 ]t , and V n = [ V n 1 , V n 2 , ..., V n N−1 ]t , is considered; then followed by the definition of the Runge-Kutta method in the following pattern Un+1 = Un + 1 6 [K11 + 2K12 + 2K13 +K14] , (8) V n+1 = V n + 1 6 [K21 + 2K22 + 2K23 +K24] , (9) where, K11 = ∆tF (Un, V n) , (10) K21 = ∆tG (Un, V n) , (11) K12 = ∆tF ( Un + K11 2 , V n + K21 2 ) , (12) K22 = ∆tG ( Un + K11 2 , V n + K21 2 ) , (13) K13 = ∆tF ( Un + K12 2 , V n + K22 2 ) , (14) K23 = ∆tG ( Un + K12 2 , V n + K22 2 ) , (15) K14 = ∆tF (Un +K13, V n +K23) , (16) K24 = ∆tG (Un +K13, V n +K23) , (17) H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 760 with ∆t representing the temporal step size. In addition, and without further delay, the application of the Runge-Kutta method requires the evaluation of eight vector functions at every computational time step; thus, we only report the resulting numerical results in the subsequent section to suppress the length of the paper. 4. Numerical experiments and discussion The present section considers certain test models to portray the compliance and compe- tence of the devised scheme on the initial-value problem of the coupled system of Burger’s equations. However, to assess the accuracy of the proposed numerical results, we resort to the engagement of the norm two L2 and norm infinity L∞ error estimators, which take the following explicit definitions L2(u) = √√√√[ h N−1∑ i=1 |uni − Un i | 2 ] , (18) L∞(u) = max 1≤i≤N−1 |uni − Un i | . (19) In addition, we are set to establish a comparative assessment between the obtained/proposed computational results, via the devised approach with the available exact solutions in the literature. Moreover, we will be assessing the proposed scheme with other numerical schemes that were once proposed for the governing model. Lastly, all the computational simulations of the study will be performed with the aid of Maple 17. Example 1. Let us make consideration of the governing coupled system of Burger’s equa- tion (1) with η = 2, while generalized values for γ and ξ are assumed. In fact, the exact hyperbolic solution of the coupled equation is reported by Soliman [3] as follows u(x, t) = D (1− tanh(B(x− 2Bt))) , v(x, t) = D ( 2ξ − 1 2γ − 1 − tanh(B(x− 2Bt)) ) , (20) with D = 0.05 and B = D(4γξ − 1)/(4γ − 2). Moreover, sufficient initial (at t = 0) and boundary (at x = c and x = d) data follow accordingly from the above exact solution. Hence, we have established in Tables 1 and 2, a comparative analysis between the proposed solution via the devised approach, and that obtained by the Chebyshev spectral collocation method [5]. Moreover, the earlier stated error norm estimators are deployed for the present assessment, which are, L2 and L∞; while fixing D = 0.05 at two different time levels, that is, t = 0.5 and t = 1.0, respectively. Indeed, it can be noted from these tables that the obtained results through the proposed approach outperformed the challenging methodology. Additionally, the numerical solution obtained by 5-point non-central MOL and the exact solution of u(x, t) and v(x, t) for a fixed time t = 0.5 and for various values of h are H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 761 γ ξ t L2 L∞ Present Ref [5] 0.1 0.3 0.5 6.35232 ×10−6 4.49177 ×10−6 4.16×10−5 1 1.019395×10−5 7.20821×10−6 8.23×10−5 0.3 0.03 0.5 6.38257 ×10−6 4.51316×10−6 4.59×10−5 1 4.105386×10−5 2.902946×10−5 9.16×10−5 Table 1: Comparison of error norms for the solution pair u(x, t) when N = 10, c = −10, d = 10, ∆t = 0.1, and D = 0.05 of Example 1. γ ξ t L2 L∞ Present Ref [5] 0.1 0.3 0.5 1.49091 ×10−6 1.05423×10−6 4.99×10−5 1 1.019395×10−5 7.20821×10−6 9.92×10−5 0.3 0.03 0.5 3.513746 ×10−5 2.484593 ×10−5 1.81×10−4 1 4.105384 ×10−5 2.902945×10−5 3.62×10−4 Table 2: Comparison of error norms for the solution pair v(x, t) when N = 10, c = −10, d = 10, ∆t = 0.1, and D = 0.05 of Example 1. given in Tables 3 and 4. Moreover, Tables 5 and 6 further made use of L∞ as an error norm estimator to establish a comparative study between the present computational results and those obtained using various approaches, including the Fourier pseudo-spectral method (FPM) [2], Chebyshev spectral collocation (CSC) [5], and the cubic B-spline collocation (CBC) [18]. Indeed, among all the comparing methods, the currently devised approach happens to proffer better results with high accuracy in contrast with the aforementioned approaches in [2, 5, 18]. In addition, we have shown in figures 1 and 2 the graphical display of the obtained computational and exact solution fields u(x, t) and v(x, t), while fixing the following values: N = 10, ∆t = 0.1, t = 0.5, γ = 1, ξ = 2, and D = 0.1. ∆t t h u(x, t) L∞present Exact 0.1 0.5 2 0.0527592618 0.0527547700 4.4918×10−6 1 0.0581890811 0.0581832978 5.7833×10−6 0.6666666667 0.0599622104 0.0599560894 6.1209×10−6 0.5714285714 0.0604642722 0.0604580646 6.2076×10−6 0.5 0.0608393804 0.0608331108 6.2696×10−6 0.3921568627 0.0614032873 0.0613969292 6.3581×10−6 Table 3: Comparison between MOL and exact solutions for the solution pair u(x, t) when c = −10, d = 10, γ = 0.1, ξ = 0.3, and D = 0.05 using various mesh sizes of Example 1. H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 762 ∆t t h v(x, t) L∞present Exact 0.1 0.5 2 0.02775582427 0.0277547700 1.0542×10−6 1 0.0331857077 0.0331832978 2.4100×10−6 0.6666666667 0.034958886 0.0349560894 2.7965×10−6 0.5714285714 0.0354609639 0.0354580646 2.8993×10−6 0.5 0.0358360849 0.0358331108 2.9742×10−6 0.3921568627 0.0364000124 0.0363969291 3.0833×10−6 Table 4: Comparison between MOL and exact solutions for the solution pair v(x, t) when c = −10, d = 10, γ = 0.1, ξ = 0.3, and D = 0.05 using various mesh sizes of Example 1. Method t γ ξ L∞ Present (RK4) 0.5 0.1 0.3 4.49177 ×10−6 0.3 0.03 4.51316×10−6 1 0.1 0.3 7.20821×10−6 0.3 0.03 2.902946×10−5 FPM [2] 0.5 0.1 0.3 9.619×10−4 0.3 0.03 4.310×10−4 1 0.1 0.3 1.153×10−3 0.3 0.03 1.268×10−3 CSC [5] 0.5 0.1 0.3 4.38×10−5 0.3 0.03 4.58×10−5 1 0.1 0.3 8.66×10−5 0.3 0.03 9.16×10−5 CBC [18] 0.5 0.1 0.3 4.16×10−5 0.3 0.03 4.59×10−5 1 0.1 0.3 8.25×10−5 0.3 0.03 9.18×10−5 Table 5: Comparison of norm L∞ errors for various approaches at different times for the solution pair u(x, t) of Example 1. H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 763 Method t γ ξ L∞ Present (RK4) 0.5 0.1 0.3 1.05423×10−6 0.3 0.03 2.48459 ×10−5 1 0.1 0.3 7.2082×10−6 0.3 0.03 2.90294×10−5 FPM [2] 0.5 0.1 0.3 3.33×10−4 0.3 0.03 1.14×10−3 1 0.1 0.3 1.16×10−3 0.3 0.03 1.63×10−3 CSC [5] 0.5 0.1 0.3 4.99×10−5 0.3 0.03 1.81×10−4 1 0.1 0.3 9.92×10−5 0.3 0.03 3.62×10−4 CBC [18] 0.5 0.1 0.3 1.48×10−4 0.3 0.03 5.72×10−4 1 0.1 0.3 4.77×10−5 0.3 0.03 3.61×10−4 Table 6: Comparison of norm L∞ errors for various approaches at different times for the solution pair v(x, t) of Example 1. Figure 1: Numerical and exact solutions for u(x, t) of Example 1 when N = 10, ∆t = 0.1, γ = 1, ξ = 2, and D = 0.1. H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 764 Figure 2: Numerical and exact solutions for v(x, t) of Example 1 when N = 10, ∆t = 0.1, γ = 1, ξ = 2, and D = 0.1. H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 765 • We compare use 3-point formula, 5-point formula, and 7-point formula results. (i) The effect of increasing points formulas . ∆t = 0.01, t = 0.1, N = 50 Z-point formula u(x, t) L2 L∞ 3-point formula 1.25330 ×10−6 1.98163 ×10−6 5-point formula 3.7714 ×10−7 5.9630×10−7 7-point formula 1.21805 ×10−6 1.92591 ×10−6 Table 7: Influence of increasing points formulas for the solution pair u(x, t) of Example 1 when γ = 0.1, ξ = 0.3. ∆t = 0.01, t = 0.1, N = 50 Z-point formula v(x, t) L2 L∞ 3-point formula 8.2861×10−7 1.31015×10−6 5-point formula 1.6555 ×10−7 2.6176×10−7 7-point formula 7.9487 ×10−7 1.25680×10−6 Table 8: Influence of increasing points formulas for the solution pair v(x, t) of Example 1 when γ = 0.1, ξ = 0.3. Finally, we conclude that the results obtained by using the proposed method are better than other methods introduced in Tables 1-6 at different values of time. Furthermore, from Table 7 and Table 8, we observe that the accuracy of the numerical solutions at each field u(x, t), v(x, t) with the 5-point formula is more efficient. Example 2. Let us again make consideration of the governing coupled system of Burgers’ equation (1) when η = −2, and γ = ξ = 5 2 . In addition, the exact solution of the coupled model for 0 ≤ x ≤ 1 was analytically reported as follows [16] u(x, t) = v(x, t) = µ [1− tanh (1.5µ(40(x− 0.5)− 3µt))] (21) where µ is a constant, which is indeed arbitrary. More so, the related initial and boundary data follow accordingly from the above analytical solution. Accordingly, we have reported in Table 9 the acquired computational results for different fixed values of µ. In addition, we have shown in figures (3) and (4) the comparison between the acquired computational solution through MOL and the analytical solution, on the other hand, using the following data µ = 0.05, 0.1 and ∆t = 0.005 at t = 0.05. Notably, it can be observed from Table 9 that a decrease in the value of µ increases the accuracy. Besides, it is also noted that the error gets smaller as time grows. H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 766 t u(x, t) = v(x, t) Absolute Error Numerical results Exact µ = 0.05 0.005 0.0693850254 0.0693855034 5.7161 ×10−5 0.01 0.0129083018 0.0129085548 5.6108 ×10−5 µ = 0.1 0.005 0.1674079059 0.1674103611 8.1566 ×10−4 0.01 0.0042990998 0.0042990998 1.3146×10−4 Table 9: Numerical solution with the exact solution of Example 2 when γ = ξ = 5 2 , µ = 0.05, 0.10, N = 11, c = 0, d = 1, and ∆t = 0.001. Figure 3: The behavior of the numerical and analytical solution for the 5-point formula of Example 2 for N = 11, ∆t = 0.005, t = 0.005, γ = ξ = 5 2 . H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 767 Figure 4: The behavior of the numerical and analytical solution for the 5-point formula of Example 2 for N = 11, ∆t = 0.005, t = 0.005, γ = ξ = 5 2 . • 3-point formula (i) The effect of increasing node points N . ∆t = 0.001, t = 0.005 N u(x, t)=v(x, t) Absolute Error Numerical results Analytical results 11 0.0693318373 0.0693855034 5.36661×10−5 21 0.0863297722 0.0864966205 1.668483 ×10−4 31 0.0901130213 0.0902545100 1.414887 ×10−4 51 0.0925043358 0.0926179683 1.136325 ×10−4 61 0.0930218434 0.0931283850 1.065417 ×10−4 71 0.0933731026 0.0934746418 1.015391 ×10−4 Table 10: The effect of increasing node points N of the numerical and analytical solution with µ = 0.05 of Example (2). H. O. Bakodah et al / Eur. J. Pure Appl. Math, 17 (2) (2024), 753-771 768 (ii) The effect of increasing time. ∆t = 0.01, N = 51 t u(x, t)=v(x, t) L2 L∞ 0.010 6.306876×10−5 4.5040103×10−4 0.050 1.5838520×10−4 1.13109658×10−3 0.1 4.7700242 ×10−4 3.40647865×10−3 0.2 1.91607952×10−3 0.01368354479 Table 11: The effect of increasing time when µ = 0.1 of the numerical and analytical solution of Example (2). • We compare by using 3-point formula, 5-point formula, and 7-point formula results. (i) The effect of increasing points formulas . ∆t = 0.001, t = 0.050 Z-point formula u(x, t) = v(x, t) L2 L∞ 3-point formula 6.26752 ×10−6 4.475902×10−5 5-point formula 5.02258 ×10−6 3.586836×10−5 7-point formula 2.89533 ×10−5 2.067682 ×10−4 Table 12: Influence of increasing points formulas of Example (2) when µ = 0.1, N = 51. It can be seen that the accuracy of the numerical solutions with the 5-point formula is more efficient. But when N = 11 and t = 0.009, the 5-point formula will be the best also: ∆t = 0.001, t = 0.009 Z-point formula u(x, t) = v(x, t) L2 L∞ 3-point formula 9.365749 ×10−5 3.1062676 ×10−4 5-point formula 9.192620 ×10−5 3.0488472×10−4 7-point formula 1.3584093 ×10−3 4.5053339×10−3 Table 13: Influence of increasing points formulas of Example (2) when µ = 0.1, N = 11. REFERENCES 769 5. Conclusions In conclusion, a coupled system of one-dimensional viscous Burger’s equation has been computationally tackled through the application of the 5-point non-central formulas in- corporated in MOL. The devised approach was further assessed by considering two test models of interest. Further, the obtained computational results when compared with their exact counterparts were found to be in perfect agreement. Indeed, the values of L2 and L∞ norms that were utilized as error estimators were discovered to be negligible; meaning, the devised scheme was reliable and effective. Based on the results, it is noted that in some cases the 5-point formula effect in significant improvement over the 3- and 7- points formulas. Notably, the acquired computational scheme concerning the proposed scheme was found to be better with the non-central point, as against the central point that gave results with relatively higher errors. As a result, we suggest that the proposed methods should be applied to other types of viscous Burger’s equations. Furthermore, we will address a higher order of viscous Burger’s equations in our future work. Disclosure statement No potential conflict of interest was reported by the authors. References [1] Mubaraki A, Kim H, Nuruddeen R, Akram U, and Akbar Y. Wave solutions and numerical validation for the coupled reaction-advection-diffusion dynamical model in a porous medium. Theor Phys, 74:125002, 2022. [2] Rashid A and Ismail A. A fourier pseudospectral method for solving coupled viscous burgers equations. Comput Meth Appl Math, 9:412–420, 2009. [3] Soliman A. The modified extended tanh-function method for solving burgers-type equations. Physica A: Stat Mech Appl, 361:394–404, 2006. [4] Sharaf AA and Bakodah HO. A good spatial discretisation in the method of lines. Appl Math Comput, 171:1253–1263, 2005. [5] Khater AH, Temsah R, and Hassan M. A chebyshev spectral collocation method for solving burgers’-type equations. J Comput Appl Math, 222:333–350, 2008. [6] Onarcan AT and Hepson OE. Higher order trigonometric b-spline algorithms to the solution of coupled burgers’ equation. AIP Conference Proceedings, 926, 2018. [7] Jin-Bing C, Xian-Guo G, and Zhi-Jun Q. New finite-gap solutions for the coupled burgers equations engendered by the neumann systems. Chinese Phys B, page 090403, 2010. REFERENCES 770 [8] SE Esipov. Coupled burgers equations: a model of polydispersive sedimentation. Phys Rev E, 52:3711, 1995. [9] Shakeri F and Dehghan M. The method of lines for solution of the one-dimensional wave equation subject to an integral conservation condition. Comput Math Appl, 56:2175–2188, 2008. [10] Ge-JiLe H, Javid K, Khan S, Raza M, Khan M, and Qayyum S. Double diffusive convection and hall effect in creeping flow of viscous nanofluid through a convergent microchannel: a biotechnological applications. Comput Methods Biomech Biomed, 24:1326–1343, 2021. [11] Bakodah HO. Non-central 7-point formula in the method of lines for parabolic and burgers’ equations. IJRRAS, 8:328–336, 2011. [12] Bakodah HO and Banaja MA. The method of lines solution of the regularized long- wave equation using runge-kutta time discretization method. Math Probl Eng, 2013:1– 8, 2013. [13] Khan M, Kadry S, Chu Y, and Waqas M. Modeling and numerical analysis of nano- liquid (titanium oxide, graphene oxide) flow viscous fluid with second order velocity slip and entropy generation. Chin J Chem Eng, 31:17–25, 2021. [14] Banaja MA and Bakodah HO. Runge-kutta integration of the equal width wave equation using the method of lines. Math Probl Eng, 2015:1–10, 2015. [15] Sarker P and Chakravarty UK. A generalization of the method of lines for the numer- ical solution of coupled, forced vibration of beams. Math Comp Simul, 170:115–142, 2020. [16] Abazari R and Borhanifar A. Numerical study of the solution of the burgers and coupled burgers equations by a differential transformation method. Comp Math Appl, 59:2711–2722, 2010. [17] El Sadat R and Ali MR. Application of the method of lines for solving the kdv-burger equation. BISKA NTMSCT, 12:39–51, 2017. [18] Mittal RC and Arora G. Numerical solution of the coupled viscous burgers’ equation. Commun Nonlinear Sci Num Simul, 16:1304–1313, 2011. [19] LeVeque RJ. Finite difference methods for ordinary and partial differential equa- tions: steady-state and time-dependent problems. Society for Industrial and Applied Mathematics, 2007. [20] El-Shiekh R.M. and Al-Nowehy AG.A.A.H. Symmetries. Reductions and different types of travelling wave solutions for symmetric coupled burgers equations. Int J Appl Comput Math, 8, 2022. REFERENCES 771 [21] Bak S, Kim P, and Kim D. A semi-lagrangian approach for numerical simulation of coupled burgers’ equations. Commun Nonlinear Sci Num Simul, 69:31–44, 2019. [22] Srivastava VK and Tamsir M. Crank-nicolson semi implicit approach for numerical solutions of two-dimensional coupled nonlinear burgers’ equations. Int J Appl Mech Eng, 7, 2012. [23] Srivastava VK, Tamsir M, Bhardwaj U, and Sanyasiraju Y. Crank-nicolson scheme for numerical solutions of two dimensional coupled burgers’ equations. Int J Sci Eng Res, 2:1–7, 2011. [24] Srivastava VK, Awasthi MK, and Singh S. An implicit logarithmic finite-difference technique for two dimensional coupled viscous burgers’ equation. AIP Adv., 3:122105, 2013. [25] Schiesser WE. Method of lines solution of the korteweg-de vries equation. omp Math Appl, 28:147–154, 1994. [26] Cicek Y, Gucuyenen KN, Bahar E, Gurarslan G, and Tangolu G. A new numerical algorithm based on quintic b-spline and adaptive time integrator for coupled burger’s equation. Comput Meth Diff Equ, 11:130–142, 2023. [27] Zhang Y, Lin J, Reutskiy S, Sun H, and Feng W. The improved backward substitution method for the simulation of time-dependent nonlinear coupled burgers’ equations. Results Phy, 18:103231, 2020. [28] Zhou Y, Li C., and Stynes M. A fast second-order predictor-corrector method for a nonlinear time-fractional benjamin-bona-mahony-burgers equation. Numer Algor, 95:693–720, 2024.