Meshless Technique for the Solution of Time-Fractional Partial Differential Equations Having Real-World Applications

Department of Basic Sciences, University of Engineering and Technology, Peshawar, Pakistan Department of Mathematics, Shaheed Benazir Bhutto Women University, Peshawar, Pakistan Section of Mathematics, International Telematic University Uninettuno, Corso Vittorio Emanuele II, 39, 00186 Roma, Italy Department of Mathematics, University of Swabi, Khyber Pakhtunkhwa, Pakistan Renewable Energy Research Centre, Department of Teacher Training in Electrical Engineering, Faculty of Technical Education, King Mongkut’s University of Technology North Bangkok, 1518 Pracharat 1 Road, Bangsue, Bangkok 10800, Thailand School of Mathematics and Information Science, Henan Polytechnic University, Jiaozuo 454000, China


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][2][3][4][5][6][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][9][10][11]. Many researchers did considerable work to find the analytic solution of FPDEs [12][13][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][16][17][18][19][20]. There are two widely used def-initions 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][26][27][28][29][30][31][32][33][34][35]. In this paper, multiquadric (MQ) RBF is used to approximate solution. MQ-RBF is defined by where r ij = ∥z i − z j ∥, i, j = 1, 2, ⋯, N, N is the number of collocation points and c is the shape parameter. Furthermore, z = z 1 in one dimensional, and z = ðz 1 , z 2 Þ in twodimensional 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][38][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 order α of a function vðz, tÞ is defined by [9] In this paper, we have tackled the following two cases of Caputo derivative: We emphasize on the following time-fractional PDE The boundary conditions (BCs) are where vðz, tÞ is the solution, ∂ α v/∂t α is the Caputo fractional derivative of order α, ψðz, tÞ 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 LðvÞ. Here, LðvÞ is a linear operator defined by where a, b, and c are functions of z or constants and Δ and ∇ denote Laplacian and gradient operator, respectively.

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., 0 < α < 1. In Section 3, we consider the second case, i.e., 1 < α < 2. 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.

Formulation of the Method for Case I
In this section, we take 0 < α < 1 and find ∂ α v/∂t α , by using Caputo derivative. Finite difference scheme is applied to approximate the 1 st 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 α ∈ ð0, 1Þ is defined by Taking derivative at t = t n+1 , we get This implies Now, we use finite difference scheme to approximate ∂ n+1 v/∂t n+1 , as follows: Journal of Function Spaces where dt is the time step size. Then where a α = dt −α /Γð2 − αÞ and b k = ðk + 1Þ 1−α − k 1−α , k = 0, 1, ⋯, n. Finally, we can write in more precise form as where φðr ij Þ are the RBFs, k:k is Euclidian norm, and λ i ' s are the unknown constants.
We can write Eq. (14) in matrix form as v n+1 = Aλ n+1 ð15Þ Provided that collocation matrix A must be nonsingular to ensure the invertibility of matrix A. This depends on the choice of RBF and the location of mesh points. Matrix A is invertible for distinct mesh points. The shape parameter has an important effect on condition number [40]. Once we find the constants λ n+1 i , we can find solution v from Eq. (14). Putting the value from Eq. (15) in Eq. (13), we get the following form Here, where g n+1 1 and g n+1 2 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 t n .    Journal of Function Spaces this matrix depend on the constant κ = dt/h ς , where dt is the time step, h 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 v n at time t n . We state a well-known theorem of Fasshauer [41] (see [42] for proof): provided that h χ,Ω ≤ h 0 . Here, Theorem 2. [41] Suppose Ω ⊆ ℝ s is open and bounded and satisfies an interior cone condition. Suppose Φ ∈ C 2k ðΩ × ΩÞ is symmetric and strictly conditionally positive definite of order m on ℝ s . Denote the interpolant to f ∈ N ϕ ðΩÞ on the ðm − 1Þ -unisolvent set χ by P f . Fix α ∈ ℕ s 0 with |α | ≤k. Then, there exist positive constants h 0 and C (independent of z, f , 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 k ∈ ℕ and |α | ≤k, we have whenever f ∈ N ϕ ðΩÞ and N ϕ ðΩÞ represent the native space of RBFs. A considerable amount of work has gone into investigating the dependence of the constant C k on k [43]. In this work, MQ-RBF is used, so it is concluded that wherev and v are the exact and approximate solution, respectively. Now let us assume that the scheme (21) is p th order accurate in space, then Let us define the residual by ε n = v∧ n − v n , then By Lax-Richtmyer definition of stability, the scheme (21) is stable if when matrix E is normal then kEk = ρðEÞ; otherwise, the inequality ρðEÞ ≤ kEk is always true. It is assumed that the step size h is to be small enough, and the solution and IC of the given problem must be sufficiently smooth. We must have dt → 0 to keep κ = dt/h p constant. Therefore, there exist some constant C such that Since the residual ε n obeys zero IC and BCs, so ε 0 = 0. So by mathematical induction, we have Hence, the scheme is convergent.

Formulation of the Method for Case II
In this section, we take 1 < α < 2 and find ∂ α v/∂t α by using Caputo derivative. We approximate the 2 nd 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 α ∈ ð1, 2Þ is defined by Taking derivative at t = t n+1 , we get    Journal of Function Spaces This implies Now, we approximate ∂ 2 v/∂t 2 by finite difference scheme, as follows: Repeating the same process as we did in Section 2.1, we get where Now there exist v −1 for n = 0 and k = n, for which we use the second IC, i.e., Hence, we get the following value of ∂ α v/∂t α ,
Proceeding in the same way as in the Section 2.2, we get v 1 and v n+1 , respectively, where Also, We can use Eq. (41) and Eq. (43) to find the solution at any time level t n for n ≥ 1.

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 onedimensional 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 t 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. where with the exact solution v z 1 , t ð Þ= t 2 sin 2πz 1 ð Þ: ð49Þ Problem 3. In the first problem, we take the following timefractional PDE [23].
IC is given as vðz 1 , 0Þ = 0 with homogeneous BCs. In Figure 1, E ∞ 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 t = 1, c = 5:1, and dt = 0:001. In Figure 3, numerical solutions at different time levels are plotted. Table 1 is concerned with the convergence order,  Journal of Function Spaces 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.
with the exact solution where IC is considered as vðz 1 , 0Þ = 0 while the BCs are vð0, tÞ = t 1+α , vð1, tÞ = et 1+α . 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 x. In Figure 5, we have plotted numerical results for different values of t. Convergence order is calculated in Table 4, and the uniform convergence is obtained.
where    Journal of Function Spaces with the exact solution v z 1 , t ð Þ= t 2 sin z 1 ð Þ: ð55Þ Problem 4. We consider the following problem [44].
Problem 5. In this problem, we take the following timefractional PDE [44] IC is vðz 1 , 0Þ = 0 and with BCs vð0, tÞ = 0, vðπ, tÞ = 0. 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 α = 0:20. We can inspect that results obtained from IS scheme agree with that obtained from the       Table 6, the order of convergence is calculated which shows the convergence of scheme. where with the exact solution Problem 6. In this problem, we take LðvÞ = −∂ 2 v/∂z 2 1 [28] and consider ICs are given by vðz 1 , 0Þ = 0, v t ðz 1 , 0Þ = 0 and with BCs vð0, tÞ = 0, vð1, tÞ = 0. In Table 7, E ∞ 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 E abs 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, respec-tively. It can be seen that the method is uniformly convergent in time as well as space. where and the exact solution is with the ICs with the exact solution v z 1 , z 2 , t ð Þ= E α − 1 2 π 2 t α cos π 2 z 1 cos where Mittag-Leffler function (one parameter) is defined by We consider IC from vðz 1 , z 2 , 0Þ = cos ðπ/2z 1 Þ cos ðπ/2z 2 Þ, and the BCs are drawn from the exact solution. We have performed computations for c = 13, t = 1, dt = 0:01, and N = 30: In Figure 15, we have given an approximate solution for different values of α for uniform nodes, while in Figure 16, we have plotted E ∞ error norm for different types of nonuniform nodes for N = 100. We got reasonable accuracy for these cases as well.
Problem 8. In this problem, we take the following time FPDE [45].

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 0 < α < 1 and 1 < α < 2 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.