Approximation of linear one dimensional partial differential equations including fractional derivative with non-singular kernel

In this article we propose a hybrid method based on a local meshless method and the Laplace transform for approximating the solution of linear one dimensional partial differential equations in the sense of the Caputo–Fabrizio fractional derivative. In our numerical scheme the Laplace transform is used to avoid the time stepping procedure, and the local meshless method is used to produce sparse differentiation matrices and avoid the ill conditioning issues resulting in global meshless methods. Our numerical method comprises three steps. In the first step we transform the given equation to an equivalent time independent equation. Secondly the reduced equation is solved via a local meshless method. Finally, the solution of the original equation is obtained via the inverse Laplace transform by representing it as a contour integral in the complex left half plane. The contour integral is then approximated using the trapezoidal rule. The stability and convergence of the method are discussed. The efficiency, efficacy, and accuracy of the proposed method are assessed using four different problems. Numerical approximations of these problems are obtained and validated against exact solutions. The obtained results show that the proposed method can solve such types of problems efficiently.


Introduction
The integer order derivatives are local in nature, i.e., using these derivatives changes in a neighborhood of a point can be described, but using the arbitrary order derivatives we can describe changes in an interval. Namely, fractional derivatives are nonlocal in nature. That is why arbitrary order derivatives are suitable to model more physical phenomena such as electrodynamics, earthquake vibrations, diffusion process, polymers, fluid flow, elasticity, hydrology, and signal and image processing [1][2][3][4][5][6][7]. Fractional order derivative of a function is in fact its definite integral. A geometrically fractional derivative accumulates the function. The corresponding accumulation includes the integer order derivative as a special case. In the literature a large number of arbitrary order derivatives are available such as the Marchaud, Grünwald-Letnikov, Riemann-Liouville, and Caputo derivatives [8][9][10][11][12][13][14][15]. However, these derivatives have certain disadvantages, these derivatives contain kernels with singularities. It is reported in the literature that these derivatives face problems when someone tries to model nonlocal phenomenon [1,16]. Therefore, to model the non-local phenomenon accurately, recently a new fractional derivative has been defined, called the Caputo-Fabrizio (CF) derivative [17].
The research community took great interest in the CF derivative due to its non-singular kernel. The CF derivative has a large number of applications, and researchers have applied it successfully to many phenomena like signal processing [18], ground water pollution [19], groundwater flow in a confined aquifer [1,20], the mass-spring damper system [21], non-Darcian flow and solute transport [22].
In many real world problems as regards their modeling, we need to obtain their exact solution which is some times quite difficult specially for nonlinear problems. Therefore, we need some sophisticated tools to deal with such problems. The researchers increasingly have used numerical and analytical techniques to handle the problems for corresponding numerical and analytical solutions. For example a homotopy analysis method based on the Laplace transform is proposed for solving linear problems including the CF derivative in [2,23]. In [1] the authors solved the model of groundwater flow with the CF derivative using the Sumudu transform. In [24] a fractional order Fisher diffusion equation with the CF derivative is solved via some iterative method. The Crank-Nicholson scheme in [25] is applied for approximating the solution of the Allen Cahn model with CF derivative. A numerical approximation of the solution of groundwater pollution model with spacetime fractional CF derivative in [19] is obtained via the Crank-Nicholson scheme. The Laplace and Fourier transforms are utilized in [26] for obtaining the fundamental solution of advection-diffusion equation with the CF derivative. In [27] the authors applied the CF derivative to the analysis of a rock fracture process. Many other numerical and analytical techniques for the solutions of different fractional order models can be found in [28][29][30][31][32][33][34][35][36][37][38][39] and the references therein. In this work our aim is to approximate linear PDEs with the CF derivative via the Laplace transform (LT) and local meshless method. The Laplace transform is coupled with a local meshless method for avoiding the time stepping procedure. The Laplace transform will help us in obtaining the solution in less computation time and without time instability. The paper is organized as follows: In Sect. 2 some basic definitions are given; In Sect. 3 the Laplace transform and contour integration methods are discussed; Sect. 3.1 shows how to construct the differentiation matrices via a local meshless method; Sect. 4 discusses the convergence of the method; Sect. 5 discusses the stability of the method; Sect. 6 contains the numerical tests, where theory and experiment are compared.

Preliminaries
Definition 2.1 The LT of a piecewise continuous function f (t), t > 0 is defined as ( 1 )

Definition 2.3
If n ∈ N, and 0 ≤ α ≤ 1 then the LT of the CF fractional derivative is defined as [2,24] using n = 0 we get and for n = 1 we have

Numerical scheme
In order to validate our method let us consider a linear one dimensional partial differential equation with the CF derivative for q -1 < α + n ≤ q: the boundary and initial conditions are and where L is the linear differential operator and B is the boundary operator. Applying the LT to Eq. (6) and Eq. (7), we have and Bû(x, s) =ĝ(x, s).
From Eq. (8), we have on simplification we get Bû(x, s) =ĝ(x, s), (12) whereĜ(x, s) iŝ We obtain the solution u(x, t) of (6)-(7) by representing it as an integral along a smooth curve as where σ 0 ∈ R is the convergence abscissa, it is large enough and the contour is a suitably selected line 0 perpendicular to x axis connecting σ -ι∞ and σ + ι∞. Then in Eq. (13), u(x, t) is the inverse Laplace ofû(x, s), with the condition that the 0 lie to the right of all the singularities of the transformû(x, s). For our purposes, however, we assume that u(x, s) may be continued analytically in an appropriate way. We shall want to take for 0 a deformed contour in the set δ φ = {0} ∪ {s = 0 : | arg s| < φ}. The deformed contour will show asymptotic behavior like the couple of lines in the left half complex plane with Im s → ±∞ and Re s → -∞, forcing the factor e st to decay in the direction of both ends of the contour . In this work we select with the parametric representation of the form [40] where Letting s = x + ιy, we notice that Eq. (14) serves as the left part of the hyperbola defined by where for Eq. (16) the asymptotes are y = ±(xδξ ) cot θ , and its x-intercept is s = δ + ξ (1 -sin θ ). Equation (15) ensures that the contour lies in δ φ = δ + φ ⊂ φ , and extends towards the left half complex plane. Using Eq. (14) in Eq. (13) we have The approximation of Eq. (17) can be obtained by employing the trapezoidal rule as The solution u k (x, t), can be obtained by solving the system of equations in (11)- (12) in parallel for the nodes s j . In order to do this we employ the meshless method in local setting for discretization of the operators L and B.

Local meshless scheme
where the vector λ i = {λ i j } n j=1 represents the expansion coefficients, φ(r) is kernel function and r = x ix j is the distance between x i and x j . i and are local and global domains, respectively. The local domain i includes the center x i and its n neighboring centers around it. Hence we have N n × n linear systems as which can be written aŝ where i is called the system matrix having entries of the form Solving the N systems in Eq. (21) we obtain the unknowns λ i = {λ i j } n j=1 . The differential operator L can be approximated as Equation (22) can be written as where the vectors ν i and λ i have orders 1 × n and n × 1, respectively. The entries of ν i are of the form solving Eq. (21) and Eq. (23) for λ i , we have hence at each center x i the operator L can be approximated via a local meshless method as where D N×N is a sparse system matrix obtained via meshless method in local setting approximating L. We can approximate the operator B in a similar way.

Accuracy and convergence
The process of approximating the solution of the problem defined in Eq. (1)-Eq. (3) involves the transformation of the given equation by employing the Laplace transform, no error occurs in this process. Then we employ the meshless method in a local setting for solving the transformed problem. The error estimate is O(γ 1 h ), 0 < γ < 1 (where h is the fill distance and is the shape parameter) for the meshless method in the local setting. Finally, we obtain the solution of the problem via the inverse Laplace transform by representing it as a Bromwhich integral (17). The integral is then approximated to high accuracy via the trapezoidal rule. While approximating the integral (17) the convergence rate is dependent on . Also in this process the convergence orders are dependent on the step k and on the chosen temporal domain [t 0 , T]. For optimal results and best convergence we choose an optimal temporal domain. A proof of the error is given in the next theorem.

Stability of the method
For the stability of the system defined in Eq. (11)-Eq. (12), we write the system in matrix form given as where Q N×N is obtained via a local meshless method. The constant of stability for the system Eq. (28) is defined as for any discrete norm · defined on R N the constant S c is finite. From Eq. (29) we may write Thus, we have From Eq. (30) and Eq. (32) we observe that the constant S c is bounded. Calculation of the pseudoinverse for the system (28) may be computationally expensive, but it confirms the stability. For square matrices we can use the MATLAB's function condest for estimating Q -1 ∞ , hence we have For our sparse matrix Q this works well in a very short time of computation (Fig.1).

Numerical experiments
Here we employ our proposed numerical scheme for the approximations different one dimensional linear PDEs. For optimization of the shape parameter we have utilized the uncertainty principle due to [41]. In our numerical experiments for generating the quadrature nodes the command μ = -M : k : M is used, also the multiquadrics (MQ) kernel is used in all experiments. Other optimal parameters involved for the contour of integration are θ = 0.15410, r = 0.13870, η = 0.10, τ = t0 T , δ = 2.0, t ∈ [t 0 , T] = [0. 5,5]. The accuracy of the method is measured via the L ∞ error given by Here u(x, t) k and u(x, t) are the numerical and exact solutions respectively.  Example 1 Here we consider a 1-D linear differential equation with the CF derivative [2] CF 0 D α+1 t u(x, t) - the initial conditions are The boundary conditions are extracted from the exact solution The results obtained for different centers n in local domain i and N in global domain are given in Table 1 for fractional order α = 0.22211 and various quadrature points. The condition number κ, L ∞ errors, shape parameter ε, and the error estimates are also given in Table 1. From Table 1 it is clear that the proposed method gives an acceptable accuracy. The approximate and exact solutions are plotted in Fig. 2(a) and Fig. 2(b), respectively, which clearly shows good agreement between them. The contour plot of absolute error is shown in Fig. 3(a), which shows a high accuracy. The plot of actual error vs error estimate  is shown in Fig. 3(b), which shows good agreement between them. From the results we conclude that the proposed method is effective.
Example 2 Here we consider a non-homogeneous linear differential equation with the CF derivative [2]: the initial condition being the boundary conditions can be obtained from the exact solution The results obtained for different centers n in local domain i and N in global domain are shown in Table 2 with fractional order α = 0.22211 and various quadrature nodes. From the observations on this table also we conclude that the results obtained using the proposed numerical scheme are acceptable. The plots of exact and approximate solutions are given in Fig. 4(a), and Fig. 4(b), respectively. The error functions for different values of α are shown in Fig. 5(a), we observe that the error decreases as we increase α. In Fig. 5 Example 3 Next the 1-D non-homogeneous diffusion equation with the CF derivative is considered:

(b)
with the forcing term The problem has the following exact solution: u(x, t) = sin(2πx) t 3 + 1 .   The initial and boundary conditions are extracted from the exact solution. The results obtained for different nodes n ∈ i , N ∈ and different quadrature nodes with α = 0.3, α = 0.7 are shown Table 3 and Table 4. A clear improvement is observed in the accuracy as we increase the spatial and quadrature nodes. The plots of exact and numerical solutions are given in Fig. 6(a) and in Fig. 6(b). In Fig. 7(a) the contour plot of the absolute error is given, which shows a high accuracy of the proposed numerical scheme. In Fig. 7(b) the plot of the actual error and error estimate is shown, which shows good agreement between them. The results ensure the convergence, stability, and efficiency of our method. Example 4 The next test problem is the 1-D non-homogeneous diffusion equation with the CF derivative: The initial condition is with the forcing term The exact solution is The initial and boundary conditions are extracted from the exact solution. The results obtained for different nodes n ∈ i , N ∈ and different quadrature nodes with α = 0.7 are given Table 5. The plots of exact and numerical solutions are given in Fig. 8(a) and in Fig. 8(b). In Fig. 9(a) the contour plot of absolute error is shown, and in Fig. 9(b) the plot of actual error vs error estimate is shown, which shows good agreement between them. For this problem also the method produced accurate and stable results.

Conclusion
In this work, we have coupled the Laplace transform and a meshless method in a local setting successfully for the approximation of linear one dimensional partial differential equations with the CF derivative. The Laplace transform reduced the problem to a time independent problem and in this way the classical time stepping procedure and hence the   time instability issues were avoided. The local meshless method produced sparse differentiation matrices and the issues of ill conditioning and shape parameter sensitivity were resolved. We discussed the convergence and stability of the method. Some numerical experiments were performed on one dimensional problems with a CF derivative for different fractional orders. Numerical results confirmed the stability and high accuracy of the proposed method. The obtained results ensured the efficiency and applicability of the method for such problems.