{'fraction_non_alphanumeric': 0.055917334912006886, 'fraction_numerical': 0.0192131747483989, 'mean_word_length': 4.161666666666667, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical solution of the general high-dimensional multi-term time-space-fractional diffusion equations

In this article, an advanced differential quadrature (DQ) approach is proposed to obtain the numerical solutions of the two- and three-dimensional multi-term time-space-fractional diffusion equation (TSFDE) on general domains. The fractional derivatives in space are firstly discretized by deriving a new class of differential quadrature (DQ) formulas with radial basis functions (RBFs) being the trial functions. Then, the original problems are converted to a group of multi-term fractional ordinary differential equations (ODEs), which are further treated by a class of high-order difference schemes based on the weighted and shifted Lubich difference operators. The presented DQ method is high accurate and convergent on arbitrary domains. Finally, several numerical examples are carried out to illustrate its accuracy and effectiveness.


Introduction
During the last ten years, fractional calculus came into the focus of interest in academic circles as an important branch of mathematics. Due to its non-locality and self-similarity, it has been proven to be very adequate to describe the dynamics phenomena with long term memory memory and hereditary effects like the anomalous transport in complex heterogeneous aquifer [1]. The application of fractional calculus has been very extensive recently, including biology, chemistry, physics, economics, and many other disciplines.
The fractional partial differential equations (PDEs) are the results of mathematical modelling based on fractional calculus, which provides a new powerful tool for the study of mathematical physics or even for the whole scientific research. Nevertheless, there raises a challenge to solve these type of equations and few of fractional PDEs can be solved by analytic techniques. As an result, numerical methods are clearly a priority for development. Until now, the numerical algorithms for the time-fractional PDEs is on the way to maturity after recent years' development [26,33,23,43,39,46,19,20], while these for high-dimensional space-fractional PDEs have to be further developed and advanced, especially for the ones defined on general domains. Regardless of the difficulties in constructing numerical algorithms for the space-fractional PDEs, many numerical methods have been invented to solve them, covering finite difference (FD) methods [32,40,21,22], spectral methods [5,44], finite element (FE) methods [15,13,45,47], and finite volume methods [44,29]. Liu et. al considered a space-fractional FitzHugh-Nagumo monodomain model by an implicit semi-alternative direction FD scheme on approximate irregular domains [28]. Qiu et al. considered a nodal discontinuous Galerkin method and solved a L-shaped domain problem [38]. Yang et.al proposed a fully discrete FE method for the two-dimensional space-fractional diffusion equation on convex domains [42]. Bhrawy and Zaky developed a spectral tau method for the multi-term time-space fractional PDEs [6]. In [11], an efficient algorithm based on FD and FE methods was addressed for the two-dimensional multi-term time-space-fractional Bloch-Torrey equations. In [16], a fully discrete FE method was given for the two-dimensional multi-term time-space fractional diffusion-wave equation. Although the numerical treatment of multi-term time-space-fractional PDEs is also an active area of research, few works can be found to solve these equations on general high-dimensional domains.
Meshless methods eliminate tedious mesh generation and reconstruction or only use easily generable meshes in a flexible manner, and can reduce computational cost, which serve as a promising alternative in dealing with structure destruction, high-dimensional crack propagation, and large deformation problems. Meshless methods have been achieved great progress, such as diffuse element method [34], reproducing kernel particle method [18], hpcloud method, element-free Galerkin method [4], meshless local Petrov-Galerkin method [2], boundary element-free method, RBF collocation methods (Kansa's methods) [17,24], point interpolation (PI) method [30], DQ methods [14,41], and so forth. They offer more advantages over mesh-dependent methods in treating the space-fractional PDEs but seldom works have bee reported. Liu et. al proposed a PI method for the space-fractional diffusion equation [31]. Cheng et. al analysed an improved moving least-squares approximation to the solution of two-dimensional space-fractional wave equation [10]. Pang et. al extended DQ method to solve the space-fractional diffusion equation and use the as Lagrange interpolating basis functions as trial functions to determine the weighted coefficients [36].
DQ method is a kind of meshless methods which uses the weighted sum of functional values at sampling points along ordinate directions to discretize derivatives [3]. Subsequently, this kind of method enjoys the advantages as high accuracy, truly mesh-free, and the adaptability to high-dimensional problems. Regarding this point, in this article, we showcase an efficient RBFs-based DQ method for the following multi-term TSFDE on general domains: (I) two-dimensional multi-term TSFDE: where 0 < θ 1 , θ 2 , . . . , θ s < 1, s ∈ N, 1 < α, β < 2, Ω ⊂ R 2 with ∂Ω being its boundary, ε ± ς (x, y) are the diffusion coefficients with ς = α, β.
The fractional derivatives in Eqs. (1.1)-(1.2) are defined in Caputo sense. More precisely, the ones tagged with "+" are the left-side derivatives, while those tagged with "-" are the right-side derivatives. For example, assuming C 1 : x = x L (y, z), C 2 : x = x R (y, z) the left and right boundaries of Ω ⊂ R 3 , in which, are defined as follows: with the Gamma function Γ(·).
The fractional derivatives in other coordinate directions can be defined in the same fashion. P θ 1 ,θ 2 ,...,θ s C 0 D t denotes the multi-term fractional derivative operator: The outline is as follows. In Section 2, some preliminaries on fractional calculus and RBFs are introduced. In Section 3 the DQ approximation of fractional derivative is proposed based on RBFs. In Section 4, using the RBFsbased DQ formulas, we construct a fully discrete DQ method for the multi-term TSFDE on general domains and its algorithm is further studied in Section 5. In Section 6, some numerical tests are presented to confirm its convergence. In the last section, a brief summary is drawn.

Preliminaries
We recall some basic preliminaries on fractional calculus and RBFs required for further discussions.

Fractional calculus
where n − 1 < α < n, n ∈ Z + , and if α = n, C a D n x u(x) = u (n) (x) and C x D n b u(x) = (−1) n u (n) (x). 3 The two frequently-used fractional derivatives are equivalent with exactness to an additive factor: Moreover, we easily realize the following properties: where C is a constant and l > [α] + 1 with [x] being the ceiling function, which outputs the smallest integer greater than or equal to x. Another point need to be noticed is that these properties are also true for the right-side fractional calculus. For more details, we refer the readers to [25] for overview.

Radial basis functions
∈ Ω ⊂ R d and || ·|| be Euclidean norm, the approximation of u(x) can be written as a weighted sum of RBFs in the form: To ensure that the interpolant properly behave at infinity, the above equations are augmented by (2.11) Putting Eqs. (2.10)-(2.11)in matrix form leads to a matrix system: where the matrix elements for The RBFs provides an effective tool for the interpolation of scattered data in multi-dimensional domains. Since they only take the space distance as independent variable, as compared to the traditional basis function defined along coordinates, RBFs have the advantages of concise form, dimension-independent and isotropic properties. In particular, in the application of solving high-dimensional PDEs, the RBFs interpolation shows great flexibility. The numerical methods constructed by RBFs not only have the simplicity in the implementation, but also avoid the tedious mesh generation, and thereby reducing the pre-processing process of algorithms. If r = ||x − x k ||, then the commonly used RBFs are listed in Table 1.

Defition 2.4. A function is completely monotonic if and only if
. ., and r ≥ 0.

The DQ approximation of fractional derivative
In this section, we derive a new class of DQ formulas for fractional derivatives by RBFs as the trial functions. The DQ formula employs the linear sum of function values at discrete points to approximate partial derivatives with regard to an independent variable. Supposing u(x) ∈ C k (Ω), k ∈ Z + , we have the following DQ formula [41]: where w m i j , i, j = 1, 2, . . . , M, are weighted coefficients. By realizing the essence of DQ formula, we propose the following DQ approximations to the left-hand fractional derivaitves: where . . , M, are the related weighted coefficients, which satisfy being the trial functions. The weighted coefficients are determined by reforming Eqs. (3.15) in matrix system in advance. Likewise, we can extend the above formulas to the right-side fractional derivatives with regard to the variables x, y and z.
The calculation of weighted coefficients comprise a major part of DQ methods. In what follows, we show how to use the aforementioned RBFs to determine these coefficients. It should be noted that for Reciprocal Multiquadrics, Inverse Quadratics and Gaussians, their interpolation matrices invertible. Therefore, ignoring the polynomial terms on the right-side of RBFs interpolation (2.10) and substitute φ k (x i ) = ϕ k (x i ) into Eq. (3.15), we obtain the weighted coefficients by solving the resulting algebraic systems for each nodal point x i . As for Multiquadrics, we prefer to (3.16) To make the problem be well posed, one more equation is required. From Eq. (2.11), it follows that µ 1 (t) = M k=0 λ k (t), which further yields By careful observation, we can consider . . , M, as the trial functions. Referring to the above discussion and C a D θ x 1 = 0, we obtain the algebraic system as follows For each node x i , the values of unknown weighted coefficients are consequently computed by w α + i = M −1 R i . We note that the DQ approximation (3.14) works because the fractional derivatives are linear operators. In other x ∈ Ω, and the weighted coefficients fulfill Eqs. (3.15), then for through the linear property of fractional derivatives. In addition, it is sure that the above convergent result is exactly true for the right-side fractional derivatives with regard to the others variables.

Discretization of time-fractional derivative
In this subsection, we introduce high-order schemes to discretize the fractional derivative in temporal direction. To this end, we define a time lattice t n = nτ, n = 0, 1, . . . , N, T = τN, N ∈ Z + , and review the relationship between Caputo and Riemann-Liouville derivatives: where m − 1 < θ < m, m ∈ Z + . Using the property (2.9), we obtain a class of high-order schemes for Caputo derivaitve by applying the weighted and shifted Lubich difference operators L θ q u(x, t n ) = 1 τ θ n k=0 ω q,θ k u(x, t n−k ), q = 1, 2, 3, 4, 5, to discretize Riemann-Liouville derivatives on the right-side of Eq. (4.19), which reads where {ω q,θ k } n k=0 are the discrete coefficients [8]. For instance, if q = 1, we have ω 1,θ where µ = 4 , and i denotes imaginary unit.
Lemma 4.1. The coefficients ω 1,θ k satisfy the following properties Proof. The proof is trivial by following Lemma 4.2 and references [8,12].
Applying the operator A θ q to the multi-term fractional derivative operator and set 0 < θ < 1, we finally obtain (4.24)
In the sequel, we employ the notations U n t n ), and U n , F n , g n are the column vectors consisting of U n i , f n i , g n i in the ascending order of subscript i, respectively. Also, we adopt ε ± α l = diag(ε ± α l ,0 , ε ± α l ,1 , . . . , ε ± α l ,M ). Then, reforming Eqs. (4.28) in a matrix system leads to where W ± α l are the weighted coefficient matrices, given by Handling the boundary conditions properly and solving the algebraic equations resulting form the scheme (4.29) at each time step, we thus obtain the desirable solutions for the considered problem on all time layers.

Algorithm and implementation
In this section, we turn our eyes to the implementation of the proposed scheme (4.29). To perform the algorithm, an important problem having to be addressed is how to compute the fractional derivatives of RBFs and determine the corresponding integral paths. It is acknowledged that the fractional derivative of a general function as RBFs is difficult to compute. Hereafter, we will be engaged in this point. For a node x k ∈ Ω and the RBF ϕ k (x), letting ξ = x l − x l − X l,L (1 + ζ)/2, one has where X l,L = min η l : (η 1 , η 2 , . . . , η d ), η i = x i , i l . Then, it can be treated by Gauss-Jacobi quadrature formula [35]. The line segment along the x l -coordinate axis from the left-side of boundary to the node x l is called the integral path, which need to be determined before the calculation of fractional derivatives.
Hereafter, we will carefully illustrate how to perform the RBFs-based DQ scheme. Using Matlab as a platform for programming, some common commands in Matlab will be used as symbols. Let nonb, boup be the index vectors of the internal and boundary nodes of Ω, respectively. For the sake of simplicity, we define U n (nonb), , boup),F n = F n (nonb) andg n = g n (nonb). Imposing the boundary condition, we obtain the following RBFs-based DQ scheme s r=1 a r ω q,θ r 0 Compute the fractional derivatives of RBFs and get W ± α l (i + 1, : Generate the discrete coefficients ω q,θ l j , j = 0, 1, 2, . . . Set U(nonb, n + 1) =Ũ n and U(boun, n + 1) =g n .
end for Return U

Numerical exmaples
Let u n j , U n j be the exact and numerical solutions at the time level n. Denoting some illustrative examples will be performed to reveal the validity and efficiency of the proposed RBFs-based DQ method. In the runs, we choose the values of shape parameters with reference to the formula c = ν/(M + 1) σ/4 or artificially choose their values, where ν > 0, σ ∈ Z + . Besides, we solve all one-dimensional problems by using the nodal distribution 6.1. One-dimensional problems Example 6.1. To test the validity of the DQ formula, we discretize the left-side fractional derivative C −1 D α x (x + 1) 3 on [−1, 1] by Grünwald-Letnikov (GL) difference operator and DQ formula, in which case, one has C −1 D α x (x + 1) 3 = 6(x + 1) 3−α Γ(4 − α) .
We choose α = 1.5 and Reciprocal Multiquadrics as the trial functions. The comparisons of numerical errors and convergent rates (Cov.rate) of these two methods are given in Table 2. It can be seen that the convergence rates of GL-operator is strictly 1, while the convergence rates of DQ formula is much higher than GL-operator, and the errors of DQ formula are less than GL-operator. Thus, our DQ formulas are effective and more efficient than GL-operator.
The exact solution is u(x, y, t) = (1 + t 2 )x 2 (1 − x) 2 y 2 (1 − y) 2 . To show the superiority of the method as much as possible, we compare the spatial accuracy of the weighted shifted GL methods (WSGL) [40] ((p, q) = (1, 0) for WSGL1 and (p, q) = (1, −1) for WSGL2), FD and FE methods [7,37], and our DQ methods based on Reciprocal Multiquadrics in term of ||u − U|| L 2 . Taking τ = 1.0 × 10 −3 , ν = 8 and σ = 2, the numerical errors at t = 1 of all above methods are presented in Table 5. In Fig. ??, the used configuration of 289 nodes and its corresponding error distribution of RMQ-based DQ method are presented. We observe that our method is more accurate than WSGL and FD method by using the same discrete parameters. Although the error magnitude of DQ method is about the same as FE method, it is sure that our method is more flexible in implementation and has a less amount of computation.
The exact solution is u(x, y, t) = t 3 x 2 y 2 . We choose Multiquadrics and Gaussian as the trial functions and report the numerical results at t = 0.5 in Table 6 by using τ = 1.0 × 10 −3 . In Fig. ??, we show the used configuration of 152 nodes and its corresponding error distribution of GA-based DQ method. It can be seen that the proposed DQ methods are convergent to the exact solution as the nodal number increases.
u(x, y, 0) = (4x 2 + y 2 − 1) 2 /10, (x, y) ∈ Ω, u(x, y, t) = 0, (x, y; t) ∈ ∂Ω × (0, T ], (6.36) on the elliptical domain Ω = {(x, y)|4x 2 + y 2 < 1} with a 1 = a 2 = a 3 = 1, θ 1 = 0.6, θ 2 = 0.7, θ 3 = 1, α = β = 1.6 and the source term where x l = − 1 − y 2 /2, x r = 1 − y 2 /2, y l = − √ 1 − 4x 2 and y r = √ 1 − 4x 2 . The exact solution is u(x, y, t) = 1+t 2 10 (4x 2 + y 2 − 1) 2 . We compare the spatial accuracy of FE method and our DQ method base on Inverse Quadratics and Gaussians in term of ||u − U|| L 2 . The numerical errors at t = 1 of FE method [27] and DQ methods are presented in Table 5 by using q = 2 and τ = 1.0 × 10 −3 . The used configurations of 112, 238 and 412 nodes, and their corresponding error distributions of IQ-based DQ method are displayed in Fig. ??. It is indicated that DQ methods produce the same error magnitude as FE method with less nodal numbers, which confirms the effectiveness of our methods.  We choose τ = 1.0 × 10 −3 and report the numerical results at t = 0.5 in Table 8, where Multiquadrics and Reciprocal Multiquadrics are used as the trial functions. Also, we highlight the used configuration of 152 nodes and its corresponding error distribution of MQ-based DQ method in Fig. ??. The results observed illustrate that the proposed DQ methods yield the numerical solutions which are in excellent agreement with the exact solution. Moreover, the errors produced by these two DQ methods do not differ significantly under the selected discrete parameters.

Conclusion
In this investigation, an advanced DQ technique was developed to obtain the numerical solutions of the multi-term TSFDEs on two-and three-dimensional general domains. The considered problem was reduced into a group of multiterm fractional ODEs by extending the DQ approximation to the fractional derivatives with the help of some familiar RBFs as the trial functions. A fully discrete scheme has been derived by applying a class of high-order difference schemes to discretizing the resulting multi-term fractional ODEs. The proposed DQ scheme produces very satisfactory accuracy both in time and space, and the numerical reuslts and comparisons with the exact solutions or the errors 13