Numerical treatment of the generalized time-fractional Huxley-Burgers’ equation and its stability examination


 In this paper, we show how to approximate the solution to the generalized time-fractional Huxley-Burgers’ equation by a numerical method based on the cubic B-spline collocation method and the mean value theorem for integrals. We use the mean value theorem for integrals to replace the time-fractional derivative with a suitable approximation. The approximate solution is constructed by the cubic B-spline. The stability of the proposed method is discussed by applying the von Neumann technique. The proposed method is shown to be conditionally stable. Several numerical examples are introduced to show the efficiency and accuracy of the method.


Introduction
Fractional calculus is a generalization of classical calculus that is concerned with operations of differentiation and integration to fractional order. The concept of the fractional operator was first raised by Marquis de L'Hôpital and G. W. Leibniz in the year 1695. From now on there have been several fundamental works on the fractional derivative and fractional differential equations, and many books were performed in this field, in which we can refer to the books of Ross and Miller [1], Samko et al. [2], and Podlubny [3]. Several numerical methods have been developed to obtain an approximate solution of fractional differential equations. These methods include finite difference method [4,5], homotopy perturbation method [6][7][8], variational iteration method [9,10], Adomian decomposition method [11,12], and spectral methods [13][14][15][16]. The generalized Huxley-Burgers' equation is in the following form [17]: where β, ν, η, γ, and δ are parameters, δ 0 > , η 0 ≥ , and γ 0, 1 ( ) ∈ . Equation (1) is an extended form of the famous Huxley and Burgers' equations. If we take β 0 = and δ 1 = , equation (1) is reduced to the Huxley equation which describes nerve pulse propagation in nerve fibers and wall motion in liquid crystals [17]: If we take η 0 = , β 1 = , and δ 1 = , equation (1) is reduced to Burgers' equation which describes the far field of wave propagation in nonlinear dissipative systems [18]: It is known that the nonlinear diffusion equations (2) and (3) play important roles in physics. They are of special significance for studying nonlinear phenomena. If we take η 0 ≠ , β 0 ≠ , and δ 1 = , equation (1) becomes the following Huxley and Burgers' equation: Equation (1) shows a prototype model for describing the interaction among reaction mechanisms, convection effects, and diffusion transport. This equation was investigated by Satsuma [19] in 1987. The generalized time-fractional Huxley-Burgers' equation whose solutions will be approximated is of the form [20]: and is subject to the conditions: and We can obtain approximate solutions of the generalized time-fractional Huxley-Burgers' equation that have applications in various fields of science and engineering. The time-fractional Huxley-Burgers' equation illustrates intercommunications between reaction mechanisms and diffusion transports and an evolutionary model depicting nerve pulse propagation [21]. Many scholars have recently worked on the fractional Burger-Huxley equation due to the necessity of finding a solution to this equation. In [22], the space-timefractional Burgers-Huxley problem was solved using a Legendre spectral finite difference method. The author in [23] looked into solitary wave solutions to this problem. In [24], the problem is solved using the finite difference method. In [20], the residual power series method is used to obtain an analytical solution to this problem.
In this work, we have proposed a collocation method for solving the time-fractional Burgers-Huxley equation using the mean value theorem for integrals and cubic B-spline basis functions. In this method, we use cubic B-splines for spatial variables and their derivatives which produce a system of fractional ordinary differential equations. The fractional derivative is simplified by using the mean value theorem for integrals and the finite difference method. B-splines have two useful qualities. One distinguishing characteristic is that the continuity conditions are inherent. As a result, when compared to other piecewise polynomial interpolating functions, the B-spline is the smoothest. B-splines also have small local support, which means that each B-spline function is only non-zero over a few mesh subintervals, resulting in a tightly banded resultant matrix for the discretization equation. B-splines have special advantages due to their smoothness and capacity to handle local occurrences. In combination with the collocation method, this significantly simplifies the solution procedure of differential equations. There is a great reduction of the numerical effort because there is no need to calculate integrals (as in variational methods) to build the final set of algebraic equations that replaces the provided nonlinear differential equation. The current method does not require extra effort to deal with the nonlinear parts, unlike some previous techniques that used various transformations to reduce the problem to a simpler equation. As a result, the equations are readily and elegantly solved utilizing the current method. The Caputo fractional derivative of order α which is used in this paper is defined by, where Г( ) ⋅ is the Gamma function. This paper is organized as follows. In Section 2, a proposed method depends on the cubic B-spline collocation method and the mean value theorem is derived. In Section 3, we obtain the initial vector that is required to start the iterative process of the proposed method. In Section 4, the stability is theoretically discussed by using the von Neumann technique. In Section 5, we illustrate three numerical examples that are introduced to illustrate the efficiency and accuracy of the proposed method. Finally, the conclusion is given in Section 6.

Derivation of the method
To approximate u x t , ( ) by cubic B-spline collocation method over finite elements, let the region R a b , ] be discretized by a set of points R ij which are the vertices of a grid of points x t , be cubic B-splines with nodes at at these knots are given as follows: x h The collocation method for approximately solving equation (5) consists in seeking approximations ] which is spanned by the linearly independent set of cubic B-spline functions , , , having the form [25]: where w t i ( ) are the unknown time-dependent quantities. The quantities w t i ( ) are chosen such that U x t , ( ) satisfies equation (5) and its conditions (6) and (7) at the points x t , Then U x t , ( ) satisfies the following equations: Substitute equation (9) and the values in Table 1 into equation Using the Caputo fractional derivative definition, the time-fractional derivative can be discretized as By using a piecewise technique, equation (14) becomes , the mean value theorem for integrals [26] will be applied as follows: This enables us to write the Caputo fractional derivative in the following form: Numerical treatment of the generalized time-fractional Huxley-Burgers' equation  439 with a local truncation error Substituting equation (16) in equation (13), to obtain the following system: From equation (18) we obtain the following: For simplicity, equation (19) can be rewritten as a recurrence relationship: ) ( ) + × + matrix system, when we add the boundary conditions. From equations (20) and (21), we obtain the following system of linear equation: ,a n d .

Initial vector
The initial condition helps us to obtain the initial vector w w w w w , , , , , Using the values in Table 1 into equation (23) gives The last system (24) consists of N 1 ( ) + equations and N 3 ( ) + unknowns. To solve the system of linear equation (24) we need two additional equations. These equations are obtained as follows: By using Table 1, in equation (25) we obtain Combining equations (24) and (26) we obtain the following system of linear equation: where which enables us to obtain the initial vector w w w w w , , , , , .

Stability analysis
The von Neumann technique is suitable for linear problems with constant coefficients. To study the stability of our numerical scheme by using the von Neumann technique, we must linearize the nonlinear terms where q 1 2 = − , h is the element size, φ is the mode number, and ξ j is the amplification factor at time level j. Inserting the latter expression for w i j into the system (  where Ω , and ϕ φh = . For j 1 = in equation (31), we obtain thus we obtain the following: where Ω Ε Η Ζ If we put j 2 = in equation (31), we obtain , and τ α 2,1 are positive quantities. Using inequality (34) implies For j 3 = in equation (31), we obtain , and τ α 3,1 are positive quantities, using inequalities (34) and (35), we obtain Thus, by the same method, we can prove that ξ ξ w g , j 1 ≥ , then the proposed numerical scheme is conditionally stable such that h is chosen to be small enough.

Application and discussion
In this section, we employ the suggested method to solve some examples, and we will show that the method produces a good approximation. Our proposed scheme's accuracy is measured by computing the l 2 error norm and maximum absolute error for several choices for α. Error norms are defined as The computations associated with the experiments were performed in the Mathematica software package on a PC, CPU 2.8 GHz.
Where the source term on the right-hand side is given by The exact solution to this problem is given by The numerical results are presented in Tables 2 and 3, which show a comparison between the approximate and exact values at different values of γ and α with h 0.01 = and k 0.01 = . Tables 4 and 5 show a comparison between the maximum absolute error and l 2 error norm at different values of γ and different time levels. Table 6 presents the maximum absolute errors at different values of N . From our numerical results, one can conclude that the numerical solutions are in a very close agreement with the exact solution. Figure 1a demonstrates that the numerical and exact results are indiscriminately same for t 0.25 = and are shown in Figure 2.
and the initial conditions Where the source term on the right-hand side is given by 2 sin 2π Γ 3 4π sin 2π 1 sin 2π sin 2π sin 2π .
The exact solution to this problem is given by The numerical results are presented in Tables 7 and 8 The exact solution to this problem is given by The numerical results are presented in Table 12, which show a comparison between the approximate and exact values with α 0.75 = , h 0.01 = , and k 0.01 = . Table 13 shows a comparison between the maximum absolute error and l 2 error norm at different time levels. Table 14 presents the maximum absolute errors at different values of N . Table 15 illustrates the comparative performance of our method and that of [29]. Also, the results conclude that the errors in the proposed method are less than or similar to those of [29]. From our numerical results, one can conclude that the numerical solutions are in very good agreement with the exact solution. Figure 5a shows that the exact and numerical results are similar for t 0.25 = and  are shown in Figure 6.

Conclusion
The proposed numerical method is presented for solving the generalized time-fractional Huxley-Burgers' equation. We discuss the stability analysis of the proposed scheme by applying the von Neumann technique. The proposed scheme is shown to be conditionally stable. The obtained results using this method are more acceptable than others. Our results illustrate the proposed method's efficiency and accuracy for solving the nonlinear generalized time-fractional Huxley-Burgers' equation.