An efficient numerical scheme based on Lucas polynomials for the study of multidimensional Burgers-type equations

We propose a polynomial-based numerical scheme for solving some important nonlinear partial differential equations (PDEs). In the proposed technique, the temporal part is discretized by finite difference method together with θ -weighted scheme. Then, for the approximation of spatial part of unknown function and its spatial derivatives, we use a mixed approach based on Lucas and Fibonacci polynomials. With the help of these approximations, we transform the nonlinear partial differential equation to a system of algebraic equations, which can be easily handled. We test the performance of the method on the generalized Burgers–Huxley and Burgers–Fisher equations, and oneand two-dimensional coupled Burgers equations. To compare the efficiency and accuracy of the proposed scheme, we computed L∞, L2, and root mean square (RMS) error norms. Computations validate that the proposed method produces better results than other numerical methods. We also discussed and confirmed the stability of the technique.


Introduction
Nonlinear partial differential equations (PDEs) are used to model many physical phenomena arising in sciences and engineering. As a result of their considerable applications and popularity, much attention has been devoted to develop an accurate and efficient numerical method for solving PDEs. Consider a one-dimensional nonlinear parabolic partial differential equation of the form where α, β, γ are real parameters, m is a positive integer, and f (Y ) is a nonlinear function. The initial and boundary conditions are Y (a, t) = g 1 (t), Y (b, t) = g 2 (t). ( When f (Y ) = Y (1 -Y m )(Y m -), Eq. (1) defines the generalized Burger-Huxley (GBH) equation Equation (4) describes the interaction between convection, reaction, and diffusion processes. When α = 0, β = 1, and m = 1, Eq. (4) reduces to the Huxley equation investigating wall motion in liquid crystallography and propagation of pulse in nerve fibers. For γ = 0, m = 1, and β = 1, Eq. (4) reduces to the Burger equation used for analysis of nonlinear wave propagation, aspect of turbulence, traffic flows, and shock waves [1]. Similarly, when f (Y ) = Y (1 -Y m ), Eq. (1) becomes the generalized Burger Fisher (GBF) equation The generalized Burger-Fisher equation has a wide application in fluid mechanics, gas dynamics, plasma physics, number theory, elasticity, and heat conduction problems [2]. Equation (5) is a highly nonlinear model, which includes a combination of reaction, convection, and diffusion mechanisms. When γ = 0 and m = 1, Eq. (5) reduces to the Fisher equation having applications in population biology, chemistry, and biological sciences such as spreading of bacterial colonies, spread of reaction fronts in chemically bistable systems, and switching in nonlinear optics [1]. Next, we consider the following two-dimensional coupled viscous Burger equations: Z t + μ{Z ξξ + Z ηη } + ν(ZZ ξ ) + β(YZ) ξ + γ {YZ ξ + ZZ η } = 0, (ξ , η) ∈ , with initial and boundary conditions where and ∂ represent the domain and its boundary, respectively, and μ, η, α, and β are arbitrary constants. The system was introduced by Esipov [3] to study a model of polydispersive sedimentation. This system has numerous applications in science and engineering such as gas dynamics, viscous flow of turbulence, shock waves, sedimentation of particles in fluid suspension, elasticity, and heat conduction problems [4,5]. Many numerical methods have been applied to approximate solutions of these equations. For example, the finite difference method [6], spectral method [7], differential quadrature method [8,9] and Adomian decomposition method [10]. Khattak [11] used a meshfree collocation method, whereas Zhu and Kang [12] applied B-spline interpolation for solution the Burger-Fisher equation. Celik [13] studied the Haar wavelet method for solving GBH equation. Dehghan [14] worked on a mixed collocation and finite difference method. Haq et al. [15] numericallly solved the Burger-Huxley equation using a meshless method of line. Zhang et al. [16] used the local discontinuous Galerkin method, whereas in [17] the author proposed a pseudospectral method for approximation of GBF equation.
Wasim et al. [18] used a hybrid B-spline collocation technique for approximation of GBH and GBF equations. Mittal and Tripathe [2] proposed a cubic B-spline technique.
The modified Burger equation has been investigated in [19,20] using a meshless method and hybrid Haar wavelet finite difference method. The author of [21] obtained approximate solution of coupled Burger equations using Adomian-Pade technique. Khater et al. [5] explained the cubic-spline collocation method for solving the coupled Burger equations. Recently, Mittal and Jiwari [22] obtained an approximate solution of onedimensional coupled Burger equation with the help of the differential quadrature method.
Dehghan et al. [23,24] proposed a mixed finite difference and Galerkin method and multisymplectic box method for numerical study of Burgers equations. Oruc et al. [25,26] applied a unified finite difference Chebyshev wavelet approach for time fractional Burger equations. The same authors studied the Chebyshev wavelet method for approximation of coupled Burgers equations [27]. Ali et al. [4] applied a meshfree collocation method based on the Crank-Nicolson method for time discretization and radial basis function for space discretization to solve two-dimensional coupled Burger equations, whereas the meshless method of radial basis functions (RBFs) and local RBFs were described in [28,29] for approximate solutions of Burger-type equations. In [30] a multiscale variational algorithm was combined with the Kriging element-free Galerkin method to produce the discontinuous solutions of Burgers type equations. Srivastava et al. [31] studied a fully implicit finite difference scheme for solving two-dimensional coupled viscous Burger equations.
In this work, we compute numerical solutions of the generalized Burger-Huxley, Burger-Fisher, and coupled Burger equations using mixed Lucas and Fibonacci polynomials combined with finite differences. The main advantage of the proposed scheme is that the higher-order derivatives can be easily computed using relation of Lucas and Fibonacci polynomials. Moreover, the proposed scheme produces better accuracy for small number of collocation points, which reduces the computational cost. These polynomials have considerable applications in the area of ordinary differential equations. For example, Elhameed and Youssri [32,33] described connection between Chebyshev and Lucas polynomials and obtained accurate solutions of boundary value problems. In [34][35][36] the author implemented the Lucas polynomials for solutions of fractional and coupled fractional differential equations in the Caputo sense. Mostefa [37] proposed the Lucas sequence for approximation of integro-differential equations. Cetin [38] obtained numerical solution of higher-order differential equations using the Lucas polynomial approach.
Farshid et al. [39] proposed a Fibonacci polynomial approach for numerical solution of Volterra-Fredholm integral differential equations. Bayku [40] presented a hybrid Taylor-Lucas polynomial method and obtained numerical solution of delay difference equations.
Oruc [41,42] for the first time applied these polynomials for solution of time-dependent partial differential equations called a mixed Lucas and Fibonacci polynomial technique.
The rest of the paper is organized as follows. In Sect. 2, we discuss description of solution methodology. In Sect. 3, describe the stability of the method. In Sects. 4 and 5, we present numerical experiments followed by conclusion of the paper.

Methodology
In this section, we give a description of the proposed method for two different cases. The suggested technique will be tested by some examples.

Case.1 Nonlinear PDEs
For solution of Eq. (1) using the proposed technique, we discretize the time derivative of the equation with finite differences and apply the θ -weighted scheme to its spatial part to get where Y n = Y (ξ , t n ), t n = nδt, n = 1, 2, . . . , N , δt is the time step. The nonlinear term in Eq. (9) is linearized using the lagging method given as Using Eq. (10) in Eq. (9), we get Approximating Y n (ξ ) by Lucas polynomials is as follows: where C n k are unknown coefficients to be computed, and L k (ξ ) are the Lucas polynomials defined by [41] To determine C n k , we use the collocation method. At collocation points ξ i = a + iδξ , Eq. (12) can be written as Putting the values from Eq. (13) into Eq. (11), we obtain the following system of N equations: The primes in Eq. (14) represent differentiation with respect to ξ , which allows us to replace L k by Fibonacci polynomials [41]: where F k (ξ ) are the Fibonacci polynomials defined as [41].
and D is differentiation matrix given by [41] where d is the square matrix of order N defined by Substituting the values from Eq. (15) into Eq. (14), we can write Similarly, Eqs. (2) and (3) can be transformed to In matrix form, Eqs. (16)- (17) can be written as where H, G, and A are square matrices of order N with components given by Solution of Eq. (18) give required unknowns C, and hence a solution of problem (1) can be obtained with the help of Eq. (13).

Case 2: Coupled PDEs
To construct a scheme for coupled Burger Eqs. (6)- (7), discretizing the temporal and spatial parts in a similar way as discussed in case 1, we have For θ = 1/2, these equations become the well-known Cran-Nicolson scheme with accuracy O(δt 2 ) [43]. For nonlinear terms, we use the lagging method Using Eq. (24) in Eqs. (22)- (23) and waiving the error terms, we get Now we approximate Y n and Z n by Lucas polynomials as follows: where λ n km and C n km are unknown coefficients, and ξ i = η i = a + (i -1) dξ i , dξ = dη, are the regular collocation points or the Chebyshev-Gauss-Lobatto (CGL) collocation points In these equations, "*" represents the usual product. Similarly, boundary conditions given in Eq. (8) take the form considering the relation between the Lucas and Fibonacci polynomials where F(ξ ) and D have the same meaning as before. Putting values from Eq. (32) into Eqs. (28)- (29), the matrix form of Eqs. (28)-(31) can be written as and Equations (33)- (34) can be written as where Since Y n = AC n and Z n = Aλ n , we get If M 1 , M 2 are fully ranked, then M -1 1 and M -1 2 exist [44,45]. In our computation, these matrices are fully ranked, which is checked using Matlab command rank(M i ). Therefore system (35)-(36) can be solved for unknown C n and λ n , and a solution of original coupled equations can be obtained from Eq. (27).

Numerical stability analysis
To check the stability of the proposed technique, we use the matrix method. For this purpose, first rewrite Eq. (18) as follows: If Y denotes the approximate solution, and u denotes the exact one, then the error is defined as Substituting the values from Eq. (39) into Eq. (40), we get Here M = WH -1 GW -1 is known as the implication matrix. Scheme (39) is stable if the matrix M satisfies the Lax-Richtmyer stability condition

Numerical examples
In this section, we apply the proposed technique to some problems. Efficiency of the scheme is tested by comparing the obtained results with exact and approximate solutions available in the literature. We study the accuracy of the method in the form of the root mean square (RMS), L 2 and L ∞ norms, and the time convergence rate given by .
Example 1 Consider Eq. (4) with β = 1 and exact solution [2] Y (ξ , t) = 2 where We obtain the approximate solution the proposed method taking special domains [0, 1] and [-10, 20]. Initial and boundary conditions are taken from the exact solution. In this example, we discuss various cases for different values of the parameters α, γ , , and m appearing in Eq. (4). We compute the solution for both regular and CGL collocation points. We compare the computed results in the form of error norms with those available in the literature [2,7]. From comparison it is clear that the present method gives better accuracy or comparable results with those available in the literature. We can observe from the comparison that the proposed method produces slightly more accurate results for CGL collocations points than for regular points. We also report and show in tables the comparison carried out for stability of the scheme in the form of spectral radius.  Table 1 it is clear that the accuracy of the solution increases with time where the system remains stable, that is, ρ(M) < 1. The computed results are compared with those of the collocation cubic B-spline method [2]. From the comparison it is obvious that the present technique gives better accuracy than the collocation cubic B-spline method.
Case 1.2: In this case, we take the parameters are α = γ = m = 1 and = 0.5 with nodal points N = 10. The solution has been computed at different time levels T = 15, 30, 60, and 120 with step size δt = 0.01. For accuracy of the scheme, different error norms were calculated and compared with the error norms of cubic B-spline [2] and tabulated in Tables 2     Table 2 we observe that for small times, the results of cubic B-spline are better than those of the present method, but as time increases, the accuracy of the proposed technique increases. We can also observe from the table that the spectral radius ρ(M) < 1 for all time levels, which shows the stability of the scheme.  Table 4. From the table it is clear that the present method gives an excellent solution in comparison to available techniques. We can observe from the table that the present method gives the same accuracy for all values of the parameter m, whereas the accuracy of spectral and cubic B-spline methods [2,7] suddenly decreases, which reflects the feasibility of the proposed scheme. The solution is also been computed using regular grid point and shown in Table 6.
Case 1.4: In this case, we take α = γ = m = 1 and = 0.01. The solution is computed over the domain [-10, 20] with step size δt = 0.003. The error norms and spectral radius  are computed for different time levels T = 0.1, 1, 5, and 10, which are presented in Tables 5 and 6. From the table we can observe that the present method gives better results than those available in the literature even for large space and time domains. We can easily understand from these tables that as the domain increases, the spectral radius remains ρ(M) < 1, which shows the stability of the proposed scheme. Overall, it is obvious that the present method gives better results and is flexible to implement.
Example 2 Consider Eq. (5) with β = 1. In this case, the exact solution [2] is given by where w 1 = -αm 2(1 + m) , The  Tables 7 and 9. From Table 7 we can notice that the proposed scheme is stable in the given domain and gives a better accuracy than the collocation cubic B-spline method [2].
Case 2.2: Now we take m = 2 and compute the solution for different time levels T = 5, 10, 15, and 20. The error norms along with spectral radius were obtained at each time level, and the results were compared with existing numerical solutions described in Tables 8 and 9. It is obvious from the table that the accuracy increases with time and converges toward the true solution. It is also clear from the table that in all the three error   norms the proposed method gives excellent results at each time level as compared to available results in the literature. The table also shows that the value of spectral radius remains less than 1, which clarifies the stability of the scheme.
with initial and boundary conditions [2] Y (ξ , 0) = sin(πξ ), 0 ≤ x ≤ 1, The solution is computed for different time levels T = 0.1, 0.3, 0.6, and 0.9 and various values of β = 2 -1 , 2 -4 , 2 -6 , and 2 -9 , respectively. The solution profile for different time levels is plotted in Fig. 1, which shows the same pattern reported in [2]. Due to unavailability of exact solution, we studied the accuracy via the stability of the method and noticed that as the value of β approaches zero, the solution diverges, and the scheme becomes unstable. From Figs. 1(a)-1(c) it is clear that the solution profile shows a proper behavior and the scheme is stable, that is, ρ(M) < 1, and as the spectral radius increases, the scheme goes to an unstable region, and the solution diverges, as shown in Fig. 1(d). We can also observe from the figure that for a fixed value of β, the graph tends to zero as time increases, and for various values of β, the curve propagates to the right and follows a sharp decay. Thus the present method shows the proper behavior of Eq. (45) for various values of β and t. The same behavior of this equation was illustrated by Mahanty and Sharma [1] and Mittal [2].
Example 4 Consider Eq. (1) with γ = 0 and α = m = 1 to obtain with initial and boundary conditions [1] We The same pattern of the problem has been reported by [8] and [1]. The accuracy of solution was discussed by means of the stability of the scheme. Example 5 In this case, we consider the coupled one-dimensional Burger Eq. (7) by taking μ = -1, ν = 2, and γ = 0, which leads to The exact solution is [47,48] where a o is an arbitrary constant, and A = 1 2 a 0 ( 4αβ-1 2α-1 ). The initial and boundary conditions are extracted from the exact solution. The numerical solution was computed for different values of α, β, and T = 0.5, 1 in the domain [-10, 10] with a 0 = 0.05 and dt = 0.001. For comparison, the error norms were computed and shown in Table 10. The stability of the scheme was discussed in terms of the spectral radius shown in the table. From the table it is clear that the proposed scheme is stable and produces better results even in large domains. The solution profile of Y and Z for T = 1 is plotted in Fig. 3. From the figure we can easily notice the betterment of the proposed technique.
with exact solution [22] Y (ξ , t) = e -t sin(ξ ), Z(ξ , t) = e -t sin(ξ ).  The initial and boundary conditions are extracted from the exact solution. The numerical solution was obtained for time levels T = 0.5, 1 in the domain [-π, π]. The obtained results were compared in the form of error norms with those of the differential quadrature method [22] and are shown in Table 11. The rate of convergence using CGL collocation points is shown in Table 12. From Table 11 we can observe that the present method is stable and produces a better solution than the available techniques. In Fig. 4 the numerical and exact solutions of Y and Z are plotted for T = 1. The figures reflect a good agreement of the obtained numerical result with exact solution.
Example 7 In this case, we take Eq. (7) with ν = 0, α = 0, β = 0, γ = 1, and μ = -1/Re, which leads to the following two-dimensional coupled Burger equation The exact solution is   computed for Re = 1, 10, and 100 and nodal points N = 100, 400, that is, N = (10 × 10) and (20 × 20), and compared with the error norms obtained by Arshad [4] using the meshfree technique presented in Table 13. From the table we notice that the present results are more accurate when Re = 1, whereas the accuracy decreases as Re increases with increasing spectral radius. The solution and error plots are shown in Fig. 5, which shows a kink-like behavior for Re = 50.
Example 8 Finally, we study two-dimensional coupled Burger equations (59)-(60) in the domain [0, 1] × [0, 1] with the following initial and boundary conditions taken from [4]: Y (ξ , η, 0) = sin(πξ ) sin(πη), Z(ξ , η, 0) = sin(πξ ) + sin(2πξ ) sin(πη) + sin(2πη) , The problem was solved using the proposed technique for nodal points N = 100, 400, that is, N = (10 × 10) and (20 × 20) at time t = 0.01, and Re = 1. Due to the nonavailability of the exact solution, the obtained results were compared at different collocation points with the numerical solution by the meshfree method [4] and finite element technique [49] shown in Table 14. From the table it is clear that the proposed method produces almost the same results as those of existing methods. The solution profile for different time levels t = 0, 0.01, 0.05 at fixed values of η are plotted in Fig. 6. From the figure we can observe that as the time increases, Z(ξ , η, t)) moves from the negative part to the positive one, and the graphs tend to zero. Similarly, a 3D plot of the solution is shown in Fig. 7.

Conclusion
In this paper, we studied a numerical method based on the Lucas polynomials and computed solutions of three different models, including the generalized Burger-Huxley equation, generalized Burger-Fisher equation, and one-and two-dimensional nonlinear coupled Burger equations. The dependent variable is approximated by the Lucas polynomials, whereas the Fibonacci polynomials are used for its derivatives. We discussed the stability of the proposed scheme in the form of spectral radius. For comparison of the proposed method, we computed the error norms in different domains and compared the results with