Numerical Inverse Laplace Transform for Solving a Class of Fractional Differential Equations

Abstract: This paper discusses the applications of numerical inversion of the Laplace transform method based on the Bernstein operational matrix to find the solution to a class of fractional differential equations. By the use of Laplace transform, fractional differential equations are firstly converted to system of algebraic equations then the numerical inverse of a Laplace transform is adopted to find the unknown function in the equation by expanding it in a Bernstein series. The advantages and computational implications of the proposed technique are discussed and verified in some numerical examples by comparing the results with some existing methods. We have also combined our technique to the standard Laplace Adomian decomposition method for solving nonlinear fractional order differential equations. The method is given with error estimation and convergence criterion that exclude the validity of our method.


Introduction
Over the years, researchers have been attracted to study the scientific problems modelled in fractional differential equations due to their constant appearance in the different disciplines of mathematical sciences and engineering such as fluid mechanics, viscoelasticity, mathematical physics, mathematical biology, system identification, control theory, electrochemistry and signal processing [1][2][3][4][5][6].Several analytical and numerical techniques have been developed to solve such kind of equations in the literature.Among these methods, Li and Sun [7] derived the generalized block pulse operational matrix to find the solution of fractional differential equations in terms of block pulse function.A truncated Legendre series together with generalized Legendre operational matrix is used to solve fractional differential equations by Saadatmandi and Dehgan in [8] and they also have presented shifted Legendre-tau method for finding the solution of fractional diffusion equations with variable coefficients in [9].Doha et al. [10] used the shifted Jacobi operational matrix of fractional derivatives applied together with spectral tau-method for solving fractional differential equations.Kazem et al. [11] constructed an orthogonal fractional order Legendre function based on Legendre polynomials to solve fractional differential equations.Mokhtary et al. [12] provided the operational tau method based on Muntz-Legendre polynomials for solving fractional differential equations.Bernoulli wavlet operational matrix of fractional order integration has been derived to approximate the numerical solution of fractional differential equations in [13].Fractional-order Lagrange polynomials have been proposed to solve the fractional differential equations in [14].Albadarneh et al. [15] adopted the fractional finite difference method for solving linear and nonlinear fractional differential equations.
Garrappa [16] provided a detailed survey on the two methods, i.e., product integration rules (PI) and Lubich's fractional linear multi step methods (FLMMs) for solving fractional differential equations.Dehghan et al. [17] adopted homotopy analysis method to solve linear fractional partial differential equations, a meshless approximation strategy for solving fractional partial differential equations based on radial basis function is used in [18].Haar wavlets have been employed to obtain the solutions of boundary value problems for linear fractional partial differential equations by Rehman and Khan [19].Li et al. [20] solved the linear fractional partial differential equations based on operational matrix of fractional Bernstein polynomials.Yang et al. [21] discussed the solution of fractional order diffusion equations within the negative Prabhakar kernel comprise with the Laplace transform and the series solutions in terms of general Mittag-Leffler functions.In [22], a new factorization technique has been adopted for nonlinear differential equations involving local fractional derivatives and found the exact solutions for nonlinear local fractional FitzHugh-Nagumo and Newell-Whitehead equations.Cesarano [23] proposed the Hermite polynomials based operational method to solve fractional diffusive equations.From the last few decades, the Laplace transform method has become popular and adopted by many researchers to solve differential and integral equations.Since then it is necessary to find the inverse Laplace transform for finding the solution in its original domain.There exist a number of analytical and numerical methods for inverting a Laplace transform.For details one can refer [24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40][41].The Laplace transform method reduces the differential or integral equation into a system of algebraic equations due to the Heaviside's operational method.Further, Hasio and Chen [42] extended the idea of Heaviside's operational method to operational matrix of integration.In the literature, few papers reported the use of operational matrices for inversion of Laplace transform: Chen et al. [39] obtained the walsh operational matrices and applied to distributed system, Wu et al. [40] adopted Haar wavlet operational matrices for numerical inverse Laplace transform of certain functions, Aznam and Hussin [38] modified Haar wavlet operational matrices by using generalized block pulse functions, Babolian and Shamloo [24] used operational matrices based on piecewise-constant block pulse function to invert the Laplace transform and solved Volterra integral equations, Maleknejad and Nouri [25] improved the operational matrices based on block pulse functions and Shamloo et al. [41] adopted this method to solve first kind of integral equations.The main purpose of using an operational matrix is that it converts the differential or integral equations into a system of algebraic equations that is simple and easy to solve.
Bernstein polynomials and its operational matrix is used to solve many differential, integral and integro-differential equations [43][44][45][46][47][48][49].But these are not adopted with Laplace transform for inverting the Laplace transform.A Bernstein operational matrix of integration has been developed to find the numerical inverse Laplace transform of certain functions in our paper [50].The proposed method expresses the solution of equations in terms of truncated Bernstein expansion and then using its operational matrix of integration, numerical inverse Laplace transform is obtained.The operational matrix of integration of Bernstein polynomials is easily calculated using a single formula of integration rather than Haar or block pulse function where the order of matrix is taken too large, i.e., 8, 16, 128.Here in our method, we achieve the accuracy using a matrix of order 6 or 7.In the present research, our aim is to find the numerical solutions of linear fractional ordinary and partial differential equations, nonlinear fractional differential equations and system of fractional differential equations using numerical inverse Laplace transform method based on Bernstein operational matrix of integration developed in [50].At first, for linear problems, we transform the fractional differential equations to system of linear algebraic equations using Laplace transform then the numerical approach for calculating the inverse Laplace transform is used to retrieve the time-domain.Here, we extend our numerical approach to solve some nonlinear fractional differential equations together with a well-known iterative method, i.e., a Laplace Adomian decomposition method [51] (briefly explained in the Section 4).One more advantage of using our proposed method is that there is no need of any fractional order matrix of integration for numerical inversion to solve fractional order differential equations.
The paper is organized as follows: In Section 2, some necessary definitions and preliminaries of the fractional calculus theory and Laplace transform are given.Section 3 reveals the basics of Bernstein polynomials, derivation of operational matrix of integration and function approximation in Bernstein polynomials.In Section 4, the proposed method is explained and presented the application to a class of fractional differential equations.Section 5 represents the error estimation and convergence analysis.In Section 6 the illustrative examples are given to show the applicability of the method.Section 7 is refers to conclusions.

Preliminaries
In this section, we give some basic definitions and properties of a Laplace transform [36] and fractional calculus as follows [1,2,4]: where s is known as Laplace variable.
Definition 2. A real function f (t), t > 0 is said to be in the space C µ if µ ∈ R, there exists a real number p > µ and the function f Definition 3. The Riemann-Liouville fractional integral of order α ≥ 0 for a function f (t) is defined as where Γ(.) denotes the Gamma function.
Definition 4. The Riemann-Liouville fractional derivative of order α > 0 for a function f (t) is defined as Definition 5.The Caputo fractional derivative of order α > 0 is defined as where n is an integer, t > 0, and f (t) ∈ C n 1 .
Definition 6.The Laplace transform of a function f (x, t), denoted as F(x, s), t 0 is defined by where s is the transformed parameter.
Definition 7. The Caputo fractional partial derivative operator of order α > 0 is defined as where n is an integer, t > 0.
Property 1.The Laplace transform of the Caputo fractional derivative D α f (t) can be found as

Bernstein Polynomials and Function Approximation
To explain the operational matrices of Bernstein polynomials, we need to give some basic definitions and properties following [52].Definition 8.The Bernstein basis polynomials of degree n are defined over the interval [0, 1]: Some useful results of Bernstein polynomials are as follows: The Bernstein polynomials form a partition of unity i.e., ∑ n i=0 b i,n (t) = 1.Bernstein polynomials have one more important property that these polynomials are not orthogonal therefore to conquer this difficulty Bernstein polynomials of degree n are orthonormailized using Gram-Schmidt orthonormalization procedure and denoted as B i,n (t), i = 0., 1 . . .n.
We can find the function approximation based on Bernstein polynomials.A function f (t), f (t) ∈ L 2 [0, 1] can be expressed in terms of orthonormal Bernstein polynomials [47]: where c i = f , B i,n and ., .denotes the standard inner product.
If the infinite series is truncated at n = k, then the approximate solution can be expressed as where T .Here, the operational matrix of integration of orthonormal Bernstein polynomials is introduced which depends on the integral property of basis vector, i.e., suppose we have a column vector φ(t) = [φ 0 (t), φ 1 (t), . . .φ k (t)] where φ 0 (t), φ 1 (t), . . .φ k (t) are the basis functions orthogonal on some interval [a, b], then the property states that t a . . .
where A m k+1 is the operational matrix of integration of φ(t) which is a constant matrix of order (k + 1) × (k + 1).Now adopting this property on vector B(t), we get where I k+1 is the operational matrix of integration of Bernstein polynomials defined as Therefore

The Method for Numerical Inverse Laplace Transform
In this section, we describe the algorithm proposed in [50] by considering the linear time varying system where u(t) is the unit step function.
Here, we convert this differential equation to integral equation Performing Laplace transform on both sides of (11), we get We can rewrite (12) as Here, we use result from Laplace theory [36] This result can be explained as the integration in time-domain is corresponding to multiplication of 1/s in s − domain.
We can say that here three domains are introduced: one is time-domain, by performing Laplace transform we move from time-domain, f (t) to s − domain F(s), or F 1 s , and then to matrix domain where we define the functional F(I k+1 ) (say) defined on the space of matrices.
The integration in time domain of the inverse Laplace function is corresponding to the definition of the functional F defined onto the space of Bernstein operational matrices.Therefore (13) becomes To solve the integral Equation ( 11), we approximate Also where defined by Therefore, the integral Equation ( 11) becomes Thus, the unknown vector C T is calculated, where F(I k+1 ) can be taken from Equation (15).Consequently, the approximate solution f k (t) can be obtained in terms of Bernstein polynomials by substituting the unknown vector C T into (5).
In all these computations for finding the solution, we used simple algebraic operations of matrices, which have been computed using MATLAB R2014a on Intel R Core TM i3 processor (2328M)(2nd Gen.).

Application to a Class of Fractional Differential Equations
Here, we present the fundamental importance of proposed method by applying it to linear fractional differential equations, linear partial fractional differential equations and nonlinear fractional differential equations.

Application to Linear Fractional Differential Equations
Consider the linear ordinary fractional differential equation with initial conditions Here D α denotes the Caputo fractional order derivative and M( f (t)) represents the linear operator, and may also contain the other derivatives than order α.
To solve the initial problem, Laplace transform is applied to Equation (20) and we attain Therefore, . Hence, f (t), i.e., the solution can be obtained by finding the inverse Laplace transform of H(s) using the above described procedure.

Application to Linear Partial Fractional Differential Equations
Consider the linear partial fractional differential equation of the form with initial condition Taking Laplace transform to both sides of Equation ( 22), we have where F(x, s) and G(x, s) denote the Laplace transform of f (x, t) and g(x, t) respectively.The above equation can be written as Now Equation ( 23) has become a second order ordinary differential equation in F(x, s), that can be easily solved by any classical method for variable x, while keeping s as Laplace variable.The obtained solution F(x, s) can be inverted to f (x, t) using our developed technique.

Application to Nonlinear Fractional Differential Equations
In view to solve, nonlinear fractional order differential equations using standard Laplace adomian decomposition method, we briefly recall the Laplace adomian decomposition method here.
Consider the nonlinear fractional order differential equation with initial condition Here M( f (t)) represents the linear operator which may include other derivatives than order α and N( f (t)) be the nonlinear operator of f (t).The standard Laplace Adomian Decomposition Method (LADM) procedure begins with taking Laplace transform to this nonlinear equation and using the properties of Laplace transform, we have where a = ∑ m−1 i=0 s α−i−1 f (i) (0).The method describes the series solution as Therefore, truncated series solution takes the form and the nonlinear term is decomposed as follows: where A n 's are the Adomian polynomials, given by Consequently, Equation ( 25) leads to following recurrence relation In general The new development is that we find the values of f 0 (t), f 1 (t) . . .by finding inverse Laplace transform using our proposed technique, i.e., Bernstein operational matrix, as described above.

Error Estimation via RK45 Method
We have investigated error estimation of the proposed method [24].The error function e k (t) of the truncated Bernstein expansion f k (t) is defined as From the following theorem the absolute error bound can be estimated: Theorem 1.Let f (t) be the function defined on [0, 1], then the upper bound for the errors in truncated Bernstein expansion can be estimated.
Proof.Suppose that f k (t) and f k+1 (t) are two consecutive approximate solutions of given differential equation.If the error sequence is increasing or decreasing, then we use the triangle inequality to the approximate solutions at k and k + 1 and we achieve where . Hence, it can be said that the error can be bounded from above.One of the absolute errors e k (t) or e k+1 (t) are bounded by f k+1 (t) − f k (t) ∞ , if the error sequence is monotone.The computable error bound is represented for 0 ≤ β < 1.But, the solution diverges when the Bernstein series diverges for β > 1.

Convergence Analysis
In order to show the effectiveness of our method, a residual function of a linear time varying system in the Banach space for the values of k is adopted to interpret the convergence of the Bernstein polynomials solution as described in [47,53].Suppose f k (t) are the approximate solution of (10), we write the residual function as The Bernstein polynomial-based numerical solution or the residual function can be expressed in terms of Taylor series expansion as: Now, we desire to prove that the residual function sequence is convergent in Banach space and satisfying the condition: For the convergence, we prove here that {R k+1 (t)}, k = 0, 1, . . .∞ is a Cauchy sequence.
Let us take that Therefore, the given condition can be written as and this can be easily shown as follows.Let us write explicitly as functions of Bernstein polynomials If we define we have Here we used the following property of Bernstein polynomials so that From where, we get Since the function in brackets is a decreasing function in Analogously, by using the properties of Bernstein polynomials, we can also estimate and since Bernstein polynomials are partition of unit there follows Let us combine Equation (36) with the condition in Equations ( 42) and ( 43), then we get and analogously This can be written as Here lim Thus, we generalize the inequality and attain For k, s ∈ N and k ≥ s, to prove that the sequence is Cauchy sequence, we take where that proves the residual function sequence is Cauchy sequence and convergent.Therefore the approximate solution is convergent.

Illustrative Examples
Here, we present some examples to demonstrate the applicability of the presented technique.The relative errors and L ∞ -errors are plotted at the different values of k.Also the error estimation is discussed in each example.
Example 1.Consider the following fractional differential equation [12] with the exact solution f (t) = t 2 √ t.
The relative error for different values of k are plotted in Figure 1 which shows that the method is giving good results as compared to analytic solution.Using the error estimation method, we estimated the upper bound of errors and calculate Therefore, it is clear from the data that the error is estimated and bounded above by f 6 − f 7 ∞ .Example 2. Consider the following fractional differential equation with exact solution f (t) = t 3 .
The relative errors for different values of k are plotted in Figure 2. Using the error estimation method, we estimated the upper bound of errors and calculated with the exact solution f (t) = t 2 .
In Figure 3, we plotted the relative errors for k = 5 and k = 6.To show the efficiency, absolute errors compared to fractional finite difference method (FFDM) [15] are given in Table 1 which clearly states that our method gives more accurate results.Using the error estimation method, we also estimate the upper bound of errors and calculate e 6 ∞ = 8.17 × 10 −4 , e 7 ∞ = 2.97 × 10 −5 and f 6 − f 7 ∞ = 8.17 × 10 −4 .Therefore max{ e 6 ∞ , e 7 ∞ } ≤ f 6 − f 7 ∞ that follows the estimated upper bound of error.Example 4. Consider the following mathematical model, which is developed for a micro-electro mechanical system instrument, that has been designed to measure the viscosity of the fluids that are encountered during oil exploration [54]: In Figure 4, we plotted the relative errors for k = 5 and k = 6, taking β = 1.To show the efficiency, absolute errors compared to cubic spline method [54] are given in Table 2 which clearly states that our method gives more accurate results than the method in [54].Using the error estimation method, we also estimate the upper bound of errors and calculate e 6 ∞ = 2.61 × 10 −4 , e 7 ∞ = 7.43 × 10 −5 and f 6 − f 7 ∞ = 1.90 × 10 −4 .Therefore e 7 ∞ ≤ f 6 − f 7 ∞ that follows the estimated upper bound of error.
subject to the conditions where the exact solution of system at α = β = 1 is f 1 (t) = e t sin t and f 2 (t) = e t cos t.
The relative error for different values of k are plotted in Figures 5 and 6.Using the error estimation method, we estimate the upper bound of errors for f 1 (t) and calculate e 6 ∞ = 1.59 × 10 −3 e 7 ∞ = 6.18 × 10 −5 We observe that the result max{ e 6 ∞ , e 7 ∞ } ≤ f 6 − f 7 ∞ is satisfied here.Similarly for the function f 2 (t) the results are reported as The numbers clearly show that the error is bounded from above, i.e., max{ e 6 ∞ , e 7 ∞ } ≤ f 6 − f 7 ∞ .To show the efficiency, absolute errors compared to variation iteration method (VIM) [55]  are given in Table 3 and 4 which excludes that our method gives more accurate results than the method in [55].
The exact solution is not known.
Here, we first convert the fractional partial differential equation to ordinary differential equation as described in Section 4, then the solution is obtained by our proposed method for α = 1.5, k = 6 presented in Table 5 and have compared with the method VIM (Table 6) and Inverse Multiquadric-Radial Basis Functions (IMQ-RBF) (Table 7) which shows that the solution is in good agreement with VIM and IMQ-RBF [18].
with exact solution f (t) = t 3 .
The relative error for different values of k are plotted in Figure 7 which show that our method gives a very close solution to the analytic solution.Using the error estimation method, we estimated the upper bound of errors and calculated  Hence, we can see that here f 6 − f 7 ∞ bounds the error e 6 ∞ .We also analyzed the L ∞ -error at different values of k, i.e., k = 5, 6, 7, 8 for examples which clearly exclude that we achieve the best results at k = 6, as shown in Figures 8-11

Conclusions
Enormous efforts and advances have been conducted to obtain the numerical solutions of fractional differential equations.An operational matrix of orthonormal Bernstein polynomials is derived to find the inverse Laplace transform in [50].Here, the practical use of our proposed method is discussed for finding the solutions of some fractional order ordinary differential equations including the mathematical model of instrument (MEMS) and partial differential equations (particularly wave equation) that converts the problem to system of linear algebraic equations.The accuracy of the method is illustrated while comparing the solutions with some existing methods like VIM, FFDM, cubic spline and IMQ-RBF.We have also combined our method with Laplace adomian decomposition method that is advantageous to solve nonlinear fractional differential equations.Finally, we have analyzed the solution of each illustrative example at different values of k, i.e., at k = 5, 6, 7, 8 and with the help of graphs relative errors are shown.It is also observed from the plots of supremum norm error that at k = 6, we achieve the best accurate result compared to others.

Figure 8 .
Figure 8.The L ∞ − error for Example 1 (left) and Example 2 (right) for different values of k.

Figure 9 .
Figure 9.The L ∞ − error for Example 3 (left) and Example 4 (right) for different values of k.

Figure 10 .
Figure 10.The L ∞ − error for f 1 (t) (left) and f 2 (t) (right) at different values of k for Example 5.

Figure 11 .
Figure 11.The L ∞ − error for at different values of k for Example 7.

Table 1 .
Comparison of absolute error in Example 3.

Table 2 .
Comparison of Absolute error in Example 4.

Table 3 .
Comparison of Absolute errors in f 1 (t) in Example 5.

Table 4 .
Comparison of Absolute errors in f 2 (t) in Example 5.

Table 5 .
Solution by proposed method for different values of x and t in Example 6.

Table 6 .
Solution by variation iteration method (VIM) for different values of x and t in Example 6.

Table 7 .
Solution by IMQ-RBF for different values of x and t in Example 6.