A Meshless Method Based on the Laplace Transform for the 2D Multi-Term Time Fractional Partial Integro-Differential Equation

: In this article, we propose a localized transform based meshless method for approximating the solution of the 2D multi-term partial integro-differential equation involving the time fractional derivative in Caputo’s sense with a weakly singular kernel. The purpose of coupling the localized meshless method with the Laplace transform is to avoid the time stepping procedure by eliminating the time variable. Then, we utilize the local meshless method for spatial discretization. The solution of the original problem is obtained as a contour integral in the complex plane. In the literature, numerous contours are available; in our work, we will use the recently introduced improved Talbot contour. We approximate the contour integral using the midpoint rule. The bounds of stability for the differentiation matrix of the scheme are derived, and the convergence is discussed. The accuracy, efﬁciency, and stability of the scheme are validated by numerical experiments.


Introduction
Recently, the theory fractional calculus has gained significant attention in the field of engineering and other sciences because of its various applications in modeling numerous phenomena. For example numerous phenomena in the mathematical biology, physics, and engineering fields can be described by fractional integro-differential equations (FIDEs), fractional partial integro-differential equations (FPIDEs), and fractional partial differential equations (FPDEs). In particular, several phenomena give rise to fractional partial integro-differential equations such as viscoelastic phenomena, signal processing, and fluid mechanics (see [1][2][3][4] and the references therein).
Finding exact or numerical solutions of FPIDEs with weakly singular kernels is an important task. Due to the possible singularities of the kernel function at the origin [5], sharp changes will occur in the solution, so the exact/analytic solution may be difficult to obtain [6]. Therefore, the alternate way is to develop an accurate numerical scheme. The approximation of FPIDEs and partial integro-differential equations (PIDEs) has been considered by many researchers, such as the authors in [7][8][9][10], who used the finite difference (FDM) and finite element (FEM) methods for the approximation of PIDEs. The authors in [11,12] used spline collocation methods for approximating PIDEs of the parabolic and hyperbolic type, respectively. Huang [13] solved parabolic PIDEs using the time discretization scheme. In [3], the B-spline solution of FPIDEs was found. The authors in [4] used the reproducing kernel method for the approximation of FPIDEs. A backward Euler difference scheme was constructed for the approximation of the partial integro-differential equation with multi-term kernels [14]. Other valuable work on integro-differential equations can be found in [15][16][17][18][19][20][21][22][23][24][25][26] and the references therein. Recently, the meshless methods have attracted researchers and become the primary tool for interpolating multidimensional scattered data.
In the literature, we can find a large number of meshless methods developed for the numerical treatment of different PDEs or PIDEs such as the authors in [27], who proposed the RBF-FD method for the approximation of PIDEs. In [28], a local method was developed for the approximation of PIDEs. In [29], the nonlinear PIDEs were approximated via the RBF and theta method. Similarly, the authors in [30] proposed a local method with the optimal shape parameter for PIDEs.
In this article, the Laplace transform and localized meshless method are combined for the approximation of the solution of FPIDEs with a weakly singular kernel. In the literature, we can find some valuable work on the Laplace transform coupled with other methods in [31][32][33][34][35][36][37][38][39] and the references therein. We consider an FPIDE of the form [40]: subject to the initial and boundary conditions: Θ is a given function, L = ∆, and B is the boundary operators. D γ τ is the Caputo derivative of order γ defined as [41,42]: in particular for p = 1, we have: and the αth integral (α) U(τ) is defined as [42]: The Laplace transform of U(τ) is defined by: and the Laplace transform of D γ τ is given by: while the Laplace transform of (α) U(τ) is given by: Equation (1) involves multi-term time fractional derivatives, which are helpful in modeling complex physical systems such as the physical multi-rate phenomenon [43], fractional Zener model [44], heavily damped motion, Newtonian fluid [42,45], and anomalous relaxation process [46]. The integral term represents the viscosity part of Equation (1) [40]. Equation (1) can be obtained from the partial integro-differential equation: by replacing the first order time derivative by a linear combination of fractional derivatives of different orders [47]. For m = 1, Equation (1) becomes the single term FPIDE of the form: and for α = 1, Equation (1) becomes the fractional multi-term diffusion equation:

Proposed Numerical Scheme
Taking the Laplace transform of Equations (1) and (2), we have: B{ U(ζ, s)} = h(ζ, s); combining the like terms, we obtain the following system: where: where I is the identity operator and L and B are the governing and the boundary differential operators. In order to solve the system given in (11) and (12), first, we employ the local meshless method to discretize the operators L and B. When we are done with the discretization of these two operators, the system of Equations (11) and (12) is solved in parallel for each point s along some suitable path Γ in the complex plane ( see, e.g., [16,34]). Finally, the solution of the original problem (1) and (2) is obtained using the inverse Laplace transform. The local meshless method for the discretization of the given differential operators L and B is described in the next section.

Localized Meshless Method
In the localized meshless method for a given set of nodes {ζ i } N i=1 ⊂ Ω, the local meshless approximate of the function U(ζ) has the form [48]: where φ(r) is a radial kernel, r = ζ i − ζ h is the distance between ζ i and ζ h , λ i = {λ i h } n h=1 is the vector of expansion coefficients, Ω is the global domain, and Ω i is a local domain containing ζ i , and n nodes around it. Hence, we obtain N linear systems of order n × n given by: the elements of the system matrix are found by solving each n × n system. Next, the operator L U(ζ) is approximated by: The above Equation (15) can be expressed as: where λ i is an n-column vector and ν i is an n-row vector with entries: solving Equation (14), for λ i , we have, From Equation (18), we use λ i in Equation (16) and get, where, Thus, the local meshless approximation for the differential operator L at each center ζ i is given as: The differentiation matrix D is sparse and has order N × N containing n number of non-zero entries, where n ∈ Ω i . The matrix D approximates the linear differential operator L. The approximation for the boundary differential operator B can be done in the same way.

Numerical Inversion of the Laplace Transform
Following the discretization by the local meshless method of the linear differential and boundary operators L and B, respectively, the system (11) and (12) is solved in parallel for each point s along some suitably chosen path in the complex plane. Finally, we get the solution of the problem (1) and (2) using the inversion formula: where σ 0 ∈ R is called the converging abscissa and Γ is an initially appropriately chosen line connecting σ − i∞ to σ + i∞. This means all the singularities of U(ζ, s) lie in the half plane Res < σ. The approximation of the integral (22) is hard because of the slow decaying transform U(ζ, s) and the highly oscillatory exponential factor e sτ . To handle these issues, we use the strategy suggested by Talbot [49]. He suggested the deformation of the contour of integration Γ. In particular, he suggested that the contour Γ be deformed in such a way that its real part starts and ends in the left half plane, and it encloses all the singularities of the transform U(ζ, s). Cauchy's theorem allows such a deformation provided that U(ζ, s) has no singularities on the contour [49]. On such contours, the exponential factor decays rapidly, which makes the integral in Equation (22) suitable for approximation using the midpoint or trapezoidal rule [49][50][51]. We consider a Hankel contour with the parametric form given by [50]: where Res(±π) = −∞, and s(ϑ) is defined as: where the parameters σ, γ, µ, and ν are to be described by the user. From Equations (22) and (24), we have: We use the M-panel midpoint rule with uniform spacing k = 2π M to approximate the integral in Equation (25) as:

Convergence and Accuracy
In order to approximate the solution of FPIDEs using our proposed numerical scheme, the Laplace transform and local meshless method are used. In our numerical scheme, we employ the Laplace transform to the time dependent equation, which eliminates the time variable, and this process causes no error. Then, the local meshless method is utilized for approximating the time independent equation.
The error estimate for local meshless method is of order O(η 1 h ); 0 < η < 1; is the shape parameter; and h is the fill distance [52]. In the process of approximating the integral in Equation (25), we achieve the convergence at different rates depending on the path Γ. While approximating the integral in Equation (25), the convergence order relies on the step k of the quadrature rule and the time domain [t 0 , T] for Γ. In order to achieve high accuracy, we need to search for the most favorable values of the parameters involved in Equation (24). The authors [50] obtained the most suitable values of the parameters given as: with the corresponding error estimate as: The optimal contour is supposed to pass neither too close to nor too far from the singularities, in which case e sτ becomes too large. Moreover, to achieve the desired accuracy, the quadrature points are required to extend far enough into the left half plane. However, their contributions must not be less than the required accuracy.
The necessary steps of our method are presented in the following Algorithm 1.

Algorithm 1:
1: Input: The computational domain, the fractional order derivative in [0,1], the final time, the contour of integration, the initial shape parameter and the other parameter of the given model, the inhomogeneous function, and other conditions. 2: Step 1: Apply the Laplace transform to the problem (1) and (2), and obtain the time independent problem (11) and (12). 3: Step 2: Discretize the linear differential operator L and boundary operator B using (21). 4: Step 3: Solve the system of Equations (11) and (12) in parallel for each point s along the contour of integration Γ given in (18). 5: Step 4: Compute the approximate solution using (25). 6: Output: The approximate solution is U k (ζ, τ).

Stability Analysis
In order to investigate the system (11) and (12) stability, we represent the system in discrete form as: where Y N×N is the sparse differentiation matrix obtained using the local meshless method described in Section 2.1. For the system (27), the constant of stability is given by: where the constant C is finite for any discrete norm . defined on R N . From (28), we may write: Similarly, for the pseudoinverse Y † of Y, we can write: Thus, we have: We can see that Equations (29) and (31) confirm the bounds for the stability constant C. Calculating the pseudoinverse for approximating the system in Equation (27) numerically can be difficult computationally, but it ensures the stability. MATLAB's function condest can be used to estimate Y −1 ∞ in the case of square systems; thus, we have:

Numerical Results and Discussion
The proposed Laplace transform based local meshless method is tested on 2D linear multi-term FPIDEs. In our numerical experiments, we utilized the multiquadric (MQ) radial kernel defined by φ(r, ε)= 1 + (εr) 2 . For the optimal shape parameter, the uncertainty principle due to [53] (e.g., in RBF methods, we cannot have both good accuracy and good conditioning at the same time) is utilized. We performed our experiments in MATLAB R2019a on a Windows 10(64 bit) PC equipped with an Intel(R) Core(TM) i5-3317U CPU @ 1.70 GHz and with 4 GB of RAM. We chose L = T = 1, and d ρ = 1, for ρ = 1, 2, 3, . . . , m. Let U(ζ, τ) be the exact solution and U k (ζ, τ) be the numerical solution. To validate the theoretical results, we used the L ∞ error defined as: In our numerical experiments, we consider the problem (1) with initial condition Θ = 0, and the inhomogeneous term is given as: The problem is solved with two types of boundary conditions, the Dirichlet boundary conditions generated from the exact solution given by: U(x, y, τ) = τ 4 5 sin(πx) sin(πy), and the Robin boundary conditions given as: We solve the problem in the square domain, nut-shaped domain, and L-shaped domain.

Square Domain
In the first test, the problem is solved in square domain [0, 1] 2 with Dirichlet boundary conditions. The domain is discretized with regularly distributed nodes, as shown in the Figure  1. Then, the proposed scheme is applied to the 2D multi-term FPIDE. The exact and approximate solutions of the problem are presented in Figure 2a,b. The computational results obtained for various points N ∈ Ω, n ∈ Ω i , and quadrature points M along the contour Γ are given in Table 1. The L ∞ error, shape parameter ε, error estimate, condition number κ, and the computational time (CPU (s)) are shown in Table 1. The obtained results ensure the efficiency and stability of the method.

Nut-Shaped Domain
In the second test, the problem is solved in the nut-shaped domain with Dirichlet boundary conditions. The domain is discretized with regularly distributed nodes, as shown in the Figure 3a. Then, the proposed scheme is applied to the 2D multi-term FPIDE. The exact and the numerical solutions of the problem are presented in Figure 3b. The computational results obtained for various points N ∈ Ω, n ∈ Ω i , and quadrature points M along the contour Γ are shown in Table 2. The L ∞ error, shape parameter ε, error estimate, condition number κ, and the computational time (CPU (s)) are shown in Table 2. The obtained results ensure the efficiency of the proposed method for problems defined in irregular domains.

L-Shaped Domain
In the third test, the problem is solved in the L-shaped domain with Dirichlet and Robin boundary conditions. The domain is discretized with regularly distributed nodes, as shown in the Figure 4a. Then, the proposed method is applied to the 2D multi-term FPIDE. The exact and the numerical solutions of the problem are presented in Figure 4b. The computational results obtained for various points N ∈ Ω, n ∈ Ω i , and quadrature points M along the contour Γ with Dirichlet boundary conditions are shown in Table 3. The L ∞ error, shape parameter ε, error estimate, condition number κ, and the computational time (CPU (s)) are shown in Table 3. Figure 5 shows the absolute error obtained using Robin boundary conditions. The obtained results ensure the efficiency of the proposed method for problems defined in irregular domains.

Conclusions
We successfully coupled the Laplace transform and localized meshless method for the approximation of the solution of the multi-term 2D FPIDE. The time stepping procedure was avoided via the Laplace transform, and the issues arising due to dense differentiation matrices were resolved via the localized meshless method. For the contour integration, we utilized the recently introduced improved Talbot's contour. The convergence and stability of the method were discussed. To validate the numerical scheme and check its efficiency, the numerical experiments were carried out in the square, nut-shaped, and L-shaped domains. From the results obtained, it was observed that the proposed numerical scheme is efficient and has better accuracy compared to other available work. It was observed that the improved Talbot's method is computationally more useful than other available methods. The results led us to the conclusion that the proposed method is capable of solving FPIDEs without time instability in less computation time.