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

Authors: Kamran; Raheel Kamal; Gul Rahmat; Kamal Shah
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, we present an efficient method for the numerical investigation of three-dimensional non-integer-order convection-diffusion 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. Then, 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. The integral is then evaluated to high accuracy by a quadrature rule. The Laplace transform is used to avoid the classical time marching procedure. The radial basis functions are important tools for scattered data interpolation and for solving partial differential equations (PDEs) of integer and non-integer order. The 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. The 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. The stability and convergence of the method are discussed. Two test problems are considered to validate the numerical scheme. The numerical results highlight the efficiency and accuracy of the proposed method.


---

## Body

## 1. Introduction

Fractional calculus is a fast growing field and has attracted the research community due to its applications in engineering and other sciences [1–3]. Fractional-order operators are non-local, i.e., they include the history from the initial states to the current state. Therefore, fractional operators are often applied to model systems to describe the influence of memory effects [4–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–11]. Therefore, efficient and accurate numerical methods for solving physical models involving fractional-order derivatives become an important research area [12–16].The 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. The fractional-order 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]. The 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–30].The 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. They 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. They 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. They also proposed two fully discrete schemes in [35] for the fractional sub-diffusion equation with second-order 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. The 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. The 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. The 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. The LT is used to avoid the stability restrictions, which are commonly encountered in time stepping procedures. The 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. The main idea of the LRBF method is the collocation on overlapping sub-domains of the whole domain. The overlapping sub-domains 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. The 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 fractional-order CDE, and the fractional derivatives are taken in Caputo sense as follows:(1)D0ctαUξ¯,t=ΔUξ¯,t+Vξ¯⋅∇Uξ¯,t−wξ¯Uξ¯,t+hξ¯,t,ξ¯∈Ω,t>0,satisfying(2)Uξ¯,0=U0,ξ¯∈Ω,and(3)ℒbUξ¯,t=gξ¯,t,ξ¯∈∂Ω,where Ω is the domain and ∂Ω is its boundary. wξ¯ is a function of ξ¯ and Vξ¯ is a vector. The forcing term hξ¯,t is sufficiently smooth and gξ¯,t and U0 are given continuous functions. ℒb is the boundary differential operator, and D0ctαUξ¯,t is the Caputo fractional derivative of order α of the function Uξ¯,t which is defined as(4)D0ctαUξ¯,t=1Γγ−α∫0t∂γUξ¯,s/∂sγt−sα+1−γds,γ−1<α<γ,γ∈ℕ,t>0.The paper is organized as follows. In Section2, 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.


## 2. Time Discretization via Laplace Transform

Here, we apply LT to given model defined in (1)–(3) to avoid the time marching method. The LT of piecewise continuous function Ut defined on 0,∞ is defined by(5)U∧s=LUt=∫0∞e−stUtdt.Also, the LT of the Caputo derivativeD0ctα is given as(6)LD0ctαUt=sαU∧s−∑i=0m−1sα−i−1Ui0.Applying the LT to equations (1)–(3), we get(7)sαU∧ξ¯,s−sα−1U0−ℒU∧ξ¯,s=h∧ξ¯,s,ℒbU∧ξ¯,s=g∧ξ¯,s,which can be rearranged as(8)sαI−ℒU∧ξ¯,s=U0+h∧ξ¯,s,ξ¯∈Ω,(9)ℒbU∧ξ¯,s=g∧ξ¯,s,ξ¯∈∂Ω,where I denotes the identity operator, ℒ=Δ+Vξ¯⋅∇−wξ¯I is the governing differential operator, and ℒ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 ℒ and ℒ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. The next section is about the approximation of the spatial operators via the LRBF method.


## 3. 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. This method provides a linear system that is sparse and well conditioned. LetΘ=ξ¯ii=1N be the set of centers in Ω¯=Ω∪∂Ω. For each node ξ¯i, we select n nearest neighbor centers ξ¯ji,j=1,2,…,n. The index i 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 Ωi=ξ¯jij=1n,i=1,2,…,N. The approximation of U∧ξ¯i,s via the LRBF method is given as(10)U∧ξ¯i,s=∑j=1nλjiϕξ¯i−ξ¯ji,and by applying the collocation method in Ωi, we have(11)U∧ξ¯1i,sU∧ξ¯2i,s...U∧ξ¯ni,s=ϕξ¯1i−ξ¯1iϕξ¯1i−ξ¯2i...ϕξ¯1i−ξ¯niϕξ¯2i−ξ¯1iϕξ¯2i−ξ¯2i...ϕξ¯2i−ξ¯ni..................ϕξ¯ni−ξ¯1iϕξ¯ni−ξ¯2i...ϕξ¯ni−ξ¯niλ1iλ2i...λni.Denote the matrix in equation (11) by Gi, where ϕξ¯i−ξ¯ji is a radial kernel, . is the Euclidean norm, and λjij=1N are the unknown coefficients. There 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(12)ϕξ¯i−ξ¯ji=1+Spξ¯i−ξ¯ji2,where Sp 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.Figure 1
A typical stencil forn=7.At each stencil, we obtain ann×n matrix of the form(13)U∧i=Giλi,i=1,2,…,N,where U∧i=U∧ξ¯1i,U∧ξ¯2i,…,U∧ξ¯niT and λi=λ1i,λ2i,…,λniT. In [46], the invertibility of the system matrix Gi based on strictly positive definite kernel functions for distinct nodes was proved. The coefficients λi in (13) can be obtained as(14)λi=Gi−1U∧i,fori=1,2,…,N.The next step in the LRBF scheme is as follows:(15)ℒU∧ξ¯i=∑j=1nλjiℒϕξ¯i−ξ¯ji,and from the above equation, we can write(16)ℒU∧ξ¯i=λi⋅ϒi=ϒi⋅λi=ϒiWi−1U∧i.Hence, we have(17)ℒU∧ξ¯i=Ψiu∧i,where Ψi=ϒiWi−1 and ϒi=ℒϕξ¯i−ξ¯ji is a 1×n vector. Equation (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]. Thus, the differential operator ℒ is approximated at each center ξ¯i via the LRBF method as(18)ℒU∧≡Ψu∧.The matrixΨN×N is sparse with n non-zero entries and N−n zero entries. We can similarly approximate the operator ℒb as(19)ℒbU∧≡Ξu∧,and substituting equations (18) and (19) in equations (8) and (9), we get(20)sαI−ΨU∧ξ¯,s=u0+h∧ξ¯,s,ξ¯∈Ω,(21)ΞU∧ξ¯,s=g∧ξ¯,s,ξ¯∈∂Ω.Solving the above system defined in equations (20) and (21) for each point s, we will obtain the approximate solution U∧ξ¯,s.


### 3.1. Algorithm for Optimal Shape Parameter

The MQ kernel function defined in (12) includes a parameter Sp 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 (CN) of the interpolation matrix Gi may be used. We utilized the uncertainty principle [47] for a decent estimation of the shape parameter Sp. In RBF methods for ill-conditioned matrices, better accuracy is achieved. Using the uncertainty principle, when (CN) is kept approximately in the range 1012<CN<1016, smallest error occurs. The matrix Gi can be expressed as Q,S,V=svdGi, where Sn×n is diagonal matrix which contains the singular values of Gi. Qn×n, Vn×n are orthogonal matrices, and the condition number CN can be computed as CN=GiGi−1=maxS/minS.The following algorithm in MATLAB can be used to search for optimal shape parameter (Algorithm1).Algorithm 1: Algorithm for optimal shape parameter.
Step 1: setCN=1Step 2: select1012<CN<1016Step 3:whileCN>CNmaximumandCN<CNminimumStep 4:Q,S,V=svdGiStep 5:CN=maxS/minSStep 6:ifCN<CNminimum,Sp=Sp−SpIncreamentStep 7:elseCN>CNmaximum,Sp=Sp+SpIncreamentreturnSp.After getting an optimal shape parameterSp, the svd is used to compute Gi−1=QSVT−1=VS−1QT (see [48]). Hence, the weights Ψi in (18) can be computed.


## 3.1. Algorithm for Optimal Shape Parameter

The MQ kernel function defined in (12) includes a parameter Sp 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 (CN) of the interpolation matrix Gi may be used. We utilized the uncertainty principle [47] for a decent estimation of the shape parameter Sp. In RBF methods for ill-conditioned matrices, better accuracy is achieved. Using the uncertainty principle, when (CN) is kept approximately in the range 1012<CN<1016, smallest error occurs. The matrix Gi can be expressed as Q,S,V=svdGi, where Sn×n is diagonal matrix which contains the singular values of Gi. Qn×n, Vn×n are orthogonal matrices, and the condition number CN can be computed as CN=GiGi−1=maxS/minS.The following algorithm in MATLAB can be used to search for optimal shape parameter (Algorithm1).Algorithm 1: Algorithm for optimal shape parameter.
Step 1: setCN=1Step 2: select1012<CN<1016Step 3:whileCN>CNmaximumandCN<CNminimumStep 4:Q,S,V=svdGiStep 5:CN=maxS/minSStep 6:ifCN<CNminimum,Sp=Sp−SpIncreamentStep 7:elseCN>CNmaximum,Sp=Sp+SpIncreamentreturnSp.After getting an optimal shape parameterSp, the svd is used to compute Gi−1=QSVT−1=VS−1QT (see [48]). Hence, the weights Ψi in (18) can be computed.


## 4. Stability

For the stability of systems (20) and (21), we write the system as(22)AU∧ξ¯,s=b.AN×N is sparsely obtained via the LRBF method. We define a constant for the aforementioned system as(23)Q=supU∧≠0U∧ξ¯,sAU∧ξ¯,s,called the constant of stability. Also, using norm . on RNQ is finite. From (23), we have(24)A−1≤U∧ξ¯,sAU∧ξ¯,s≤Q.Similarly, we have(25)A†=supv≠0A†vv,where A† is the pseudoinverse of A, and thus we have(26)A†≥supv=AU∧≠0A†AU∧ξ¯,sAU∧ξ¯,s=supU∧≠0U∧ξ¯,sAU∧ξ¯,s=Q.Equations (24) and (26) confirm the bounds for Q. The calculation of the pseudoinverse for the approximation of the system defined in (22) may be hard computationally, but it validates the stability. The MATLAB function condest provides the estimation for A−1∞, and thus(27)Q=condestA′A∞.


## 5. 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(28)Uξ¯,τ=12πi∫σ−i∞σ+i∞esτU∧ξ¯,sds.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. Therefore, 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]. The exponential factors on these contours decay rapidly making the integral in equation (28) appropriate for solving via trapezoidal or midpoint rule [50–52]. In this work, we consider a contour having the following parametric form:(29)Λ:s=sβ,−π≤β≤π,where Res±π=−∞, and(30)sβ=Mϑβt,where(31)ϑβ=−ϱ+ξβcotρβ+ζιβ,and the parameters ϱ,ρ,ξ, and ζ must be selected by the user. From equations (28) and (30), we get(32)Uξ¯,τ=12πi∫ΛesτU∧ξ¯,sds=12πi∫−ππesβτU∧ξ¯,sβs′βdβ.TheM-panel midpoint rule with k=2π/M is used to evaluate the integral in equation (32) numerically as(33)Ukξ¯,τ=1Mi∑j=1MesjτU∧ξ¯,sjs′j,forsj=sβj,sj′=s′βjandβj=−π+j−12k.


### 5.1. 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. The error estimate of the LRBF method is of order Oμ1/Sph, 0<μ<1, where h and Sp 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 t0,T, and the path Λ. The error analysis of improved Talbot’s method is based on the following theorem.Theorem 1 (see [50]).
Letβj be defined as in (33). If f:D⟶ℂ is an analytic function in the set(34)D=β∈ℂ:−π<Reβ<πand−d<Imβ<c,when c,d>0, then(35)∫−ππfβdβ−2πM∑j=1Mfβj=B−θ+B+ϑ,where(36)B+ϑ=12∫−π−π+ιϑ+∫−π+ιϑπ+ιϑ+∫π+ιϑπ1+ιtanMβ2fβdβ,B−θ=12∫−π−π−ιϑ+∫−π−ιϑπ−ιϑ+∫π−ιϑπ1−ιtanMβ2fβdβ,∀0<ϑ<cand0<θ<dandMeven; for odd M, we can replace tanMβ/2 with −cotMβ/2, and if fβ is real valued, that is, fβ¯=fβ¯, then(37)Bϑ=B+ϑ+B−ϑ=Re∫−π+ιϑπ+ιϑ1+ιtanMβ2fβdβ,and by analyzing the behavior of tangent function in complex, this can be bounded as(38)Bϑ≤4πCexpcM−1,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). The authors [50] have obtained the optimal parameters given aswith the error estimate(39)errorest=Uξ¯,t−Ukξ¯,t=Oe−1.3580M.


## 5.1. 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. The error estimate of the LRBF method is of order Oμ1/Sph, 0<μ<1, where h and Sp 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 t0,T, and the path Λ. The error analysis of improved Talbot’s method is based on the following theorem.Theorem 1 (see [50]).
Letβj be defined as in (33). If f:D⟶ℂ is an analytic function in the set(34)D=β∈ℂ:−π<Reβ<πand−d<Imβ<c,when c,d>0, then(35)∫−ππfβdβ−2πM∑j=1Mfβj=B−θ+B+ϑ,where(36)B+ϑ=12∫−π−π+ιϑ+∫−π+ιϑπ+ιϑ+∫π+ιϑπ1+ιtanMβ2fβdβ,B−θ=12∫−π−π−ιϑ+∫−π−ιϑπ−ιϑ+∫π−ιϑπ1−ιtanMβ2fβdβ,∀0<ϑ<cand0<θ<dandMeven; for odd M, we can replace tanMβ/2 with −cotMβ/2, and if fβ is real valued, that is, fβ¯=fβ¯, then(37)Bϑ=B+ϑ+B−ϑ=Re∫−π+ιϑπ+ιϑ1+ιtanMβ2fβdβ,and by analyzing the behavior of tangent function in complex, this can be bounded as(38)Bϑ≤4πCexpcM−1,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). The authors [50] have obtained the optimal parameters given aswith the error estimate(39)errorest=Uξ¯,t−Ukξ¯,t=Oe−1.3580M.


## 6. Numerical Experiments

To check the performance of the proposed transformed LRBF method, we have used three error measures. They are the absolute errorError1, the maximum absolute error Error2, and the relative error Error3 which are defined as follows:(40)Error1=Uξ¯i−Ukξ¯i,Error2=maxerrorabs,Error3=∑i=1NUξ¯i−Ukξ¯i2∑i=1NUξ¯i2,where Uξ¯,τ and Ukξ¯,τ are the analytical and numerical solutions, respectively. We test our proposed method using two examples.


### 6.1. Problem 1

Consider equation (1) with exact solution as [54](41)Ux,y,z,t=eγ1t−γ2x+γ2y+γ2z,Ω=−1,13,t>0,where γ2=γ1/6,γ1=0.1,Vξ=−γ2,−γ2,−γ2T, wξ=0, and hξ¯,t=0. The accuracy of the problem is investigated at t=1 in cube domain −1,13 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. The numerical solution in rectangular area with different values of α is given in Figures 2(a)–2(c). In Figures 3(a) and 3(b), the exact and numerical solutions in the xz-plane for y=−1,−0.5,0,0.5,1 are shown; it is observed that the exact and numerical solutions are in better agreement. Figures 4(a)–4(c) show the plots of absolute errors for different values of the fractional order α with N=3375, n=81, and M=26. From the figures, we see that the method has produced accurate results. In Figures 5(a)–5(c), the contour slice plots of the absolute errors for various values of α in rectangular area using N=3375, n=81, and M=26 are presented. We also plot the contour graphs of absolute error for various values of α with N=3375,n=81,M=26, and z=0.5 in 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.Table 1
The computational results for problem 1 withN in the domain −1,13.
n=80,M=28N⟶12521672917283375α=0.1Error28.00 × 10−33.89 × 10−44.38 × 10−44.49 × 10−44.88 × 10−4Error31.60 × 10−31.21 × 10−41.30 × 10−41.40 × 10−41.55 × 10−4α=0.5Error28.00 × 10−31.00 × 10−31.10 × 10−31.10 × 10−31.00 × 10−3Error31.80 × 10−33.47 × 10−44.08 × 10−44.32 × 10−44.36 × 10−4α=0.9Error28.00 × 10−35.78 × 10−46.02 × 10−45.97 × 10−45.61 × 10−4Error31.80 × 10−31.92 × 10−42.23 × 10−42.35 × 10−42.31 × 10−4Table 2
The computational results for problem 1 withN in the domain −1,13.
M=26,N=1728n⟶4050608385α=0.1Error21.20 × 10−31.50 × 10−34.60 × 10−34.48 × 10−47.15 × 10−4Error32.00 × 10−41.81 × 10−44.55 × 10−41.40 × 10−45.87 × 10−4α=0.5Error21.60 × 10−31.50 × 10−34.60 × 10−31.10 × 10−37.00 × 10−4Error34.58 × 10−44.58 × 10−45.75 × 10−44.32 × 10−45.74 × 10−4α=0.9Error21.40 × 10−31.40 × 10−34.60 × 10−35.97 × 10−45.97 × 10−4Error32.78 × 10−42.73 × 10−44.64 × 10−42.35 × 10−45.50 × 10−4Figure 2
Approximate solution in rectangular area using various values ofα with N=1728, n=40,andM=26.
(a)(b)(c)Figure 3
Exact and numerical solutions of problem 1 withα=0.1 and y=−1,−0.5,0,0.5,1.
(a)(b)Figure 4
Absolute error in rectangular area using different values ofα with N=3375, n=81,andM=26.
(a)(b)(c)Figure 5
Contour slice of absolute error in rectangular area using different values ofα with N=3375, n=81,andM=26.
(a)(b)(c)Figure 6
Contour plots of absolute error in rectangular area using different values ofα with N=225, n=81,M=26, and y=0.5.
(a)(b)(c)Figure 7
(a) Plot of absolute error withα=0.9, N=225, n=81,M=26, and y=−0.5. (b) Plot of absolute error with α=0.95, N=400, n=25,M=26, and y=−0.5.
(a)(b)


### 6.2. Problem 2

Consider equation (1) with exact solution given as [55](42)Ux,y,z,t=tα+3sinπxsinπysinπz,t>0,with wξ=1,Vξ=1,1,1T and the forcing term hξ¯,t is specified to satisfy the given equation and the exact solution. The initial and boundary data are extracted from the exact solution. The accuracy of the problem is investigated in cube domain 0,13 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 n=83. The numerical solutions in rectangular area for different values of α are presented in Figures 8(a)–8(c). The slice plots of exact and numerical solutions using α=0.85,N=2744,n=85,andM=26 are shown in Figures 9(a) and 9(b). Exact and numerical solutions with α=0.1,N=4096,n=85,M=28, and z=0.4,0.6,0.8 are shown in Figures 10(a) and 10(b), respectively. Absolute error in rectangular area using various values of α with N=1000,n=75,andM=28 is presented in Figures 11(a)–11(c). The contour slice plots of absolute error in rectangular area using various values of α with N=4096,n=85,andM=28 are shown in Figures 12(a)–12(c). The contour graphs of absolute error using various values of α with N=4096,n=85,M=28, and z=0.5 are presented shown in 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.Table 3
The computational results for problem 1 withN in the domain 0,13.
n=80,M=28N⟶12521672917283375α=0.1Error21.90 × 10−38.07 × 10−46.66 × 10−47.23 × 10−48.39 × 10−4Error33.20 × 10−37.92 × 10−47.89 × 10−45.38 × 10−46.61 × 10−4α=0.5Error21.90 × 10−37.98 × 10−46.60 × 10−47.11 × 10−48.22 × 10−4Error33.10 × 10−37.85 × 10−47.81 × 10−45.31 × 10−46.47 × 10−4α=0.9Error21.80 × 10−37.82 × 10−46.48 × 10−46.90 × 10−47.92 × 10−4Error33.00 × 10−37.71 × 10−47.65 × 10−45.19 × 10−46.22 × 10−4[55]————2.35 × 10−3—Table 4
The computational results for problem 1 withN in the domain 0,13.
n=83,N=4913M⟶2022242628α=0.1Error27.14 × 10−47.14 × 10−47.15 × 10−47.15 × 10−47.15 × 10−4Error35.86 × 10−45.87 × 10−45.87 × 10−45.87 × 10−45.87 × 10−4α=0.5Error26.98 × 10−46.99 × 10−47.00 × 10−47.00 × 10−47.00 × 10−4Error35.72 × 10−45.74 × 10−45.74 × 10−45.74 × 10−45.74 × 10−4α=0.9Error26.66 × 10−46.72 × 10−46.72 × 10−46.72 × 10−46.72 × 10−4Error35.44 × 10−45.49 × 10−45.50 × 10−45.50 × 10−45.50 × 10−4[55]————2.35 × 10−3—Figure 8
Approximate solution in rectangular area using different values ofα with N=1000, n=75,andM=28.
(a)(b)(c)Figure 9
Slice plots of exact and approximate solutions of problem 2 withα=0.85,N=2744,n=85,andM=26.
(a)(b)Figure 10
Exact and approximate solutions of problem 2 withα=0.1 and z=0.4,0.6,0.8.
(a)(b)Figure 11
Absolute error in rectangular area using different values ofα with N=1000, n=75,andM=28.
(a)(b)(c)Figure 12
Contour slice of absolute error in rectangular area using different values ofα with N=4096, n=85,andM=28.
(a)(b)(c)Figure 13
Contour plots of absolute error in rectangular area using different values ofα with N=256, n=85,M=28, and z=0.5.
(a)(b)(c)Figure 14
(a) Plot of absolute error withα=0.9, N=256, n=25,M=28, and z=0.2. (b) Plot of absolute error with α=0.95, N=625, n=38,M=28, and z=0.8.
(a)(b)


## 6.1. Problem 1

Consider equation (1) with exact solution as [54](41)Ux,y,z,t=eγ1t−γ2x+γ2y+γ2z,Ω=−1,13,t>0,where γ2=γ1/6,γ1=0.1,Vξ=−γ2,−γ2,−γ2T, wξ=0, and hξ¯,t=0. The accuracy of the problem is investigated at t=1 in cube domain −1,13 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. The numerical solution in rectangular area with different values of α is given in Figures 2(a)–2(c). In Figures 3(a) and 3(b), the exact and numerical solutions in the xz-plane for y=−1,−0.5,0,0.5,1 are shown; it is observed that the exact and numerical solutions are in better agreement. Figures 4(a)–4(c) show the plots of absolute errors for different values of the fractional order α with N=3375, n=81, and M=26. From the figures, we see that the method has produced accurate results. In Figures 5(a)–5(c), the contour slice plots of the absolute errors for various values of α in rectangular area using N=3375, n=81, and M=26 are presented. We also plot the contour graphs of absolute error for various values of α with N=3375,n=81,M=26, and z=0.5 in 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.Table 1
The computational results for problem 1 withN in the domain −1,13.
n=80,M=28N⟶12521672917283375α=0.1Error28.00 × 10−33.89 × 10−44.38 × 10−44.49 × 10−44.88 × 10−4Error31.60 × 10−31.21 × 10−41.30 × 10−41.40 × 10−41.55 × 10−4α=0.5Error28.00 × 10−31.00 × 10−31.10 × 10−31.10 × 10−31.00 × 10−3Error31.80 × 10−33.47 × 10−44.08 × 10−44.32 × 10−44.36 × 10−4α=0.9Error28.00 × 10−35.78 × 10−46.02 × 10−45.97 × 10−45.61 × 10−4Error31.80 × 10−31.92 × 10−42.23 × 10−42.35 × 10−42.31 × 10−4Table 2
The computational results for problem 1 withN in the domain −1,13.
M=26,N=1728n⟶4050608385α=0.1Error21.20 × 10−31.50 × 10−34.60 × 10−34.48 × 10−47.15 × 10−4Error32.00 × 10−41.81 × 10−44.55 × 10−41.40 × 10−45.87 × 10−4α=0.5Error21.60 × 10−31.50 × 10−34.60 × 10−31.10 × 10−37.00 × 10−4Error34.58 × 10−44.58 × 10−45.75 × 10−44.32 × 10−45.74 × 10−4α=0.9Error21.40 × 10−31.40 × 10−34.60 × 10−35.97 × 10−45.97 × 10−4Error32.78 × 10−42.73 × 10−44.64 × 10−42.35 × 10−45.50 × 10−4Figure 2
Approximate solution in rectangular area using various values ofα with N=1728, n=40,andM=26.
(a)(b)(c)Figure 3
Exact and numerical solutions of problem 1 withα=0.1 and y=−1,−0.5,0,0.5,1.
(a)(b)Figure 4
Absolute error in rectangular area using different values ofα with N=3375, n=81,andM=26.
(a)(b)(c)Figure 5
Contour slice of absolute error in rectangular area using different values ofα with N=3375, n=81,andM=26.
(a)(b)(c)Figure 6
Contour plots of absolute error in rectangular area using different values ofα with N=225, n=81,M=26, and y=0.5.
(a)(b)(c)Figure 7
(a) Plot of absolute error withα=0.9, N=225, n=81,M=26, and y=−0.5. (b) Plot of absolute error with α=0.95, N=400, n=25,M=26, and y=−0.5.
(a)(b)


## 6.2. Problem 2

Consider equation (1) with exact solution given as [55](42)Ux,y,z,t=tα+3sinπxsinπysinπz,t>0,with wξ=1,Vξ=1,1,1T and the forcing term hξ¯,t is specified to satisfy the given equation and the exact solution. The initial and boundary data are extracted from the exact solution. The accuracy of the problem is investigated in cube domain 0,13 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 n=83. The numerical solutions in rectangular area for different values of α are presented in Figures 8(a)–8(c). The slice plots of exact and numerical solutions using α=0.85,N=2744,n=85,andM=26 are shown in Figures 9(a) and 9(b). Exact and numerical solutions with α=0.1,N=4096,n=85,M=28, and z=0.4,0.6,0.8 are shown in Figures 10(a) and 10(b), respectively. Absolute error in rectangular area using various values of α with N=1000,n=75,andM=28 is presented in Figures 11(a)–11(c). The contour slice plots of absolute error in rectangular area using various values of α with N=4096,n=85,andM=28 are shown in Figures 12(a)–12(c). The contour graphs of absolute error using various values of α with N=4096,n=85,M=28, and z=0.5 are presented shown in 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.Table 3
The computational results for problem 1 withN in the domain 0,13.
n=80,M=28N⟶12521672917283375α=0.1Error21.90 × 10−38.07 × 10−46.66 × 10−47.23 × 10−48.39 × 10−4Error33.20 × 10−37.92 × 10−47.89 × 10−45.38 × 10−46.61 × 10−4α=0.5Error21.90 × 10−37.98 × 10−46.60 × 10−47.11 × 10−48.22 × 10−4Error33.10 × 10−37.85 × 10−47.81 × 10−45.31 × 10−46.47 × 10−4α=0.9Error21.80 × 10−37.82 × 10−46.48 × 10−46.90 × 10−47.92 × 10−4Error33.00 × 10−37.71 × 10−47.65 × 10−45.19 × 10−46.22 × 10−4[55]————2.35 × 10−3—Table 4
The computational results for problem 1 withN in the domain 0,13.
n=83,N=4913M⟶2022242628α=0.1Error27.14 × 10−47.14 × 10−47.15 × 10−47.15 × 10−47.15 × 10−4Error35.86 × 10−45.87 × 10−45.87 × 10−45.87 × 10−45.87 × 10−4α=0.5Error26.98 × 10−46.99 × 10−47.00 × 10−47.00 × 10−47.00 × 10−4Error35.72 × 10−45.74 × 10−45.74 × 10−45.74 × 10−45.74 × 10−4α=0.9Error26.66 × 10−46.72 × 10−46.72 × 10−46.72 × 10−46.72 × 10−4Error35.44 × 10−45.49 × 10−45.50 × 10−45.50 × 10−45.50 × 10−4[55]————2.35 × 10−3—Figure 8
Approximate solution in rectangular area using different values ofα with N=1000, n=75,andM=28.
(a)(b)(c)Figure 9
Slice plots of exact and approximate solutions of problem 2 withα=0.85,N=2744,n=85,andM=26.
(a)(b)Figure 10
Exact and approximate solutions of problem 2 withα=0.1 and z=0.4,0.6,0.8.
(a)(b)Figure 11
Absolute error in rectangular area using different values ofα with N=1000, n=75,andM=28.
(a)(b)(c)Figure 12
Contour slice of absolute error in rectangular area using different values ofα with N=4096, n=85,andM=28.
(a)(b)(c)Figure 13
Contour plots of absolute error in rectangular area using different values ofα with N=256, n=85,M=28, and z=0.5.
(a)(b)(c)Figure 14
(a) Plot of absolute error withα=0.9, N=256, n=25,M=28, and z=0.2. (b) Plot of absolute error with α=0.95, N=625, n=38,M=28, and z=0.8.
(a)(b)


## 7. 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. The main advantage of this method is obtaining the solution at final time and avoiding the time stepping procedure. The LRBF produced sparse differentiation matrices due to which the ill-conditioning issues of the collocation matrices and the shape parameter’s sensitivity were resolved. The convergence of the method was discussed, and the bounds for the stability constant of the fully discrete system were derived. The accuracy of the method was tested by means of two examples, and the results were compared with those from other schemes. The 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.



---

*Source: 4640467-2021-10-08.xml*