Numerical study of 1D and 2D advection-diffusion-reaction equations using Lucas and Fibonacci polynomials

In this work, a numerical scheme based on combined Lucas and Fibonacci polynomials is proposed for one- and two-dimensional nonlinear advection–diffusion–reaction equations. Initially, the given partial differential equation (PDE) reduces to discrete form using finite difference method and θ-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\theta -$$\end{document} weighted scheme. Thereafter, the unknown functions have been approximated by Lucas polynomial while their derivatives by Fibonacci polynomials. With the help of these approximations, the nonlinear PDE transforms into a system of algebraic equations which can be solved easily. Convergence of the method has been investigated theoretically as well as numerically. Performance of the proposed method has been verified with the help of some test problems. Efficiency of the technique is examined in terms of root mean square (RMS), L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_2$$\end{document} and L∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_\infty $$\end{document} error norms. The obtained results are then compared with those available in the literature.

where g is source term depends on space and time, , ∂ are spatial domain and boundary of the domain, respectively, α and γ are positive constants, ∇ = ∂ ξ + ∂ η and = ∂ ξξ + ∂ ηη represents Laplacian operator. When α = 0 and β = 1 Eq. (1) becomes heat equation which is used to model thermal conductivity over a solid body. This equation is also used to model other real-world phenomenons like Black-Scholes model; it is applied for price estimation, and in chemical process it is designed in the form of diffusion equation [18]. In the form of advection-diffusion equation it describes heat transfer in a solid body and dissipation of salt in ground water. In image processing it handles image at different scales [28]. The solution of this equation plays a significant role to examine behavior of such physical processes. Different numerical methods have been used for solution of heat equation by many researchers. A comparative study between the classical finite difference and finite element method was investigated by Benyam Mebrate [25]. B-spline and finite element method have been used by Darbel et al. [12]. Hooshmandasl et al. [18] used wavelet methods to solve one dimensional heat equations, while Dehghan [13] used second order finite difference scheme for two-dimensional heat equations. In [34], the authors applied haar wavelet method for two-dimensional fractional reaction-diffusion equations. In [2], the authors applied fifth-kind Chebyshev polynomials for the spectral solution of convection-diffusion equation.
When α = 1 and β = 1/Re ('Re' is Reynolds number), then Eq. (1) becomes well-known Burger equation which was experienced in the field of turbulence, shock wave theory, viscous fluid flow, gas dynamics, cosmology, traffic flow, quantum field, and heat conduction [27]. The low-kinematic viscosity shocks and the relation between cellular and large-scale structure of the universe have been described by one-and threedimensional Burger equations. When traffic is treated as one-dimensional incompressible fluid, then the density wave in traffic flow which changes from non-uniform to uniform distribution is described by Burger equation [10]. The Burger equation was first introduced by Bateman in viscous fluid flow, which was then extended by Burger in (1948) to examine turbulence phenomena and that is why it is known as Burger equation [27]. Due to wide applications of Burger equation many numerical methods have been implemented to study behavior of the model. One-dimensional Burger equation has been solved using various techniques in [23,35]. Mittal and Jiwari [27] implemented differential quadrature method for solution of Burger-type equation. Similarly El-sayed and Kaya [20] solved two-dimensional Burger equation using decomposition method. Liao [24] used fourth-order finite difference technique for the study of two-dimensional Burger equation. Oruc [33] applied meshless pseudo-spectral method to modified Burgers equations. The same author in [30] studied three-dimensional fractional convection-diffusion equation using meshless method based on multiple-scale polynomials. In this work, we study afore mentioned models by using combination of Lucas and Fibonacci polynomials. These polynomials can be directly obtained as a special case from the work done in [4,6]. These polynomials are nonorthogonal and do not require domain and problem transformation which is important point of the proposed scheme. Also the higher order derivative of unknown functions can be easily approximated via Lucas and Fibonacci Polynomials. Second, it is straightforward and produces better accuracy for less number of nodal points. Many researchers applied these polynomials for the solution of fractional differential equations (FDEs) such as Elhameed and Youssri [3], who applied Lucas polynomial in a Caputo sense to FDEs. Moreover they computed the solution of fractional pantograph differential equations (FPDEs) using generalized Lucas polynomials [39]. Other polynomials applied for approximation of FDEs are studied in [5,7]. Cetin [11] used Lucas polynomial approach to study a system of higher order differential equation, whereas Bayku [9] applied hybrid Taylor-Lucas collocation technique for delay differential equations. Mostefa [29] obtained solution of integro differential equation using Lucas series. Farshid et al. [26] applied Fibonacci polynomials for solution of Voltera-Fredholm integral equations. Omer Oruc [31,32] applied a combined Lucas and Fibonacci polynomials approach for numerical solution of evolutionary equation for the first time. Recently, Ali et al. [8,15,17] applied Lucas polynomials coupled with finite differences and obtained accurate solution of various classes of one-and two-dimensional PDEs. In this paper, we implement the proposed method to one-and two-dimensional Burger and heat equations. The simulation is carried out with the help of MATLAB 2013 using Intel core-i7 machine with 4GB RAM. The error bound of the scheme is also investigated in this work. The paper is organized as follows: In Sect. 2 we define basic definitions and important concept which will be used in this work. In Sect. 3 methodology of the proposed scheme is formulated. Error analysis of the method is described in Sect. 4. Numerical experiments are presented in Sect. 5 followed by conclusion of the paper.

Basic concepts and definition
In this section, we describe some necessary definitions and concepts required for our subsequent development. Definition 2.1 [22] Lucas polynomials are the generalization of well-known Lucas numbers which are generated by linear recursive relation as follows: with L 0 (ξ ) = 2 and L 1 (ξ ) = ξ , by putting ξ = 1, Eq (2) gives Lucas numbers.
Function Approximation. Let Y (ξ ) be square integrable on (0, 1) and suppose that it can be expressed in terms of the Lucas polynomials given as follows: Similarly the mth order derivative of the function Y (ξ ) can be approximated in terms of finite Lucas series is given as follows: where D is (N + 1) × (N + 1) differentiation matrix such that: in which d is square matrix of order N × N defined as: For example, if we choose N = 3, then we have and for m = 2 the second-order derivative of L(ξ ) can be obtained using Eq. (6) where k = [0, 1, 2, 3] after element wise multiplication we get Now assume a two-dimensional continuous function Y (ξ, η) can be written in terms of Lucas polynomials as follows:

Solution methodology
Consider the following evolutionary equation: where £ is differential operator and g(x, t) is a given smooth function. The initial and boundary conditions are given as follows: where B is boundary operator. To approximate Eq. (10), first we define time discretization given by the following: where δt = T /M is time step size for variable t and T is final time. Now applying finite difference scheme to temporal part and θ −weighted scheme to spatial part of Eq. (10), one can write where Y n+1 (x) = Y (x, t n+1 ) and so forth. By putting θ = 0.5, the scheme in Eq. (12) represents Crank-Nicolson which is O(δt 2 ) accurate in time.
For discretization of spatial domain = [a, b] we use regular grid point defined as follows: For (n + 1) time level and at the collocation point x = (ξ i , η j ) Eq. (9) can be written as follows: Using Eq. (13) in Eq. (12), we have The above equation can also be written as The boundary conditions (11) Matrix form of Eqs. (14) and (15) is given by For k, m = 0, . . . , N elements for the above matrices can be obtained in the following: where H, G and B are (N + 1) 2 order square matrices. The unknown coefficient vector C can be obtained by solving Eq. (16). Once the values of unknown coefficient are computed the solution of the problem under consideration can be obtained from Eq. (9).

Truncation error estimate
To study the error estimate of the proposed scheme, we follow the approach of Abd-Elhameed and Youssri [3].
Proof Consider the absolute error between exact and approximate solution Then the truncated term is given as It is shown in [3] where ϑ is well-known golden ratio. Therefore, Eq. (22) implies that where κ = Pϑ. The above inequality can also be written as Here, (N + 1, κ) is the incomplete gamma function and (N + 1) is complete gamma function [36]. In integral form Eq. (23) can be written as As, exp(−t) < 1, therefore, for all t > 0, we have This completes the proof.

Numerical examples
In this section, one-and two-dimensional Burgers' and diffusion equations have been solved using the proposed method. Performance of the technique is examined by computing L 2 , L ∞ , and root mean square (RMS) error norms for different collocation points N and time T . The obtained results are then compared with available results in the literature.
with initial and boundary conditions and exact solution [18].
Comparison of the above equations with Eqs. (10) and (11) gives Applying the technique discussed in Sect. 2, Eq. (14) takes the following form: where L k (ξ ) represents second-order derivative of Lucas polynomials which can be obtained using Eq. (6). The matrices in Eq. (17)- (19) takes the following form: The problem has been solved for different values of nodal points N, and time T. Computed solutions are compared with the results provided by Chebyshev wavelet method in the form of various error norms which are shown in Table 1. It is clear from the table that proposed scheme gives excellent results than cited work. Numerical convergence and Cpu time have been reported in Table 2 Table 3 and better results noticed than that for small time level. Solution profile and absolute errors are presented in Fig. 1 which show that exact and approximate solution are well matched with each other showing efficiency of the proposed scheme.
Example 5.2 Consider Eq. (1), with α = 0, β = 1, γ = −2 and g(ξ, t) = 0, defined as homogeneous heat equation given as follows Initial and boundary conditions are extracted from exact solution. The problem has been solved by adopting the procedure discussed in Sect. 2. The results are computed for different values of N and T . RMS, L 2 , and L ∞ error norms have been calculated and comparison is made with available results in literature [18] for different values of T, N = 16, δt = 0.001 that are shown in Table 4. From the table, it is straightforward that the results achieved using the proposed method are better than those available in literature which show proficiency of the method used. For convergence in space the results obtained are shown in Table 5 for different values of dξ showing that the solution converges as the value dξ decreases. Graphs of solutions along with absolute errors for T = 1 and N=15 are plotted in Fig. 2. The figure shows that exact and approximate solutions are closely matched.
There are two cases. Case 1: In this case the exact solution is given as [16] Y (ξ, η, t) = 1 Initial and boundary conditions are extracted from the exact solution. The nonlinear part in Eq. (29) can be linearized by the following formula [8]: Applying the technique discussed in Sect. 2 Eq. (14) takes the following form: The matrices H, G and B in Eq. (16) for k, m = 0, 1, . . . , N , are given as follows: The approximate solution of the problem has been computed for different values of viscosity β, time T , and nodal points N . Error norms L 2 , L ∞ and RMS have been calculated and the obtained results compared with the results of Haar wavelet [16] and differential quadrature [19]. The results are presented in Tables 6 and 7. From these tables it is obvious that proposed method works pretty well than those available in the literature. The solution profile and error plots for T = 2 , δt = 0.001, β = 1, and nodal points [20 × 20] are shown in Fig. 3 showing efficiency of the proposed technique.

Case 2:
In this case exact solution of Eq. (29) is given as [19] Y Here also initial and boundary conditions are extracted from the exact solutions. The method discussed in previous example is applied to solve this example for different value of viscosity β. Here also the error norms L 2 , L ∞ and RMS have been computed and are compared with the results of meshless collocation method [19] available in literature for different values of N . The obtained results presented in Table 8 which shows that the present method gives better results than those available in literature. The graph of numerical and exact solution are shown in Fig. 4 which shows efficiency of the current technique.
Example 5.4 Finally, consider the case when α = 0, β = 1, γ = 0 and g(ξ, η, t) = 0; Then Eq. (1) takes the following form: The initial and boundary conditions are extracted from exact solution [37] Y (ξ, η, t) = sin(πξ ) sin(πη)   Tables  9 and 10. From these tables it is clear that the results got using the proposed technique are comparable with those of multiquadric (MQ) RBFs and better than those of wendland (WL) RBFs [37]. The solution profile is plotted in Fig. 5 when T = 0.2 showing efficiency of the method suggested.

Conclusion
In this paper, we studied a numerical technique based on Lucas and Fibonacci polynomials. First, we discretized temporal part of PDEs by finite difference and spatial part by θ − weighted scheme with θ = 1/2 (Crank Nicolson method). Thenceforth, the unknown functions are expanded by Lucas series while their derivatives are replaced by Fibonacci polynomials. Performance and convergence of the method are investigated by several   test problems including one-and two-dimensional linear and nonlinear equations. The results are compared with exact as well numerical results available in the literature. The comparison of the results justify demonstrates efficiency and applicability of the proposed methodology.