On the Numerical Approximation of Three-Dimensional Time Fractional Convection-Diffusion Equations

In this paper, we present an efficient method for the numerical investigation of three-dimensional non-integer-order convectiondiffusion equation (CDE) based on radial basis functions (RBFs) in localized form and Laplace transform (LT). In our numerical scheme, first we transform the given problem into Laplace space using Laplace transform. %en, the local radial basis function (LRBF) method is employed to approximate the solution of the transformed problem. Finally, we represent the solution as an integral along a smooth curve in the complex left half plane. %e integral is then evaluated to high accuracy by a quadrature rule. %e Laplace transform is used to avoid the classical time marching procedure. %e radial basis functions are important tools for scattered data interpolation and for solving partial differential equations (PDEs) of integer and non-integer order. %e LRBF and global radial basis function (GRBF) are used to produce sparse collocation matrices which resolve the issue of the sensitivity of shape parameter and ill conditioning of system matrices and reduce the computational cost. %e application of Laplace transformation often leads to the solution in complex plane which cannot be generally inverted. In this work, improved Talbot’s method is utilized which is an efficient method for the numerical inversion of Laplace transform.%e stability and convergence of the method are discussed. Two test problems are considered to validate the numerical scheme.%e numerical results highlight the efficiency and accuracy of the proposed method.


Introduction
Fractional calculus is a fast growing field and has attracted the research community due to its applications in engineering and other sciences [1][2][3]. Fractional-order operators are non-local, i.e., they include the history from the initial states to the current state. erefore, fractional operators are often applied to model systems to describe the influence of memory effects [4][5][6]. Fractional-order operators generalize the notion of integer-order operators [7]. Fractional-order operators can be used to model many processes such as the signal process and anomalous diffusion [8][9][10][11]. erefore, efficient and accurate numerical methods for solving physical models involving fractional-order derivatives become an important research area [12][13][14][15][16].
e CDE is a mixture of the equations of diffusion and convection and explains physical phenomena in which particles, electricity, or other physical quantities are transmitted within a physical structure through two procedures: convection and diffusion. Fractional-order CDEs are regarded as the extended form of ordinary CDEs. Fractional-order CDEs can express physical problems more accurately as compared to ordinary CDEs. e fractionalorder CDE has applications in many phenomena such as the mass transport [17], environmental science [18], gas transport through heterogeneous soil and gas reservoirs [19], contaminant or liquid transport in complex media and heat conduction [20], and global weather production, in which an initially discontinuous profile is propagated by diffusion and convection, the latter with a speed of c [21]. e solution of fractional-order CDE has been investigated by many authors, and they have developed many robust analytical and numerical methods in this regard. For example, the authors in [22] have studied the analytical solution of fractional-order diffusion and convection equations. In [23], the authors have utilized the Adomian decomposition method for solving the fractional-order CDE. Momani and Yıldırım in [24] have obtained the approximate analytical solution of CDE via He's homotopy perturbation method. Since the analytical solution of fractional-order problems in most of the cases is difficult to obtain, the researchers have developed numerous numerical methods for the investigation of the solutions of fractional-order problems such as the finite difference method [25], finite element method [26], finite volume method [27], and meshless methods [28][29][30]. e authors in [31] studied the solution CDEs of fractional via the discontinuous Galerkin method. A box-type scheme in [32] was developed for the sub-diffusion equations of fractional order. Hristov in [33] studied the solution of fractional sub-diffusion equation via an integral method. Zhai et al. [34] proposed a compact finite difference scheme for the numerical simulation of fractional-order 3D CDE. In [8], the authors developed a numerical scheme for anomalous sub-diffusion equation. ey have studied the stability and convergence of the scheme using the energy method. Chen et al. [9] proposed a numerical method with high accuracy for variable-order fractional anomalous sub-diffusion equation. ey discussed the convergence and stability using Fourier analysis. In [25], a compact difference scheme was developed for the fractional-order modified anomalous sub-diffusion and wave-diffusion equations. Zeng et al. [26] solved the time fractional sub-diffusion equation by utilizing the finite difference and finite element approach. ey also proposed two fully discrete schemes in [35] for the fractional sub-diffusion equation with secondorder accuracy. In [36], the numerical scheme for the approximation of 2D sub-diffusion equation of fractional order was proposed. In [37], the authors proposed a numerical method for non-integer-order sub-diffusion equation. Gao and Sun [38] obtained the numerical approximation of fractional sub-diffusion equation via finite difference scheme. Hussain et al. [27] analyzed the solution of CDE numerically using the finite element method. In computational modeling, the well-known classical mesh-based methods such as finite difference, finite elements, and finite volume have received much interest in recent years. However, these algorithms were designed with constraints related to issues such as node selection and multidimensional scaling.
Recently, the meshless method has attracted more and more attention because it is among one of the powerful tools for handling PDEs. Especially the RBF-based meshless methods are one of the most popular types among these methods. e benefit of these techniques over other classical techniques such as the finite elements, finite differences, and boundary elements methods is that in these methods, no mesh or element is required for discretization of the domain or the boundary which is a difficult task especially in 3D problems. Instead, meshless methods only require a few scattered points, the connection of which is not important. In meshless methods, the refinement of mesh can easily be done with minimal cost. Meshless methods are easily extendible to multidimensional problems. Meshless methods have been proved successful for solving PDEs on both regular and irregular node arrangements. Meshless methods have been successfully applied to various engineering and other science problems. For example, in [28], a meshfree method based on finite differences was applied to diffusion equation. Oruç [11] applied the meshfree method based on finite differences for the numerical solution of Zakharov-Rubenchik equations. Dehghan and Shafieeabyaneh [29] utilized the RBF method based on finite difference to simulate the regularized long-wave equation (RLWE) and extended Fisher-Kolmogorov equation (EFKE) in higher dimensions. Mohebbi et al. [14] investigated the solution of 2D modified anomalous sub-diffusion equation of arbitrary order via the RBF method. e authors in [39] derived the error estimates of the numerical solution for reaction-sub-diffusion of fractional order via the meshless method. Wei et al. [30] obtained the numerical solution of 2D fractional-order diffusion equation based on the implicit RBF method. In [40], the authors approximated the solution of 2D fractional-order sub-diffusion equation via Kansa's method. Similarly, approximation of diffusion equation of variable fractional order with different boundary conditions based on local RBF method was carried out in [41]. Shivanian in [42] proposed an efficient local radial basis function method for the simulation of 2D fractional-order convection-diffusion reactions. e stability and convergence of the method was also proved theoretically. In [43], the numerical solutions of advection diffusion reaction equations on complex geometries using a local radial basis function method were studied. Golbabai and Kalarestaghi [44] proposed a local radial basis function method for dominated convection-diffusion equations.
However, in these time stepping schemes, the computations may be very expansive because each new iteration is dependent on the previous time step. An alternative way is to use the LT coupled with these numerical methods. In this article, we propose the LRBF method coupled with Laplace transform. e LT is used to avoid the stability restrictions, which are commonly encountered in time stepping procedures. e LRBF method is used to resolve the issue of ill conditioning of the differentiation matrices and the sensitivity of shape parameter in the GRBF method. e main idea of the LRBF method is the collocation on overlapping sub-domains of the whole domain. e overlapping subdomains remarkably reduce the size of collocation matrix by solving many small size matrices. Each small matrix has the same size as the number of nodes in the domain of influence of each node. e Caputo fractional derivatives are considered here because they allow traditional initial and boundary conditions to be included in the formulation of the problem. In this paper, we consider the linear fractionalorder CDE, and the fractional derivatives are taken in Caputo sense as follows: and where Ω is the domain and zΩ is its boundary. w(ξ) is a function of ξ and V(ξ) is a vector. e forcing term h(ξ, t) is sufficiently smooth and g(ξ, t) and U 0 are given continuous functions. L b is the boundary differential operator, and c 0 D α t U(ξ, t) is the Caputo fractional derivative of order α of the function U(ξ, t) which is defined as (4) e paper is organized as follows. In Section 2, the Laplace transform method for the given model is discussed. Section 3 describes the local RBF approximation of the transformed problem in Laplace space. Section 4 is about the stability of the discrete system. Section 5 is devoted to the approximation of inverse Laplace transform. In Section 6, the numerical scheme is applied to 3D CDE. Finally, in Section 7, the conclusion is drawn.

Time Discretization via Laplace Transform
Here, we apply LT to given model defined in (1)-(3) to avoid the time marching method. e LT of piecewise continuous function Also, the LT of the Caputo derivative c 0 D α t is given as Applying the LT to equations (1)-(3), we get which can be rearranged as where I denotes the identity operator, L � Δ + V(ξ) · ∇ − w(ξ)I is the governing differential operator, and L b is the boundary differential operator. In order to obtain the approximate solution of the system defined in (8) and (9), first we need to discretize the operators L and L b via LRBFs. Afterwards, systems (8) and (9) will be solved in parallel for each point s (see, e.g., [45]). In the final step, the solution of problems (1)-(3) can be obtained by expressing it as a Bromwich integral. e next section is about the approximation of the spatial operators via the LRBF method.

Space Discretization via LRBF Method
Here, we propose a LRBF for approximating the three-dimensional elliptic PDEs in Laplace space. Unlike the GRBF method, the LRBF method can be employed only for small number of nodes (sub-domain) at each nodal point instead of implementing the entire set of points. is method provides a linear system that is sparse and well conditioned.
denotes that the nodes belonging to the local domain Ω i of each node ξ i . To construct the LRBFs, the collocation technique is applied to the local domain via the LRBF method is given as and by applying the collocation method in Ω i , we have . .
Mathematical Problems in Engineering 3 Denote the matrix in equation (11) by j ‖) is a radial kernel, ‖.‖ is the Euclidean norm, and λ [i] j N j�1 are the unknown coefficients.
ere are a large number of RBFs which are used commonly including the Gaussian (GA), multiquadrics (MQ), thin-plate splines (TPS), and inverse multiquadrics (IMQ). Among these RBFs, the multiquadric (MQ) is perhaps the most popular that is used in applications [41]; therefore, in this study, we have selected the MQ, which is given as where S p denotes the shape parameter. Each node ξ i and its n nearest neighbor points are called stencils [43]. A typical stencil for n � 7 is shown in Figure 1.
At each stencil, we obtain an n × n matrix of the form where U In [46], the invertibility of the system matrix G [i] based on strictly positive definite kernel functions for distinct nodes was proved. e coefficients λ [i] in (13) can be obtained as e next step in the LRBF scheme is as follows: and from the above equation, we can write Hence, we have where (17) gives us the local form, and the global form can be obtained by expanding the matrix Ψ [i] to Ψ by adding zeros to corresponding position in each row [41]. us, the differential operator L is approximated at each center ξ i via the LRBF method as e matrix Ψ N×N is sparse with n non-zero entries and N − n zero entries. We can similarly approximate the operator L b as and substituting equations (18) and (19) in equations (8) and (9), we get Solving the above system defined in equations (20) and (21) for each point s, we will obtain the approximate solution U ∧ (ξ, s).

Algorithm for Optimal Shape Parameter.
e MQ kernel function defined in (12) includes a parameter S p which can be varied on every local domain to enhance the accuracy of numerical solution. To estimate the accuracy of approximate solution and quantify the sensitivity to perturbation of the linear system, the condition number (C N ) of the interpolation matrix G [i] may be used. We utilized the uncertainty principle [47] for a decent estimation of the shape parameter S p . In RBF methods for ill-conditioned matrices, better accuracy is achieved. Using the uncertainty principle, when (C N ) is kept approximately in the range 10 12 < C N < 10 16 , smallest error occurs. e matrix G [i] can be expressed as Q, S, V � sv d(G [i] ), where S n×n is diagonal matrix which contains the singular values of G [i] . Q n×n , V n×n are orthogonal matrices, and the condition number C N can be computed as C N � ‖G [i] ‖‖(G [i] ) − 1 ‖ � max(S)/min(S). e following algorithm in MATLAB can be used to search for optimal shape parameter (Algorithm 1).
After getting an optimal shape parameter S p , the svd is used to compute ( [48]). Hence, the weights Ψ [i] in (18) can be computed.

Stability
For the stability of systems (20) and (21), we write the system as A N×N is sparsely obtained via the LRBF method. We define a constant for the aforementioned system as called the constant of stability. Also, using norm ‖.‖ on R N Q is finite. From (23), we have Similarly, we have where A † is the pseudoinverse of A, and thus we have Equations (24) and (26) confirm the bounds for Q. e calculation of the pseudoinverse for the approximation of the system defined in (22) may be hard computationally, but it validates the stability. e MATLAB function condest provides the estimation for ‖A − 1 ‖ ∞ , and thus

Numerical Inversion of Laplace Transform
After solving the system defined in (20) and (21) for each node s, we obtain the solution of (1)-(3) via the inversion of Laplace transformation as In applying the LT method, the main hurdle is the inversion of the LT. In most of the cases, it is quite difficult to invert the LT analytically. erefore, we need some numerical techniques for the inversion of LT. In the literature, a large number of techniques are available for the inversion of LT [49], among which one of the efficient and easy methods in implementation is the method of contour integration used in combination with trapezoidal and midpoint rules. To handle the exponential factor in equation (28), contour deformation is performed. Particularly, the contour in equation (28) can be deformed to a Hankel contour [50]. e exponential factors on these contours decay rapidly making the integral in equation (28) appropriate for solving via trapezoidal or midpoint rule [50][51][52]. In this work, we consider a contour having the following parametric form: where Res( ± π) � −∞, and where ϑ(β) � −9 + ξβ cot(ρβ) + ζιβ, (31) and the parameters 9, ρ, ξ, and ζ must be selected by the user. From equations (28) and (30), we get e M-panel midpoint rule with k � (2π/M) is used to evaluate the integral in equation (32) numerically as for s j � s β j , s j ′ � s ′ β j and β j � −π + j − 1 2 k.

Convergence and Accuracy.
While obtaining the numerical solution of problems (1)-(3), we utilized the LT in combination with the LRBF method. First we apply the LT to reduce the problem to an equivalent elliptic PDE, and no error occurs in this step. In the second step, we employ the LRBF for approximating the solution of the elliptic PDE. e error estimate of the LRBF method is of order O(μ (1/S p h) ), 0 < μ < 1, where h and S p are the fill distance and the shape parameter, respectively [53]. While solving the integral in equation (32) numerically, the convergence of the quadrature rule is achieved at different time rates which rely on the step size k, the domain [t 0 , T], and the path Λ. e error analysis of improved Talbot's method is based on the following theorem.
Step 1: set C N � 1 Step 2: select 10 12 < C N < 10 16 Step 3: while C N > C N maximum and C N < C N minimum

Step 4: Q, S, V � sv d(G [i] )
Step 5: C N � (max(S)/min(S)) Step 6: if C N < C N minimum , S p � S p − S p Increament Step 7: else C N > C N maximum , S p � S p + S p Increament return S p . ALGORITHM 1: Algorithm for optimal shape parameter.

Mathematical Problems in Engineering
Theorem 1 (see [50]). Let β j be defined as in (33). If f: D ⟶ C is an analytic function in the set D � β ∈ C: −π < Reβ < π and − d < Imβ < c , (34) when c, d > 0, then where ∀0 < ϑ < c and 0 < θ < d and M even; for odd M, we can replace tan(Mβ/2) with −cot (Mβ/2), and by analyzing the behavior of tangent function in complex, this can be bounded as where M is assumed to be even and C and c are some positive constants. For odd M, the analysis is similar.
For optimal accuracy, one needs optimal parameters in equation (30). e authors [50] have obtained the optimal parameters given as

Numerical Experiments
To check the performance of the proposed transformed LRBF method, we have used three error measures. ey are the absolute error Error 1, the maximum absolute error Error 2, and the relative error Error 3 which are defined as follows: Error 2 � max error abs , where U(ξ, τ) and U k (ξ, τ) are the analytical and numerical solutions, respectively. We test our proposed method using two examples. (1) with exact solution as [54] U

Problem 1. Consider equation
where and h(ξ, t) � 0. e accuracy of the problem is investigated at t � 1 in cube domain [−1, 1] 3 with α � 0.1, 0.5, 0.9. Table 1 shows the maximum absolute error Error 2 and relative error Error 3 for several nodes N, stencil size n � 80, and M � 28. Table 2 shows the maximum absolute error Error 2 and relative error Error 3 for several stencil sizes n, N � 1728, and M � 26. It is clear that the method has approximated the problem with acceptable accuracy. We observe that the proposed method produced accurate results. e numerical solution in rectangular area with different values of α is given in Figures 2(a) Figures 6(a)-6(c). Also, the graph of absolute error for α � 0.9 with N � 225, n � 81, M � 26, and y � −0.5 is shown in Figure 7(a) and that for α � 0.95 with N � 400, n � 25, M � 26, and y � −0.5 is shown Figure 7(b). A high accuracy is observed. (1) with exact solution given as [55] 6

Problem 2. Consider equation
Mathematical Problems in Engineering U(x, y, z, t) � t (α+3) sin(πx)sin(πy)sin(πz), t > 0, with w(ξ) � 1, V(ξ) � (1, 1, 1) T and the forcing term h(ξ, t) is specified to satisfy the given equation and the exact solution. e initial and boundary data are extracted from the exact solution. e accuracy of the problem is investigated in cube domain [0, 1] 3 at t � 1 for α � 0.1, 0.5, 0.9. Table 3 shows the maximum absolute error Error 2 and relative error Error 3 for several nodes N, stencil size n � 80, and M � 28. Table 4 shows the maximum absolute error Error 2 and relative error Error 3 for several values of M, N � 4913, and stencil size          Figures 13(a)-13(c). Also, the graph of absolute error for α � 0.9 with N � 256, n � 25, M � 28, and z � 0.2 is shown in Figure 14(a) and that for α � 0.95 with N � 625, n � 38, M � 28, and z � 0.8 is shown Figure 14(b). From all the obtained results presented in figures and tables, we conclude that the proposed method is stable, accurate, and efficient.

Conclusion
In this work, we constructed a method combining LT with the LRBF method for the solution of time fractional CDE in terms of Caputo derivative. e main advantage of this method is obtaining the solution at final time and avoiding the time stepping procedure. e LRBF produced sparse differentiation matrices due to which the ill-conditioning issues of the collocation matrices and the shape parameter's sensitivity were resolved. e convergence of the method was discussed, and the bounds for the stability constant of the fully discrete system were derived. e accuracy of the method was tested by means of two examples, and the results were compared with those from other schemes. e obtained results led us to the conclusion that the proposed method is powerful and effective to find the numerical solutions of three-dimensional time-dependent fractional PDEs, so it can be also applied to a wide range of complex problems that occur in natural sciences and engineering.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
All authors contributed equally to this study.