Numerical solution of fractional partial differential equations by numerical Laplace inversion technique

In this paper, we propose a numerical method for solving fractional partial differential equations. This method is based on the homotopy perturbation method and Laplace transform. The transformed problem obtained by means of temporal Laplace transform is solved by the homotopy perturbation method. Then we use Stehfest’s numerical algorithm for calculating inverse Laplace transform to retrieve the time domain solution. The approximate solutions obtained by our proposed method are in excellent agreement with the exact solutions. It is worthwhile to note that our method is applicable to a variety of fractional partial differential equations occurring in fluid mechanics, signal processing, system identification, control robotics, etc. The utility of the method is shown by solving some interesting examples. MSC:34A08, 44A10.


Introduction
Fractional differential equations are found to be an effective tool to describe certain physical phenomena such as damping laws, rheology, diffusion processes, and so on. Several methods have been developed to solve fractional differential equations. Lin and Xu [] proposed the numerical solution for a time-fractional diffusion equation. In [], an unconditionally stable finite element (FEM) approach for solving a one-dimensional Caputo-type fractional differential equation with singularity at the boundary was presented. Kexue and Jigen [] discussed the Laplace transform (LT) method for solving fractional differential equations with constant coefficients. Jafari et al. [] applied the homotopy analysis method to obtain the solution of a multi-order fractional differential equation in the Caputo sense. Merrikh-Bayat [] developed a low-cost numerical algorithm to find the series solution of nonlinear fractional differential equations with delay. In [], the Riemann-Liouville fractional integral for repeated fractional integration was expanded in block pulse functions to yield the block pulse operational matrices for the fractional order integration. Esmaeili et al. [] developed a computational technique based on the collocation method and Muntz polynomials for the solution of fractional differential equations. In [], three different numerical methods were used to solve a singularly perturbed Able Volterra integral equation, presented by a fractional differential equation. Ibrahim [] discussed holomorphic solutions for nonlinear singular fractional differential equations. http://www.advancesindifferenceequations.com/content/2013/1/375 Homotopy perturbation method (HPM) has been applied by several researchers to solve different kinds of functional equations. This method was further developed and improved by He [] and applied to develop a coupling method for a homotopy technique [ The Laplace transform method has been applied to a wide class of ordinary differential equations (ODEs), partial differential equations (PDEs), integral equations (IEs) and integro-differential equations (IDEs). In these problems it is necessary to calculate the Laplace transform and inverse Laplace transform of certain functions. The inverse of Laplace transform is usually difficult to compute by using the techniques of complex analysis, and there exist numerous numerical methods for its evaluation [, ]. Sastre et al. [] developed an application of Laguerre matrix polynomial series to the numerical inversion of Laplace transforms of matrix functions. Laguerre matrix polynomials were introduced in [] and theorems for the expansion of matrix functions in series of Laguerre matrix polynomials can be found in [, ]. In [], the dynamical differential equations with initial conditions were converted into the model of linear operator action, in which the linear operator is just the infinitesimal generator for the solver of the differential equations, and the resolvent of the linear operator is the Laplace transform of the solver of original differential equations. In [], a method for the numerical inversion of Laplace transform on the real line of heavytailed (probability) density functions is presented. The method assumes a finite set of real values of the Laplace transform and chooses the analytical form of the approximant maximizing Shannon-entropy, so that positivity of the approximant itself is guaranteed. In [], a Laplace homotopy perturbation method is employed for solving one-dimensional non-homogeneous partial differential equations with a variable coefficient. This method is a combination of the Laplace transform and the homotopy perturbation method (LHPM). LHPM presents an accurate methodology to solve non-homogeneous partial differential equations with a variable coefficient. Sheng et al. [] proposed an application of numerical inverse Laplace transform algorithms and obtained an easy way to solve the complicated fractional-order differential equations numerically. [] proposed a new algorithm for the numerical inversion of Laplace transforms by using multi-precision computational environment and provided controlled accuracy, that is, the inversion can be carried out to yield any pre-specified number of significant digits. The fundamental collocation method was extended to handle two-dimensional transient heat conduction problems in solids in []. The method was applied in the Laplace transform domain, followed by an inversion technique to retrieve the time-domain solution.
In [], the authors developed a numerical algorithm for inverting a Laplace transform, http://www.advancesindifferenceequations.com/content/2013/1/375 based on Laguerre polynomial series expansion of the inverse function under the assumption that the Laplace transform is known on the real axis only. The main contribution of the paper is to provide computable estimates of truncation, discretization, conditioning and roundoff errors introduced by numerical computations. In the present work, we apply the Stehfest [] algorithm for numerical inversion of Laplace transform.
In this paper, the method for numerical solution of fractional partial differential equations is based on Laplace transform (LT), the homotopy perturbation method (HPM) and Stehfest's numerical algorithm for calculating inverse Laplace transform. The accuracy and efficiency of the method is verified by solving some examples of physical interest.

Homotopy perturbation technique
In this section, we describe the homotopy perturbation method [-] for a general type of the nonlinear differential equation with boundary conditions where A is a general differential operator, B is a boundary operator, f (r) is a known analytical function and is the boundary of the domain . The operator A can be divided into two parts L and N , where L is a linear operator and N is a nonlinear operator. Therefore, Eq. () can be rewritten as follows: By the homotopy technique, we define a homotopy H(r, p) : × [, ] → R as follows: or where p ∈ [, ] is an embedding parameter, and u  is an initial approximation for Eq. () with Note that the process of varying the values of p from zero to unity corresponds to that of u(r, p) from u  (r) to u(r). We assume that the solution of Eq. () can be written as a power series in p, that is, Substituting () in () and comparing the coefficients of powers of p yields a successive procedure to determine u k . Finally, by setting p =  in (), we obtain the solution of Eq. (). http://www.advancesindifferenceequations.com/content/2013/1/375

Preliminaries
In this section, we recall some basic concepts of fractional calculus [-] and Laplace transform. () Definition  A two-parameter Mittag-Leffler function is defined by the following series: where s is the transform parameter and is assumed to be real and positive. Note that the Laplace transform of Mittag-Leffler function E α,β (t) is .

Description of the method
Consider the following linear fractional partial differential equation: with the initial conditions and the boundary conditions where f k , k = , , . . . , m-, h, g  , g  , A and B are known functions and T >  is a real number and m - < α ≤ m. Now we explain the method of solution for solving initial-boundary value problem ()-().
Taking the Laplace transform of problem ()-() and using (), we obtain where (x, s) and h(x, s) denote the Laplace transform of u(x, t) and h(x, t), respectively, and Rewriting Eq. (), we have h(x, s).
In the limit p → , we note that () becomes the approximate solution for the problem of ()-() and is given by

Applying Stehfest's algorithm [] to H n (x, s), the solution u(x, t) is found to be
where p is a positive integer and Here [r] denotes the integer part of the real number r.

Numerical results
In this section, we show the efficiency and accuracy of the new Laplace homotopy perturbation method (LHPM) by applying it to several test problems.

Example  Consider the following initial-boundary value problem []
.

Conclusions
In this paper, we have developed a new numerical method for solving fractional partial differential equations. This method is based on Laplace transform, the homotopy perturbation method and Stehfest's numerical algorithm for calculating inverse Laplace transform. We demonstrate the efficiency and accuracy of the proposed method by applying it to three typical examples. It is found that the approximate solutions produced by our method are in complete agreement with the corresponding exact solutions. Moreover, in view of its simplicity, our method is applicable to a wide class of initial-boundary value problems occurring in applied sciences.