Applying Discrete Homotopy Analysis Method for Solving Fractional Partial Differential Equations

In this paper we developed a space discrete version of the homotopy analysis method (DHAM) to find the solutions of linear and nonlinear fractional partial differential equations with time derivative α (0<α≤1). The DHAM contains the auxiliary parameter ℏ, which provides a simple way to guarantee the convergence region of solution series. The efficiency and accuracy of the proposed method is demonstrated by test problems with initial conditions. The results obtained are compared with the exact solutions when α=1. It is shown they are in good agreement with each other.

Various techniques have been investigated to solve partial differential equations of fractional order, such as the homotopy analysis method (HAM) [17][18][19][20][21][22], homotopy perturbation method (HPM) [23][24][25][26], Adomian decomposition method (ADM) [27][28][29], meshless method [30][31][32][33], operational matrix [34,35] and so on. In 1992, Liao introduced the homotopy analysis method, a semi-analytical method, for solving strongly nonlinear differential equations [20]. The main advantage of HAM is that it provides great freedom to choose equation type and solution expression of related linear high-order approximation equations. HAM gives rapidly convergent successive approximations of the exact solutions, if such a solution exists, otherwise approximations can be used for numerical purposes. It is an analytical approach to get the series solution of linear and nonlinear partial differential equations. Unlike the other analytical techniques, HAM is independent of small/large physical parameters. Since HAM has many advantages in comparison to other analytical methods, it is employed to solve continuous problems. Hence after the discrete ADM method [36], the discrete homotopy analysis method (DHAM) was introduced in 2010 by Zhu et al. [37]. This method can be applied to complex problems containing discontinuity in fluid characteristics and geometry of the problem. In addition, it needs little computational cost as numerical method in comparison to HAM; as an analytical approach DHAM has similar advantages to continuous HAM. By means of introducing an auxiliary parameter one can adjust and control the convergence region of the solution series. This method should be employed for solving various differential equations to highlight its high capabilities in comparison with other numerical methods.
The well-known property of Riemann-Liouville operator J α is

Definition 4 ([39]
). Form to be the smallest integer that exceeds α > 0, the Caputo fractional derivative of u(x, t) with respect to t of order α > 0 is defined as Note that the Caputo fractional derivative is considered due to its suitable for initial conditions of the differential equations.
The relations between Riemann-Liouville operator and Caputo fractional differential operator are given as follows

Discrete Homotopy Analysis Method
Consider the following general difference equation respect to j where N is a linear or nonlinear operator, j and t denote the independent variables. Suppose that ∆x = h and the function u(x, t) = u(j∆x, t) is the discrete function and denoted by u j (t). For simplicity, we ignore all boundary or initial conditions, which can be treated in the similar way. Similarly to continuous HAM, we first construct the so-called zeroth-order deformation equation by means of the discrete HAM (DHAM) where p ∈ [0, 1] is an embedding parameter, = 0 is an auxiliary parameter, L is an auxiliary linear operator, u j,0 (t) is an initial guess of u j (t), H j (t) denotes a nonzero auxiliary function, ϕ j (t; p) is an unknown function about j, t, p. It is important that one has great freedom to choose auxiliary things in (2). Obviously, when p = 0 and p = 1, it holds ϕ j (t; 0) = u j,0 (t) = u j (0) and ϕ j (t; 1) = u j (t), respectively. Thus, as p increases from 0 to 1, the solution ϕ j (t; p) varies from initial guess u j,0 (t) to the solution u j (t). Expanding ϕ j (t; p) in Taylor series with respect to p, we have Similarly continuous HAM by Liao [21], if the auxiliary linear operator, the initial guess, the auxiliary parameter , and the auxiliary function are so properly chosen, the series (3) converges at p = 1, then we have u j (t) = u j,0 (t) + +∞ ∑ m=1 u j,m (t).
which is used in the discrete homotopy perturbation method (DHPM) [16], where as the solution obtained directly, without using Taylor series. According to Equation (4), the governing equation can be deduced from the zeroth-order deformation Equation (DHAM). Define the vector → u n = u j,0 (t), u j,1 (t), · · · u j,n (t) .
Differentiating the zeroth-order deformation Equation (2) m times with respect to the embedding parameter p and then setting p = 0 and finally dividing them by m!, we obtain the following mth-order deformation equation: where and X m = 0, m ≤ 1, 1, m > 1. It should be emphasized that it is very important to ensure that Equation (3) converges at p = 1 otherwise, the Equation (5) has no meaning.
Theorem 1 (Convergence Theorem). As long as the series (5) is convergent, where u j,m (t) is governed by the high deformation Equation (6). It must be the solution of the original Equation (1).
and it holds lim n→∞ u j,n (t) = 0.
From the mth-order deformation Equation (6) and by using the definition of X m, it yields On the other side, according to the definition (7), we have In general, ϕ j (t; p) doesn't satisfy the original Equation (1). Let denote the residual error of Equation (1). Obviously, corresponds to the exact solution of the Equation (1). According to the above definition, the Maclaurin series of the residual error ε j (t; p) with respect to the embedding parameter p is Entropy 2018, 20, 332 5 of 15 when p = 1, the above expression gives That is, according to the definition of ε j (t; p) we have the exact solution of the original Equation (1) when p = 1. Thus as long as the series is convergent, it must be the solution of the original Equation (1).

Examples
Example 1. Consider the time fractional discrete diffusion equation Discrete diffusion equation is widely used in applied sciences. For example, population growth modeled by geographical spread [40], to model ionic diffusion on a lattice [41], and so on. Moreover, the entropy production was calculated for fractional diffusion Equation [42,43].
The standard central differences D h u j (t) and D 2 h u j (t) are defined by Initial value problem (13) and (14) is the discrete form of initial value problem for diffusion equation where D α t u(x, t) is Caputo fractional derivative of order α. To solve Equation (13) by DHAM let us consider the following linear operator: with the property that where c is constant coefficients. We define the nonlinear operator as Using the above definition, we construct the zeroth-order deformation equation by Equation (2). It is obvious that when the embedding parameter p = 0 and p = 1, Equation (2) becomes respectively. Then we obtain the mth-order deformation equation for m ≥ 1 with where For simplicity, we select H j (t) = 1 in this problem. So, the approximations of u j (t) are only depend on auxiliary parameter .
Solve the above equation under the initial condition Thus the rest of components u n , n > 4 of the DHAM can be completely obtained. So, we approximate the analytical solution u j (t) = u j,0 (t) + u j,1 (t) + u j,2 (t) + u j,3 (t) + u j,4 (t) 4 + · · · Setting = −1, we get an accurate approximation solution in the following form: is the exact solution of the continuous form.  We can see that different fractional order lead to different diffusion behaviors.
In Figure 2, we show that the method has good agreement with the exact solution when = 1. Example 2. Consider the nonlinear fractional discrete Schrödinger equation with initial condition We can see that different fractional order lead to different diffusion behaviors.
In Figure 2, we show that the method has good agreement with the exact solution when α = 1. We can see that different fractional order lead to different diffusion behaviors.
In Figure 2, we show that the method has good agreement with the exact solution when = 1. Example 2. Consider the nonlinear fractional discrete Schrödinger equation with initial condition (0) = . (20) with initial condition u j (0) = e ijkh . (20) Discrete nonlinear Schrödinger equation is widely used in applied sciences. Describing the propagation of electromagnetic waves in glass fibers, one-dimensional arrays of coupled optical waveguides [18] and light-induced photonic crystal lattices [44]. Moreover, they are an established model for optical pulse propagation in various doped fibers [45,46].
Discrete nonlinear Schrödinger equations are also called lattice nonlinear Schrödinger equations [47]. The parameter ε = h −2 is called (discrete) dispersion and the parameter q is called anharmonicity. Initial value problem (19) and (20) is the discrete form of initial value problem for Schrödinger equation we set u j (t) 2 u j (t) = u 2 j (t)u j (t). By means of DHAM, we choose the linear operator: with property L[c] = 0, where c is constant. We define the nonlinear operator as we construct the zeroth-order deformation equation by Equation (2). For p = 0 and p = 1, we can write respectively. Thus, we obtain the mth-order deformation equation we can select again H j (t) = 1. Thus, the approximations of u j (t) are only depend on auxiliary parameter . Therefore the solution of the mth-order deformation Equation (23) for m ≥ 1 become Entropy 2018, 20, 332
u(x, t) = e i(kx−ωt) , x ∈ R, t > 0, is the plane wave solution of the continuous form, where k is the wave number and ω denotes the frequency. We can observe the different behaviors of the discrete fractional Schrödinger equations, with different fractional parameters.
In Figure 5, we show that the method has good agreement with imaginary part of the exact solution when α = 1.     which is the same solution obtained in [6].     We can observe the different behaviors of the discrete fractional Schrödinger equations, with different fractional parameters.
In Figure 5, we show that the method has good agreement with imaginary part of the exact solution when = 1. In Figure 6, we show that the method has good agreement with real part of the exact solution when = 1.  In Figure 6, we show that the method has good agreement with real part of the exact solution when α = 1. In Figure 6, we show that the method has good agreement with real part of the exact solution when = 1.
with initial condition u j (0) = sin jh.
Initial value problem (27) and (28) is the discrete form of initial value problem for nonlinear fractional Burgers' equation with initial condition u(x, 0) = sin x.
We take into consideration the linear operator: with property L[c] = 0, where c is constant. We can consider the nonlinear operator as Therefore, we construct the zeroth-order deformation equation by Equation (2). For p = 0 and p = 1, we can write ϕ j (t; 0) = u j,0 (t) = u j (0), ϕ j (t; 1) = u j (t), respectively. Thus, we obtain the mth-order deformation equation For simplicity, we select again H j (t) = 1. So, the approximations of u j (t) are only depend on auxiliary parameter .
The solution of the mth-order deformation equation for m ≥ 1 give rise to u j,m (t) = X m u j,m−1 (t) + J α R m → u m−1 .
when we use the initial condition (28) along with (33), we attain the first three of terms of (33) as following: u j,0 (t) = sin jh . . . and so on. Thus, we can conclude that u j (t) = u j,0 (t) + u j,1 (t) + u j,2 (t) + · · ·  We can see that the different behaviors of the discrete fractional Burgers' equations for different fractional parameters.

Discussion and Conclusions
In this paper, the discrete HAM is successfully applied to find the solutions of linear and nonlinear fractional partial differential equations with time derivative (0 < ≤ 1). In contrast to