Abstract

A new application of the hybrid generalized differential transform and finite difference method is proposed by solving time fractional nonlinear reaction-diffusion equations. This method is a combination of the multi-time-stepping temporal generalized differential transform and the spatial finite difference methods. The procedure first converts the time-evolutionary equations into Poisson equations which are then solved using the central difference method. The temporal differential transform method as used in the paper takes care of stability and the finite difference method on the resulting equation results in a system of diagonally dominant linear algebraic equations. The Gauss-Seidel iterative procedure then used to solve the linear system thus has assured convergence. To have optimized convergence rate, numerical experiments were done by using a combination of factors involving multi-time-stepping, spatial step size, and degree of the polynomial fit in time. It is shown that the hybrid technique is reliable, accurate, and easy to apply.

1. Introduction

The nonlinear reaction-diffusion equations have found numerous applications in pattern formation, in many branches of biology, chemistry, and physics [14]. Reaction-diffusion (RD) equations have also been applied to other areas of science and can be successfully modelled by the use of fractional order derivatives. [518]; for example, the RD equations are employed to describe the CO oxidation on Pt (110) [5], the study of waves on Xenopus oocytes [11], and the study of reentry in heart tissue [7, 13]. A great deal of effort has been expended over the last 10 years in attempting to find robust and stable numerical and analytical methods for solving fractional partial differential equations of physical interest. There has also been a wide variety of numerical methods, for example, finite difference techniques, finite element methods, spectral techniques, adaptive and nonadaptive algorithms, and so forth, which have been developed for RD’s numerical solution [19, 20].

The differential transform method was used first by Zhou [21] who solved linear and nonlinear initial value problems in electric circuit analysis. This method constructs an analytical solution in the form of a polynomial. It is different from the traditional higher order Taylor series method, which requires symbolic computation of the necessary derivatives of the data functions. The Taylor series method computationally takes long time for large orders. The differential transform is an iterative procedure for obtaining analytic Taylor series solution of ordinary or partial differential equations. The method is well addressed in [2226]. Recently, the application of differential transform method is successfully extended to obtain analytical approximate solutions to ordinary differential equations of fractional order. Fractional differential transform method (FDTM) is a method that Arikoglu and Ozkol [23] developed for solving linear and nonlinear integrodifferential equations of fractional order. This method solves problems with high accuracy while constructing semianalytic solutions in the polynomial forms. FDTM is based on classical differential transform method, fractional power series, and Caputo fractional derivative [27]. Arikoglu and Ozkol [23] tested their approach on several examples and the results obtained are in good agreement with the existing ones in the open literature. Momani et al. [28] presented a new generalization of the differential transform method that extended the application of the method to differential equations of fractional order. The new technique is named generalized differential transform method (GDTM) and is based on one-dimensional differential transform, generalized Taylor’s formula [29], and Caputo fractional derivative [27].

Fractional partial differential equations (FPDEs) are also an interesting and an important topic. The fractional derivatives and integrals have been occurring in many physical and engineering problems with noninteger orders. Fractional calculus is based on the definition of the fractional derivatives and integrals. They play a major role in engineering, physics, and applied mathematics. FPDEs are used to model complex phenomena since the fractional order differential equations are naturally related to the systems with memory and nonlocal relations in space and time which exist in most physical phenomena. Fractional order differential equations are as stable as their integer order counterpart. One of the fundamental equations of physics is the Schrödinger equation which describes how the quantum state of physical system changes with time. The fractional Schrödinger equation provides us with a general point of view on the relationship between statistical properties of quantum mechanical path and structure of fundamental equations of quantum mechanics [30]. In other words, the fractional quantum mechanics includes the standard quantum mechanics as a particular case. So FPDEs are obtained by generalizing differential equations to an arbitrary order. There are three popular methods for seeking approximate solutions for FPDEs which are the finite difference, finite element and spectral methods. In the literature there are many papers on these three methods. In these papers, the authors proposed the use of least squares finite element solution and fully discrete Galerkin method to solve nonlinear space fractional partial differential equations [20, 31].

The differential transform is well suited to combine with other numerical techniques, as shown by Yu and Chen [32] who applied the hybrid method to solve the transient thermal stress distribution in a perfectly elastic isotropic annular fin. Kuo and Chen [33] employed the hybrid method to solve Burger’s equation for flow systems with high Reynolds numbers. This method was also employed to analyze the dynamic response of an electrostatically actuated micro fixed-fixed beam [34].

In the current study, the hybrid generalized differential transform/finite difference method is used for solving time fractional nonlinear RD equations. The validity of the proposed approach has been confirmed by comparing the results derived in the literature using the GDTM method [19], homotopy perturbation method (HPM) [35], and fractional variational method (FVIM) [36].

There are several approaches of definitions for the fractional derivative. Among them, one is called Riemann-Liouville fractional derivatives and defined byHere is the -order Riemann-Liouville integral operator which is expressed as follows:If we use this definition, we must know the initial value of some fractional order derivative of the unknown function or we must have homogenous initial conditions. Unlike the Riemann-Liouville approach in the case of the Caputo derivative there are no restrictions on the initial conditions. Thus the following definition is used in this study:This operator is generally called “-order Caputo differential operator” [37, 38].

For Caputo derivative we have

2. Generalized Differential Transform/Finite Difference Method

We firstly introduce the main features of GDTM [19, 23, 28, 3941] according to the generalized differential transform of the th derivative of a function of one variable defined as follows:where ,    times. In (5), is the original function, is the transformed function, and the differential inverse transform of is defined as follows:In case of , GDTM reduces to the classical DTM. From definitions (2) and (3) all fundamental properties of GDTM can be obtained easily [19, 39]. Since has been proved from the definitions of fractional calculus, the fractional solutions reduce to the standard solution .

In this study, we use a hybrid method that is a combination of generalized, temporal differential transform and spatial finite difference methods to solve nonlinear fractional reaction diffusion equations.

We present a solution of a more general model of RD equationwhere is the diffusion coefficient and is a nonlinear function. We consider two different forms of which are called time fractional Fisher equation and time fractional FitzHugh-Nagumo equation.

We apply GDTM to discretize fractional order time derivative and central difference method to discretize derivatives in direction, respectively. After transforming (7) using the GDTM, we get the following equation:where is the transformed function of and is the transformed function of . The region is divided into several equal intervals and each interval has a width . The time interval of interest is discretized using a time step . After discretization of the equation we get solution at time and these results are adopted as the initial values for the next time interval. This time-stepping procedure assists in obtaining a converged solution to a desired accuracy [39]. An appendix has been added to the paper to show that the error, if any, in the method is bounded. This implies stability of the scheme and by Lax equivalence theorem it thereby implies convergence.

The new algorithm has been developed to solve the nonlinear reaction diffusion equation and our aim of this approach is to combine the flexibility of differential transform and the efficiency of finite differences. This algorithm also provides an iterative procedure to calculate the numerical solutions; therefore, it is not necessary to carry out complicated symbolic computation. On applying the differential transform method with respect to time on the equation we are basically transforming the time-evolutionary equation to an elliptic type. In essence this means that the central finite difference approximation that is subsequently used on the transformed equation is a Poisson solver. The resulting system of linear algebraic equations is then diagonally dominant and hence the Gauss-Seidel iterative method used for solving the same has assured convergence as the coefficient matrix remains nonsingular throughout the computation. The algorithm used thus succeeds in segregating the time discretization from explicitly influencing the computation in the spatial domain and this presents a situation wherein the two can be handled independent of each other in the course of computation without having to bother about the stability of the solution if the differential transform part is properly handled. The latter is achieved deftly in the differential transform part of the algorithm by using the multistepping procedure as first enunciated by Yu and Chen [32] in their phenomenal work and used subsequently by Odibat et al. [29, 39]. In summary, this means that convergence is never in doubt in the algorithm but slow convergence can be if the time and spatial discretizations are badly handled. The convergence is optimized in the paper computationally by proper selection of time step in the differential transform part of the algorithm and then the spatial step size in the finite difference part of the algorithm. Next important step in the algorithm is the decision on the number of terms to be adopted in the inverse differential transform that gives us the solution of the problem as a power series in time. The aforementioned three vital components of the algorithm have been meticulously handled and a brief summary of the numerical experiment undertaken concerning the same is presented in a table. The numerical study recommends that the combination of 5, 10, and 50 spatial step size with a time step of 0.0005 or 0.001 assures the best rate of convergence if we take minimum ten terms in the time series.

3. Illustration of Generalized Differential Transform and Finite Difference Method

To show effectiveness of the proposed numerical solution using the temporal generalized differential transform and the spatial finite difference method and to give an understandable overview of the methodology, two examples of the reaction diffusion equations will be discussed in the following section. Then our results will be compared with published work of Rida et al. [19] in which GDTM was used to solve the same equations.

Example 1. The time fractional Fisher equation is

In this example, we have the nonlinear function . The initial condition used is

Operating the generalized differential transform on (9) gives us the following equation:where is the generalized differential transform of .

Now we apply the central finite difference method to the derivatives with respect to and this gives us The initial condition on discretization yields

Equation (12) is a recurrence relation. The time series solution of the given equation is then obtained by using (12) and (13) with to obtain . Some of are recorded in Table 1.

The time series solutions of (12) with the initial condition (13) are obtained as follows:

The numerical calculation results are shown in Figures 1 and 2, respectively. Our results are in agreement with the published work of Rida et al. [19], who considered the same equation. An exact solution of the standard form of Fisher equation for isThe comparison of our results with the exact solution is shown in Figure 1, for , and quite clearly good agreement is found.

Approximate solutions are shown in Figure 2 for and .

The influence of on the function is shown in Figure 3. This figure indicates a decrease in the fractional order by choosing the fixed that corresponds to an increase in the function and also indicates a slow diffusion for the values of and and a fast diffusion for the values of , respectively. It is clearly seen that increase for with the increases in .

Numerical comparison between GDTM [19], HPM [35], FVIM [36], and hybrid method are found in Table 2 which shows hybrid method is more promising.

It is also found that the result is in complete agreement with the result of HPM [42, 43] and ADM [44] for .

We investigate convergence criteria of our solutions for different values of and . To illustrate this, we compared our results with the analytical solution in case of . Here is order of differential transformation method and denotes the number of terms to be calculated.

In Figures 4, 5, and 6 the difference between the results obtained in this study and the results of the analytical solution is of the order of . This is a pointer to the fact that there is convergence and is a restatement in numerical terms of what was shown in the Appendix.

One important observation made from the computation is that when the number of mesh points was increased, less number of terms was required in the time series solution to have convergence for a predetermined accuracy. The hybrid method of the present study gives faster convergence than other traditional methods; for example, if we take (mesh point is 50) then the solution converges for . We now consider another example.

Example 2. The time fractional FitzHugh-Nagumo equation isIn this type of equation, the nonlinear function depends on and it is . The initial condition isUsing the hybrid method on the above initial boundary value problem (IBVP), as done in the previous example, we getUsing second order finite difference method, the boundary values were obtained as follows:Table 3 presents some of the ’s.

The time series solution for the above IBVP at different times isNumerical solutions for the time fractional FitzHugh-Nagumo equation with various values are shown in Figures 7 and 8. A comparison of the results in a limiting case wherein an analytical solution exists is shown in Figure 7. The results are in close agreement with those of Rida et al. [19] for the same equation.

For , it can easily be seen that the exact solution of FitzHugh-Nagumo equation isFigure 9 is prepared to show the influence of on the function . It is clearly seen that decrease for with the decreases in .

As shown in the Table 4, our results show close agreement with the exact solution and agree with those of Rida et al. [19]. From this table it is clear that the present work gives better approximation than GDTM as we increase .

Numerical comparison between GDTM, FVIM, and hybrid method is shown in Table 5, which indicates hybrid method is more promising.

4. Conclusion

Many real physical problems can be best modelled with fractional differential equations but the fact is when the equation is nonlinear there are very few reliable methods. The numerical methods that can be used to solve fractional differential equations are known to have problems of convergence and stability. These aspects are well addressed in the paper by suggesting a new procedure that uses a combination of the generalized differential transform and central difference methods. The Appendix clearly spells out the fact that the error as a result of discretization and computation is bounded and hence implies stability of the method. Lax equivalence theorem further implies convergence of the scheme. Two time fractional nonlinear reaction-diffusion equations considered for illustration of the hybrid method highlight the usefulness of the method in obtaining the solution of IBVPs involving time fractional derivatives. The control of convergence through a judicious choice of time and spatial step sizes and also the number of terms in the time series solution spells assured convergence. The segregation of the time domain from the spatial domain in the solution method ensures the fact that problem of stability does not arise. Diagonal dominance of the coefficient matrix in the system of linear algebraic equations resulting from the use of the central difference approximation in the Poisson equation ensures the fact that the matrix remains nonsingular during iterations and hence has assured convergence. An appropriate computational decision on the number of terms to be taken in the time series solution results in a convergent solution with fast convergence. Excellent comparison of the present results with the previous works on generalized differential transform method [19] and homotopy perturbation method [35] and fractional variational iteration method [36] provides confidence in the methodology adopted for the solution of time fractional differential equations.

Appendix

Estimation of Bounds on Truncation Error

Consider the fractional differential equations (9) and (16) in a general form asThe differential transform of (A.1) at the spatially discretized points gives uswhere is the differential transformed function of and is that of . Let us further denote as and its transform by . In this notation, (A.2) reads asWe now follow Jang et al. [45] and move on to arrive at an estimate on the bounds for the truncation error in a general way by considering the Taylor series expansion of :where , is the remainder. Let denote an approximate solution to that satisfiesClearly the difference between and is of . Thus, the local error isLet us now suppose thatUsing (A.7) in (A.6), we getWe now consider more numbers of terms in Taylor expansion (A.4); that is,As done earlier, let us denote by the following expression:Again, as earlier, let us suppose thatSimilar to (A.8), we now getWe so far addressed the local error due to two different truncations in the time series. In what follows we estimate the bounds on the cumulative error that includes the error discussed above.

Let denote the solution of (A.2). The local error in relative to isSince is a better approximation than , we may assume thatIn view of (A.14), we now have Using (A.8) and noting that is quite small in (A.15), we may take to be Thus is the bound on the tolerance in the to-be-obtained solution. When using different number of terms in the Taylor series expansion, earlier we denoted the solutions using a time step by and , respectively.

The paper uses an adaptive step size in computing the results. This is because such a procedure succeeds in keeping the error bounded and ensures convergence as a consequence of Lax equivalence theorem. To see what the adaptive step size produces and to show how such a procedure keeps the error bounded we start with the premise that is the most appropriate step size for the problem. This step size is determined using the definition of inverse differential transform:In our actual calculation we will not be able to consider infinite number of terms. We consider “” terms in respect of and “” terms in respect of .

Thus To write down a simpler expression we change the summation index from to . So we have from (A.18) the following:Thus, (A.16) on using (A.19) may be written as Using yet another step size , also relation (A.20) is satisfied. Let be the solution using . So from (A.8) we now haveUsing (A.20) and (A.21), we may writefor ,Thus, the above proceedings tell us that if criterion (A.23) is satisfied, then the error is bounded. In effect, this means that the scheme is convergent in lieu of Lax equivalence theorem. In our computations has been always chosen to satisfy inequality (A.23).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors are thankful to Ondokuz Mayıs University, Samsun, Turkey, for providing financial support to carry out this work under a major research project (Grant no. pyo.fen.1901.13.003).