Abstract
In this article, radial basis function collocation scheme is adopted for the numerical solution of fractional partial differential equations. This method is highly demanding because of its meshless nature and ease of implementation in high dimensions and complex geometries. Time derivative is approximated by Caputo derivative for the values of and . Forward difference scheme is applied to approximate the 1st order derivative appearing in the definition of Caputo derivative for , whereas central difference scheme is used for the 2nd order derivative in the definition of Caputo derivative for . Numerical problems are given to judge the behaviour of the proposed method for both the cases of . Error norms are used to asses the accuracy of the method. Both uniform and nonuniform nodes are considered. Numerical simulation is carried out for irregular domain as well. Results are also compared with the existing methods in the literature.
1. Introduction
Fractional order calculus is a dynamic branch of calculus, which is concerned with the integration and differentiation of noninteger order. This branch of mathematics attracted researchers in the last few decades [1–7]. Fractional partial differential equations (FPDEs) are commonly used to model problems in the field of science, engineering, and many other fields including fluid mechanics, chemistry, viscoelasticity, finance, and physics. Some interesting applications of FPDEs can be found in [8–11]. Many researchers did considerable work to find the analytic solution of FPDEs [12–14], but it is difficult and sometimes impossible to find the analytic solution of most of FPDEs. Therefore, many researchers referred to numerical techniques to find the solution of FPDEs [15–20]. There are two widely used definitions of fractional derivatives, namely, Caputo and Riemann-Liouville. The main difference between these two operators is the order of evaluation. Many authors analyzed time-fractional partial differential equations (PDEs), for example, Wyss [14], Agrawal [21 ], Liu et al. [22], Jiang and Ma. [23], and Chang et al. [24].
Radial basis function (RBF) method has been used to find the solution of FPDEs. RBF collocation schemes are used to find the solution of PDEs, integral equations, integrodifferential equation, etc. The main idea behind the RBF method is to approximate space derivatives by RBFs which converts PDE to a system of linear equations. The solution of this system of linear equations leads to the solution of governing equation. This method is getting fame due to its meshless nature and easy to use in high dimensions and complex geometries. To utilize this advantage of RBF scheme, it is applied to time-fractional PDEs for higher dimensions and with different types of domains. In this paper, implicit scheme (IS) and Crank-Nicolson scheme (CNS) are coupled with RBF. Many authors used meshless RBF method to solve FPDEs [25–35]. In this paper, multiquadric (MQ) RBF is used to approximate solution. MQ-RBF is defined bywhere , , is the number of collocation points and is the shape parameter. Furthermore, in one dimensional, and in two-dimensional case.
Kansa [36] has applied the multiquadric radial basis function (MQ-RBF) collocation method to solve PDEs. After that, there are a lot of applications and developments of the MQ-RBF as an efficient meshless method to solve engineering problems. However, the ill-conditioned behaviour and the sensitivity to the shape parameter are the main obstacles in the Kansa’s MQ-RBF method. Many researchers have discussed the optimal shape parameter used in the MQ-RBF [37–39]. Formulation of the method flows in the following major steps:(1)Approximate time-fractional derivative by using Caputo definition(2)Approximate the space variable by RBFs(3)Substitute the values obtained from the previous two steps in the problem to get a system of linear equations
Definition 1. Caputo derivative of noninteger orderof a functionis defined by [9]
In this paper, we have tackled the following two cases of Caputo derivative:
Case I.
Case II.
1.1. Governing Equation
We emphasize on the following time-fractional PDEThe boundary conditions (BCs) areand the initial conditions (ICs) arewhere is the solution, is the Caputo fractional derivative of order , is the source term, is the bounded domain, and is the boundary. Equation (3) can be fractional diffusion, fractional wave diffusion, or fraction anomalous diffusion equation depending on . Here, is a linear operator defined bywhere , , and are functions of or constants and and denote Laplacian and gradient operator, respectively.
1.2. Main Objective of the Paper
This paper is aimed at solving FPDEs by using the combination of Caputo fractional derivative operator and RBFs. Caputo fractional derivative operator is applied to approximate the time derivative whereas RBFs are adopted to approximate the space derivatives. The organization of the rest of the paper is as follows: Section 2 is dedicated to construct meshless scheme for the first case, i.e., . In Section 3, we consider the second case, i.e., . In section 4, the numerical method is applied to different problems and comparison is made with some other methods. Section 5 is authoritative to give the concluding note of this work.
2. Formulation of the Method for Case I
In this section, we take and find , by using Caputo derivative. Finite difference scheme is applied to approximate the 1st order time derivative appearing on the right-hand side of Caputo derivative. Then, -weighted scheme is applied to the governing equation, and the value of the time derivative is also substituted.
2.1. Time-Fractional Derivative
Caputo fractional derivative for is defined by
Taking derivative at , we getThis implies
Now, we use finite difference scheme to approximate , as follows:where is the time step size.
Thenwhere and , .
Finally, we can write in more precise form as
2.2. -Weighted Scheme
Applying -weighted scheme to Eq. (3) and substituting Eq. (12) in Eq. (3), we get
Now, we employ Kansa’s method and interpolate by RBFs. We interpolate the solution at different points where are interior points while and are boundary points.where are the RBFs, is Euclidian norm, and s are the unknown constants.
We can write Eq. (14) in matrix form asor
Provided that collocation matrix must be nonsingular to ensure the invertibility of matrix . This depends on the choice of RBF and the location of mesh points. Matrix is invertible for distinct mesh points. The shape parameter has an important effect on condition number [40]. Once we find the constants , we can find solution from Eq. (14).
Putting the value from Eq. (15) in Eq. (13), we get the following formwhere
Here,where and are some known functions given in BCs. Finally, from Eq. (18) and Eq. (16), we get
We can use this scheme to find the solution at any time level .
2.3. Stability and Convergence
The scheme (21) is basically a recurrence relation which gives the values of solution at the time by using solution at time . The matrix is known as amplification matrix, and the elements of this matrix depend on the constant , where is the time step, is the distance between any two successive nodes, and is the order of spatial differential operator. Let us denote the exact solution of Eq. (3) by at time .
We state a well-known theorem of Fasshauer [41] (see [42] for proof):provided that . Here,
Theorem 2. [41] Suppose is open and bounded and satisfies an interior cone condition. Suppose is symmetric and strictly conditionally positive definite of order on . Denote the interpolant to on the -unisolvent set by. Fix with . Then, there exist positive constants and (independent of , and ) such that
Application of Theorem 2 to infinitely smooth functions such as Gaussians or generalized (inverse) multiquadrics immediately yield arbitrarily high algebraic convergence rates, i.e., for every and , we havewhenever and represent the native space of RBFs. A considerable amount of work has gone into investigating the dependence of the constant on [43].
In this work, MQ-RBF is used, so it is concluded thatwhere and are the exact and approximate solution, respectively. Now let us assume that the scheme (21) is order accurate in space, then
Let us define the residual by , then
By Lax-Richtmyer definition of stability, the scheme (21) is stable ifwhen matrix is normal then ; otherwise, the inequality is always true. It is assumed that the step size is to be small enough, and the solution and IC of the given problem must be sufficiently smooth. We must have to keep constant. Therefore, there exist some constant such that
Since the residual obeys zero IC and BCs, so . So by mathematical induction, we haveUsing the condition given in Eq. (28), we have
Hence, the scheme is convergent.
3. Formulation of the Method for Case II
In this section, we take and find by using Caputo derivative. We approximate the 2nd order time derivative (appearing in the Caputo derivative) by the central difference scheme. We then apply -weighted scheme to the governing equation.
3.1. Time-Fractional Derivative
Caputo fractional derivative for is defined by
Taking derivative at , we getThis implies
Now, we approximate by finite difference scheme, as follows:
Repeating the same process as we did in Section 2.1, we getwhere and ,.
Now there exist for and , for which we use the second IC, i.e.,
Hence, we get the following value of ,
3.2. -Weighted Scheme
Applying -weighted scheme to Eq. (3) and substituting Eq. (38) in Eq. (3), we get the following scheme for and , respectively.Proceeding in the same way as in the Section 2.2, we get and , respectively,where
Also,where
We can use Eq. (41) and Eq. (43) to find the solution at any time level for .
4. Numerical Results
This section is devoted to the numerical implementation of the schemes constructed in Section 2 and Section 3. We have applied schemes over six problems including one-dimensional and two-dimensional time-fractional PDEs. Problems with different types of domains and geometries are also included. We assess the accuracy of the method by taking different values of and . We have utilized the following error norms:
Computational orders for time and space are calculated, respectively, by the formula,
Figures are incorporated to show the performance of the method. We have applied IS and CNS and have compared the results. An attempt is made to apply the schemes for some nonuniform nodes including Chebyshev, random, Halton, and scattered data nodes. Also, numerical simulation is performed for some irregular domains. Moreover, we have compared our results with the results reported in [28, 44]. Convergence order is calculated in all the problems, and there is uniform convergence in all the problems.wherewith the exact solution
Problem 3. In the first problem, we take the following time-fractional PDE [23].
IC is given as with homogeneous BCs. In Figure 1, is plotted for different values of . It can be observed that increasing value of causes less accuracy of results. Also, it is clear that IS gives more accurate results than CNS. Figure 2 displays exact/approximate solution and absolute error at and . In Figure 3, numerical solutions at different time levels are plotted. Table 1 is concerned with the convergence order, and it shows that the scheme is convergent as proved theoretically. In Table 2, we have computed the numerical results utilizing the IS and the CNS for various values of s. It is observed from the table that the IS produced better results as compared to CNS.wherewith the exact solutionwhere IC is considered as while the BCs are , . In Table 3, IS and CNS are compared with collocation finite element scheme (CFES) [44] which indicates the admirable performance of both the schemes. Figure 4 shows the behaviour of IS scheme by plotting absolute error and approximate/exact solution. One can examine that error decays with increasing . In Figure 5, we have plotted numerical results for different values of . Convergence order is calculated in Table 4, and the uniform convergence is obtained.wherewith the exact solution
Problem 4. We consider the following problem [44].
Problem 5. In this problem, we take the following time-fractional PDE [44]
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
IC is and with BCs . Table 5 shows the performance of IS, CNS, and CFES. One can see that IS and CNS give less error than CFES. Figure 6 is aimed at showing absolute error and exact and approximate solution for . We can inspect that results obtained from IS scheme agree with that obtained from the exact solution. Figure 7 concerns with the numerical solution at different time levels. In Table 6, the order of convergence is calculated which shows the convergence of scheme.wherewith the exact solution
Problem 6. In this problem, we take [28] and consider
(a)
(b)
(a)
(b)
ICs are given by and with BCs . In Table 7, error norm is given for IS, CNS, and meshless Galerkin method (MGM) [28]. We can see that IS produces more accurate results than both CNS and MGM. Figure 8 displays and relationship between exact and numerical solution. Figure 9 is devoted to plot solution at different time levels. In Tables 8 and 9, convergence orders are calculated in time and space, respectively. It can be seen that the method is uniformly convergent in time as well as space.whereand the exact solution iswith the ICsand the BCs
Problem 7. In this problem, two-dimensional time fractional diffusion-wave equation [(40)] is considered
(a)
(b)
(a)
(b)
In this problem, we have taken , , and . We have plotted approximate solution for uniform nodes and different values of in Figure 10. Numerical simulation is performed for irregular domains as well. Numerical solution for different types of irregular domains are presented in Figures 11–14. It can be seen that the better approximate solution is achieved for the irregular domain.with the exact solutionwhere Mittag-Leffler function (one parameter) is defined byWe consider IC from and the BCs are drawn from the exact solution. We have performed computations for , , , and In Figure 15, we have given an approximate solution for different values of for uniform nodes, while in Figure 16, we have plotted error norm for different types of nonuniform nodes for . We got reasonable accuracy for these cases as well.
Problem 8. In this problem, we take the following time FPDE [45].
(a)
(b)
(c)
(d)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(d)
(a) Chebyshev nodes
(b) using Chebyshev nodes
(c) Random nodes
(d) using random nodes
(e) Halton nodes
(f) using Halton nodes
(g) Scattered data nodes
(h) using scattered data nodes
5. Conclusion
In this work, an attempt is made to propose IS and CNS schemes for the solution of time-fractional PDEs. The time derivative is defined and simplified in Caputo sense and then its value is substituted in governing equation along with replacement of space derivatives by RBFs. This paper has an edge over the existing papers in the sense that in this paper, the scheme is constructed for and and for IS and CNS. Problems are given to show the behaviour of the method. Numerical results for different values of are demonstrated to examine the effect of over solution. Results produced by IS are compared with that by CNS. Results are also compared with some other methods in the literature. This comparison clearly indicates the impressive performance of our schemes. In order to utilize the advantage of RBF collocation method for nonuniform nodes and irregular domain, numerical simulation is performed and remarkable results are achieved for nonuniform nodes and irregular domain.
Data Availability
Data will be available on request.
Conflicts of Interest
The authors declare that there are no conflicts of interest associated with this publication.
Acknowledgments
The research was supported by the National Natural Science Foundation of China (Grant No. 11501170). The first author is grateful to Shaheed Benazir Bhutto Women University Peshawar, Pakistan, for study leave.