A new approach to solving multi-order fractional equations using BEM and Chebyshev matrix

In this paper, the boundary element method is combined with Chebyshev operational matrix technique to solve two-dimensional multi-order time-fractional partial differential equations; nonlinear and linear in respect to spatial and temporal variables, respectively. Fractional derivatives are estimated by Caputo sense. Boundary element method is used to convert the main problem into a system of a multi-order fractional ordinary differential equation. Then, the produced system is approximated by Chebyshev operational matrix technique, ans its condition number is analyzed. Accuracy and efficiency of the proposed hybrid scheme are demonstrated by solving three different types of two-dimensional time fractional convection-diffusion equations numerically. The convergent rates are calculated for different meshing within the boundary element technique. Numerical results are given by graphs and tables for solutions and different type of error norms.


INTRODUCTION
Fractional order differential operators are the representative of non-local phenomena while many integer-order differential operators are mostly applied to examine local phenomena [1]; Therefore, fractional calculus can be useful to describe many of real-world problems which cannot be covered in the classic mathematical literature [2,3]. Since the next state of many systems depend on its current and historical states, there is a great demand to improve topical methods for the real life problems [4,5]. These problems happen in bioengineering [6], solid mechanics [7], anomalous transport [8], continuum and statistical mechanics [9], economics [10], relaxation electrochemistry [11], diffusion procedures [12], and complex networks [13,14], optimal control problems [15,16]. Fractional diffusion equations are largely used in describing abnormal convection phenomenon of liquid in medium. Models of convection-diffusion quantities play significant roles in many practical applications [4,5], especially those involving fluid flow and heat transfer, such as thermal pollution in river system, leaching of salts in soils for computational simulations, oil reservoir simulations, transport of mass and energy, and global weather production. Numerical methods for convection-diffusion equations described by derivatives with integer order have been studied extensively [17]. Due to the mathematical complexity, analytical solutions are very few and are restricted to the solution of simple fractional ordinary differential equations (ODEs) [18]. Several numerical techniques for solving fractional partial differential equations (PDEs) have been reported, such as variational iteration [19], Adomian decomposition [20], operational matrix of B-spline functions [21], operational matrix of Jacobi polynomials [12], Jacobi collocation [22], operational matrix of Chebyshev polynomials [23,24], Legendre collocation [25], pseudo-spectral [26], and operational matrix of Laguerre polynomials [27], Pade approximation and two-sided Laplace transformations [28]. Besides finite elements and finite differences [29], spectral methods are one of the three main methodologies for solving fractional differential equations on computers [30]. The main idea of spectral methods is to express the solution of the differential equation as a sum of basis functions and then to choose the coefficients in order to minimize the error between the numerical and exact solutions as well as possible. Therefore, high accuracy and ease of implementing are two of the main features which have encouraged many researchers to apply such methods to solve various types of fractional integral and differential equation. In this article shifted chebyshev polynomials are used. The boundary element method (BEM) always requires a fundamental solution to the original differential equation in order to avoid domain integrals in the formulation of the boundary integral equation. In addition, the nonhomogeneous and nonlinear terms are incorporated in the formulation by means of domain integrals. The basic idea of BEM is the transformation of the original differential equation into an equivalent integral equation only on the boundary, which has been widely applied in many areas of engineering such as fluid mechanics [31], magnetohydrodynamic [32], and electrodynamics [33]. In this paper the BEM is developed for the numerical solution of time fractional partial differential equations (TFPDEs) for nonhomogeneous bodies, which converts the main problem into a system of fractional ODE with initial conditions, described by an equation having a known fundamental solution. The proposed method introduces an additional unknown domain function, which represents the fictitious source function for the equivalent problem. This function is determined from a supplementary domain integral equation, which is converted to a boundary integral equation using a meshless technique based on global approximation by a radial basis function (RBF) series. The Delaunay graph mapping method can be viewed as a fast interpolation scheme. Despite its efficiency, the mesh quality for large deformation may deteriorate near the boundary, in particular, if the deformation involves large rotation, which may even lead to an invalid Delaunay graph. Furthermore, the RBF method can generally better preserve the mesh quality near the boundary but the computational cost is much higher, as the mesh size increases. In order to develop methods that are more efficient and because of their flexibility and simplicity, the Delaunay graph based RBF method (DG-RBF) were proposed [34] to overcome the difficulty of meshing and remeshing the entire structure. Thus, the pure boundary character of the method is maintained, since elements discretization and the integrations are limited only to the boundary. To obtain the fictitious source we use the Chebyshev spectral method based on operational matrix. The primary aim of this method is to propose a suitable way to approximate linear multi-order fractional ODEs with constant coefficients using a shifted Chebyshev Tau method, that guarantees an exponential convergence speed [35]. Once the fictitious source is established, the solution of the original problem can be calculated from the integral representation of the solution in the substituted problem. The outline of this paper is as follows. In Section 2, we introduce the multi-order time fractional convection diffusion equation (TFCDE) for a class of the TFPDEs as a mathematical modelling of natural phenomena, and some basic preliminaries are also given. Section 3 is devoted to applying the BEM for converting the main problem into a system of multi-order fractional ODE with initial conditions. In Section 4, the Chebyshev operational matrix (COM) of fractional derivative is obtained by applying the spectral methods to solve the generated multi-order fractional ODE. In order to demonstrate the efficiency and accuracy of the proposed method, along with the analysis of the condition number of COM, some numerical experiments are presented in Section 6 using the definitions and lemmas of Section 5. Eventually, we conclude the paper with some remarks, and add the appendix including notation table 7 to make more convenient understanding of the proposed algorithm.
In which is an integer number and is the Caputo fractional time derivative of order . The Caputo derivative [12], is employed because initial conditions having direct physical meaning can be prescribed. This derivative is defined as

IMPLEMENTATION OF BOUNDARY ELEMENT METHOD
Taking advantage of the following boundary element, the initial boundary value of the equation (1-3) is reformed into an ODE problem. Let ( , ) be the sought solution to the problem (1-3) and assume that is twice continuously differentiable in Ω. After applying Laplace operator on we have [31] where ( , ) known as an unknown fictitious source function. That is the solution of Eq. (1) could be established by solving Eq. (5) under the boundary condition (2), if convenient ( , ) is first established. This is accomplished by adhering to the following procedure.
We write the solution of Eq. (5) in the integral form. Thus, for as the normal derivative of we have [36] ( , ) = ∫ where * = ln ∕2 is the fundamental solution to Eq. (5), is the distance between any two points and also * stands for its normal derivative on the boundary. is the free term coefficient taking the values = 1 if ∈ Ω, = ∕2 if ∈ Γ, otherwise = 0; is the interior angle between the tangents of boundary at point . = 1∕2 for points where the boundary is smooth. After applying Eq. (6), to boundary points by means of Greens second identity [37], we yield the boundary integral equation wherê is a certain solution of the equation Also is the number of interior points inside Ω. Here are the coefficients that must be determined to satisfy is a set of radial basis approximating functions; are collocation points in Ω. The radial basis function method is used to map the nodes rather than that based on surface or volume ratios [34]. The algorithm is set out in the following procedure; At first we generate the Delaunay graph by using all the boundary nodes of the original mesh, and then we locate the mesh points in the graph, after that we move the Delaunay graph according to the specified geometric motion/deformation, and the final step is mapping the mesh points in the new graph according to the RBF matrix and Delaunay triangle. Procedures before the last step are exactly the same as the original Delaunay graph mapping method [38]; hence the details of these steps are skipped in this paper. The difference is in the last step, where the radial basis function interpolation is used to calculate the displacement of the internal mesh nodes from the given displacement of the Delaunay triangle nodes on the boundary, while the original Delaunay mapping method uses surface or volume ratios to calculate the displacement of inner nodes. Eq. (7) is solved numerically by using the BEM. The boundary integrals in Eq. (7) are approximated using boundary nodal points. Here = 1∕2, as we ace the smooth boundary.
The domain integral can be evaluated when the fictitious source is estimated by a radial basis function series and subsequently it is reformed to a boundary line integral using the GreenâĂŹs reciprocal identity [37]. For the sake of simplification, we use multiquadric radial basis function in practice.
internal nodals are here used to define Delaunay linear triangular elements in Ω. Therefore, after discretization and application of the boundary integral equation (7) at the boundary nodal points we have where and are × known as coefficient matrices originating from the integration of the kernel functions on the boundary elements and̄ is an × coefficient matrix originating from the integration of the kernel function on the domain elements; and are vectors containing the nodal values of the boundary displacements and their normal derivatives. Also, is the vector of the nodal values of the fictitious source at the internal nodal points. For a second order differential equation, the boundary condition is a correlation of 1 ( ) + 2 ( ) = ℎ( , ); after applying it at the boundary nodal points yields where 1 and 2 are × known diagonal matrices and ( ) = [ℎ 1 , , … , ℎ( , )] is a known boundary vector, where are boundary nodal points. Eqs. (9) and (10) can be combined and solved for and . This yields Further, differentiating Eq. (6) for points inside the domain ( = 1) with respect to and , using the same discretization and collocating at the internal nodal points, we have the following expression for the spatial derivativeŝ where thê is vector of values for and its derivatives at the internal nodal points;̂ and̂ are × known coefficient matrices originating from the integration of the kernel functions on the boundary elements and̂ is an × coefficient matrix originating from the integration of the kernel functions on the domain elements. Eliminating and from Eq. (12) using Eqs. (11) yieldŝ where The final step of the method is to apply Eq. (1) at the internal nodal points. This gives in which = and = for = 0, ..., . Now, from Eq. (13), we can write the initial conditions (3) for = 0, 1, .
where ( ) = 1 , , … , , . The above proposed procedure reduces the problem of multi-order two-dimensional time fractional PDE (1-3) to a simpler system of multi-term fractional ODE (16) with initial condition (18). The existence, uniqueness, and continuous dependence of the system (16)-(18) when S = 1 can be rigorously discussed (see e.g. Diethelm and Neville's paper [39]). In the next section, we show the implementation of Chebyshev operational matrix, as a spectral technique [30] for fractional calculus, to solve the system of initial value problem (16)-(18).

COM METHOD FOR SYSTEM OF MULTI-ORDER FRACTIONAL ODES
The Chebyshev polynomials ( ) are defined on the interval (−1, 1) [35]. Thus, by changing variable → 2 − 1, the shifted Chebyshev polynomials , ( ) of degree on the interval ∈ (0, ), with an orthogonality relation can be introduced by [30,40] , where , (0) = (−1) and , ( ) = 1. In this form, , ( ) may be generated by the following recurrence formula: where ,0 ( ) = 1 and ,1 ( ) = 2 ∕ − 1. Therefore, a given function ∈ 2 [0, 1] may be approximated by + 1 terms of shifted Chebyshev polynomials as where the coefficients are described by weight functions and suppose > 0 and the ceiling function ⌈ ⌉ denotes the smallest integer greater than or equal to , then where ( ) is the ( + 1) × ( + 1) COM of derivatives of order in the Caputo sense and is defined by [30,40]: where where 0 = 2, = 1, ≥ 1. Note that in ( ) , the first ⌈ ⌉ rows are all zero. In order to solve Eq. (16) with initial conditions (18), we approximate ( ) and ( ) in terms of shifted Chebyshev polynomials as where for = 1, ..., and Ψ Ψ Ψ, are × ( + 1) matrices that are defined as (21) and Eq. (23) can be used to write Employing Eqs. (23), (24) and (25), then the residual for Eq. (16) can be written as ( ) is a vector with respect to . If ( ) be the th component of ( ), then in a typical Tau method [35], we generate ( − + 1) linear equations with ( + 1) unknown coefficients of Ψ Ψ Ψ by applying Also, by substituting Eq. (23) into Eq. (18), and with the fact that we get a system of linear equations with ( + 1) unknown coefficients for Ψ Ψ Ψ as following Equations (27) and (28) can be rewritten in the matrix form where is an ( + 1) × ( + 1) coefficient matrix. The system of algebraic equations (29) can be easily solved for the unknown vector Ψ Ψ Ψ. Consequently, ( ) given in Eq. (23) can be calculated, which gives a solution of Eq. (16) with the initial conditions (18). Once the vector ( ) of the values of the fictitious source at the internal nodal points has been established, then the solution of Eq. (1) and its derivatives can be computed from Eq. (13). For the points = ( , ) that do not coincide with the prespecified internal nodal points, the solution could be drawn from the discretized counterpart of Eq. (6) with = 1 using the same boundary and new domain discretization. Note that here the matriceŝ ,̂ and̂ corresponding to previous internal nodes plus the additional points must be recomputed.

ERROR ESTIMATION
The convergence of the proposed method is shown by employing the following error norms, maximum error ( ∞ ), maximum relative error ( ) to assess the accuracy of the method in multi-scale problems, and root mean square ( ) to globally examine the method efficiency, where and denote the th components of the exact and approximated solutions, respectively, and denotes the number of internal points. It is not convenient to certainly determine what is the convergence rate of the proposed hybrid method; for example, for the number of nodal points and , and the size of the shifted Chebyshev polynomials , the convergence order of the method would be ( ( , )) where is a function of the convergence rate of BEM with the order and the convergence rate of COM with the order . The accuracy of the method depends on several factors, the convergence speed for BEM, the domain and boundary discretization, the shape parameters of radial basis functions, the orders of the derivatives, and condition number of COMs, as such the analysis of truncation errors for methods solving a two dimensional multi-term time fractional differential equation is not straightforward. Nonetheless, information about matrix from the algebraic system (29), particularly its condition number, will be useful. The condition number is defined by [35] such that a matrix with a large condition number is so-called ill conditioned, whereas the matrix is named well conditioned if its condition number is of a moderate size. We also suggest two -orders in the following lemmas to examine the rate of convergence for BEM and COM distinctly. The first one is directly tested by the exact solution and the effect of domain-discreitization. While the second one is addressed by comparing a sequence of numerical solutions of the ODE system (16) with different degree sizes of COMs which have been offered exponential rates of convergence accuracy for smooth problems in simple geometries [35]. Proof. When the leading terms in the spatial-discretization error are proportional to 1 and 2 , and ‖.‖ denoting the root mean square norm (32), the temporal convergence order for the COM method presented in Section 4 is estimated using -order≃ where and denote the th components, and 1 = 1∕ 1 , 2 = 1∕ 2 , and 3 = 1 3 .
Proof. When the leading terms in the error of COM are proportional to 1 , 2 and 3 , In the following section, the numerical errors are computed based on assumptions described in Lemma 5.1 and 5.2.

NUMERICAL RESULTS AND DISCUSSION
On the basis of the described procedure, some problems are solved to illustrate the efficiency and the accuracy of the proposed method. In the first example, a simple two-dimensional fractional heat-like equation is considered for two different conditions. In the second example, a nonlinear two-dimensional fractional wave-like equation is tested. In the third and fourth problems, two linear TFCDEs are solved to test the impact of external force ( ) and the final time on the convergence rate of the method. In the fifth and sixth test problems, multi-order time-fractional diffusion-wave equations in bounded homogeneous anisotropic plane bodies are solved. The condition number of system 29 is examined for each example. Since, the size of the matrix depends on the number of internal points, , and the degree of COM, , the condition number of can be compared versus and the length of distance between nodal points. However, most domains are not discretized uniformly. In this regard, suppose denotes the mean length of all the distances between the internal points and their adjacent points (e.g. see in Figure 8 ). Thus, the numerical results show that the condition number behaves as ≃ −2 ( + 1) 2 for example 6.1, and ≃ −2 ( + 1) 3 for other examples.
where the following one-parameter Mittag-Leffler function is defined as , > 0. Figure 1 demonstrates ∞ errors and versus the degree (right) for Example 6.1 (case I) when = 160, and = 0.5 at = 1.5. The convergence rate of COM is estimated that -order> 5. Furthermore, the condition numbers of , on the figure 1 (left), are shown versus the polynomial degree, , and the mean length of discrete elements, . It can be numerically deduced that the condition number behaves as ( ) ≃ −2 ( + 1) 2 ; e.g. when = 8, then ( ) ≃ 76.4 × −2 , and when = 0.3, then ( ) ≃ 10.43 × 2 . In Table 1 , numerical results are compared with the exact solutions (34) for Example 6.1, case , for fixed = 12, = 1.5, with the differential orders = 0.5 and = 0.75, and different number of nodal points; the convergence rate of BEM is algebraic ( -order> 4.4) when the number of nodal points is increased from = 40 to = 80 and it is quadratic ( -order> 2) when = 80 goes to = 160. Apart from the value of , it can be inferred that the computation cost of the second discretization for moderate and is more effective than the third one. For case (II), a similar behavior of ( ) versus and is shown in Figure 2 (left). The relative absolute error (right) with = 200, and = 12 for the final time = 1 is exhibited. Intuitively, the relative absolute errors are approached to 10 −5 , which it could be expected for = 12 and = 200 based on the information from 2 . This table shows the estimated convergence for two terms of shifted Chebyshev series for the final time = 0.5; in general, there is an improvement for errors when the degree increases, but no relationship between -order and degree is observed. It may also be concluded = 64 is more computational cost-effective in this case.     The exact solution for = 2 is found to be, This problem is solved for different with = 100 at = 0.5 and integer order = 2 to compared with the exact solution (35). The results from Table 3 offer an error improvement by increasing the number of degree of COMs. The condition numbers of the system 29 illustrate such behaviour ( ) ≃ −2 ( + 1) 3 (see Table 3 ). Due to the fact that the domain and boundary nodal points are fixed here, the numerical solutions of U are directly affected by the numerical solutions of b; in other words, affected by the accuracy of COM. However, the exact solution of the generated ODE system is not clear, and the convergence order of COMs is estimated by norm ‖.‖ for three distinct degrees with a same proportion. Interestingly, a comparison between the two columns and ‖.‖ of Table 3 suggests direct relationships, but with different speed between the approximation solution of U and b. In addition, by considering the scale of the solutions, and a comprehensive assessment of the absolute errors for distinct degrees , it could be concluded that the Chebyshev Tau method converges with an oscillating manner around the exact solution of fractional ODE system (16).   Table 4 gives ( ), the RMS error and the convergence rates are obtained by solving Example 6.3 for different values of . It indicates better errors for near to 1 than 2, that is not true for -orders.   = 150 would be sufficient to achieve the semi-equivalent errors. Importantly, by considering the results of the previous example 6.3 and Figure 4 (right plan), it may convey that the method has a better performance for the less values of . In contrast, Table 5 refuses this idea; there are irrelevant outcomes versus the  values of , although the table shows a reliable numerical convergence.   TABLE 5 The  error for the Table 6 showing the efficiency of the proposed method by -order > 2, and the condition number of matrix with the manner as −2 × ( + 1) 3 . In figure 5 , the contour plots illustrate the absolute errors distributions of the approximations of , and on the plane, including ∞ , , , at = 4 for specific nodal points and Delaunay triangulation when = 210, = 132, = 16.
in a "C-shape" made by the elimination of a circle with radius 2 = 3 and null origin from the inside of a circle with radius 1 = 5 and the same origin, and extracting the space between the lines = −1 and = 1 from the right side of the outcome (see , to find the exact solution as = ( ) (x). This problem is solved with = 216, = 741, and = 16, for = 5. In Figure 7 , the distribution of the absolute errors on the domain, ∞ , , and of , , and are illustrated. Figure 6 exhibits the behavior of the condition number matrix as −2 × ( + 1) 3 ; e.g. when = 1, ( ) ≃ 0.92 ×( + 1) 3 , and when = 10, ( ) ≃ 1202.46 × −2 . Among these six examples, an interesting point can be concluded that the error distribution into a plan depends not only on the positions of the nodal points but also on the final time solving the problem; with an increasingly asymmetric discretization, and longer computing time, the error distribution becomes more "random."

FIGURE 6
The condition number of matrix versus the polynomial degree (left) and versus the mean length of discrete element edges (right) for Example 6.6 when = 5.

CONCLUSION
Here, we have proposed a hybrid algorithm to solve two-dimensional multi-order time-fractional partial differential equations. Their general form is given in equations (1-3). The method consists of the boundary element method combined with spectral Chebyshev operational matrix. The BEM is used to transfer the corresponding time fractional PDE into a system of ODEs while COM is used to solve the system efficiently. This method is applied to the two-dimensional fractional heat-like, wavelike and diffusion-wave equations, which shows that the errors of the approximate solution decay exponentially. When the exact solution exists, comparison is made with , and the convergence rate is calculated using Lemma 5.1. When the exact solution is not in hand, the order of convergence is estimated by three approximate solutions with various degrees of Chebyshev polynomials in a same grid-point based on Lemma 5.2. By applying the assumptions of the Lemmas, the numerical results show the efficiency and convergence rate for the proposed hybrid method. Notwithstanding, it is not easy to emphasize a unique conclusion for the accuracy of the method on the ground that given the vast range of architectures, spectral methods, boundary element methods, fractional calculus, and meshing used with in such a hybrid-technique framework. In general, for multi-order

FIGURE 7
The geometry of the plane body discretization with = 216, = 741 with DG-RBF technique and Contour plots of absolute errors when = 16 for Example 6.6. two-dimensional time fractional PDE (1-3), the current method calculate for the test examples, with this range of convergence rate: 2 < -order< 4.5 for the moderate values of and . And for multi-term fractional ODE (16) with initial condition (18) the current method works well to calculate solution with a range of the convergence rate around 1 < -order< 5.5. Moreover, The condition number of matrix from linear system (29) behaves like ≃ −2 ( + 1) 2 for the problems with ⌈ ⌉ = 1, and ≃ −2 ( + 1) 3 for the problems with ⌈ ⌉ = 2. For the future direction, the authors believe that establishing new methods to examine long-term effects of memory in complex systems modeled by fractional calculus is highly required as fractional calculus is a proper mathematical tool for describing memory [14], while the proposed COM technique is not an appropriate scheme for long-term problems. It is instead efficient for problems with multi-term orders.

Algorithm in a nutshell
For the convenient, a notation table, and the proposed algorithm's description is given in this section. At the beginning, determine boundary points and divide internal nodal points into groups by the Delaunay graph, then use the RBF method to interpolate the nodal points to its new position. Hence, there is no need to optimize shape parameters for the domain discretization. After discretization, consider ∫ indicating integration on k-element on the boundary (see Figure 8 and Table 7 ), and set , = 1, ..., , = 1, ..., for the boundary points , the domain points , and the algorithm implementation. Notice, different must be considered for computing the boundary integrals at the corner points, particularly in inhomogeneous shapes, in comparing to smooth boundary [31].
The location of nodal points and relative distances for constant element discretization.

Step 14
Output the corresponding solution of the problem (1) and its derivatives by applying equation (13). Unknown field function of spatial x ∈ Ω ∪ Γ, and time Normal derivative of , , , Given coefficient functions of x and = 0, 1, ..., , , , g Given independent function of x and t B Linear operator with respect to x of order one ℎ(x, ) Given function in the boundary condition; x ∈ Γ (x) Given function in the initial condition; = 0, 1, ..., − 1 (x, ) Unknown fictitious source function * , * Fundamental solution of (5), and normal derivative of * on the boundarŷ ,̂ Particular solution of (9), and normal derivative of̂ Free term coefficient; = 1 if x ∈ Ω, = ∕2 if x ∈ Γ, else = 0, see details in book [31] Interior angle between the tangents of boundary at point x, see details in book [31] Distance between two points (or mean of all distances of internal points) Number of interior points after discretization internal nodal points; = 1, ..., Number of boundary nodal points after discretization boundary nodal points; = 1, ..., condition number of matrix of the algebraic system equations (29)