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

Authors: Mehnaz Shakeel; Iltaf Hussain; Hijaz Ahmad; Imtiaz Ahmad; Phatiphat Thounthong; Ying-Fang Zhang
Journal: Journal of Function Spaces
 (2020)

---

## 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α∈0,1 and α∈1,2. Forward difference scheme is applied to approximate the 1st order derivative appearing in the definition of Caputo derivative for α∈0,1, whereas central difference scheme is used for the 2nd order derivative in the definition of Caputo derivative for α∈1,2. 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.


---

## Body

## 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 by(1)φrij=rij2+c2,where rij=∥zi−zj∥, i,j=1,2,⋯,N, N is the number of collocation points and c is the shape parameter. Furthermore, z=z1 in one dimensional, and z=z1,z2 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(2)∂αv∂tα=1Γ1−α∫0t∂η+1v∂tη+1t−τη−αdτ,η<α<η+1,∂η+1v∂tη+1,α=η+1.Definition 1.
Caputo derivative of noninteger orderαof a functionvz,tis defined by [9]In this paper, we have tackled the following two cases of Caputo derivative:Case I.
0<α<1Case II.
1<α<2


### 1.1. Governing Equation

We emphasize on the following time-fractional PDE(3)∂αv∂tα+Lv=ψz,t,z∈Ω,η<α<η+1,t>0.The boundary conditions (BCs) are(4)vz,t=Gz,t,z∈∂Ω,t>0,and the initial conditions (ICs) are(5)vz,0=μz,vtz,0=νz,where vz,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 Lv. Here, Lv is a linear operator defined by(6)Lv=aΔv+b∇v+cv,where a, b, and c are functions of z 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: Section2 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.


## 1.1. Governing Equation

We emphasize on the following time-fractional PDE(3)∂αv∂tα+Lv=ψz,t,z∈Ω,η<α<η+1,t>0.The boundary conditions (BCs) are(4)vz,t=Gz,t,z∈∂Ω,t>0,and the initial conditions (ICs) are(5)vz,0=μz,vtz,0=νz,where vz,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 Lv. Here, Lv is a linear operator defined by(6)Lv=aΔv+b∇v+cv,where a, b, and c are functions of z 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: Section2 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.


## 2. Formulation of the Method for CaseI

In this section, we take0<α<1 and find ∂αv/∂tα, 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α∈0,1 is defined by(7)∂αv∂tα=1Γ1−α∫0t∂v∂tt−τ−αdτ,0<α<1,∂v∂t,α=1.Taking derivative att=tn+1, we get(8)∂αvz,tn+1∂tα=1Γ1−α∫0tn+1∂v∂ttn+1−τ−αdτ.This implies(9)∂αv∂tα=1Γ1−α∑k=0n∫tktk+1∂v∂ttn+1−τ−αdτ.Now, we use finite difference scheme to approximate∂n+1v/∂tn+1, as follows:(10)∂v∂t=vk+1−vkdt+odt,where dt is the time step size.Then(11)∂αvz,tn+1∂tα≈1Γ1−α∑k=0n∫tktk+1vk+1−vkdt+odttn+1−τ−αdτ,∂αv∂tα=aα∑k=0nbkvn−k+1−vn−k+odt2−α,where aα=dt−α/Γ2−α and bk=k+11−α−k1−α, k=0,1,⋯,n.Finally, we can write in more precise form as(12)∂αv∂tα=aαvn+1−vn+aα∑k=1nbkvn−k+1−vn−k,n≥1,aαv1−v0,n=0.


### 2.2.θ-Weighted Scheme

Applyingθ-weighted scheme to Eq. (3) and substituting Eq. (12) in Eq. (3), we get(13)aαvn+1+θLvn+1=aαvn−1−θLvn−aα∑k=1nbkvn−k+1−vn−k+ψn+1,n≥1,aαv0−1−θLv0+ψ1,n=0.Now, we employ Kansa’s method and interpolatevz,tn+1 by RBFs. We interpolate the solution at N different points zi∣i=1,2,⋯,N, where zi∣i∈Ω are interior points while z1 and zN are boundary points.(14)vz,tn+1=∑i=1Nλin+1φ∥zi−zj∥,j=1,2,⋯,N=∑i=1Nλin+1φrij,where φrij are the RBFs, . is Euclidian norm, and λi’s are the unknown constants.We can write Eq. (14) in matrix form as(15)vn+1=Aλn+1or(16)λn+1=A−1vn+1.Provided that collocation matrixA 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 λin+1, we can find solution v from Eq. (14).Putting the value from Eq. (15) in Eq. (13), we get the following form(17)aαAλn+1+θLAλn+1=aαAλn−1−θLAλn−aα∑k=1nbkvn−k+1−vn−k+ψn+1,n≥1,aαAλ0−1−θLAλ0+ψ1,n=0.(18)λn+1=M−1ℕλn+M−1Gn+1,n≥0,where(19)M=aαA+θLA,N=aαA−1−θLA,MGn+1=G1n+1+G2n+1,n≥1,ψ1,n=0.Here,(20)G1n+1=g1n+1,0,⋯,g2n+1,G2n+1=−aα∑k=1nbkvn−k+1−vn−k+ψn+1,where g1n+1 and g2n+1 are some known functions given in BCs. Finally, from Eq. (18) and Eq. (16), we get(21)vn+1=AM−1NA−1vn+AM−1Gn+1.We can use this scheme to find the solution at any time leveltn.


### 2.3. Stability and Convergence

The scheme (21) is basically a recurrence relation which gives the values of solution at the time tn+1 by using solution at time tn. The matrix E=AM−1NA−1 is known as amplification matrix, and the elements of 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 vn at time tn.We state a well-known theorem of Fasshauer [41] (see [42] for proof):(22)Dαfz−DαPfz≤Chχ,Ωk−∣α∣CΦxfNϕΩ,provided that hχ,Ω≤h0. Here,(23)CΦz=maxβ,γ∈ℕ0s∣β∣+∣γ∣=2kmaxw,z∈Ω∩Bx,c2hχ,ΩD1βD2γΦw,z.Theorem 2.
[41] Suppose Ω⊆ℝs is open and bounded and satisfies an interior cone condition. Suppose Φ∈C2kΩ×Ω is symmetric and strictly conditionally positive definite of order m on ℝs. Denote the interpolant to f∈NϕΩ on the m−1-unisolvent set χbyPf. Fixα∈ℕ0s with ∣α∣≤k. Then, there exist positive constants h0 and C (independent of z, f, and Φ) such thatApplication of Theorem2 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(24)Dαfz−DαPfz≤Ckhk−∣α∣vNϕΩ,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 Ck on k[43].In this work, MQ-RBF is used, so it is concluded that(25)Dαv^z−Dαvz≤Ckhk−∣α∣vNϕΩ,where v^ and v are the exact and approximate solution, respectively. Now let us assume that the scheme (21) is pth order accurate in space, then(26)vn+1=Evn+AM−1Gn+1+odt2−α+hp,dt,h→0.Let us define the residual byεn=v∧n−vn, then(27)εn+1=Eεn+odt2−α+hp,dt,h→0.By Lax-Richtmyer definition of stability, the scheme (21) is stable if(28)E≤1,when matrix E is normal then E=ρE; otherwise, the inequality ρE≤E 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/hp constant. Therefore, there exist some constant C such that(29)εn+1≤Eεn+Cdt2−α+hp,n=0,1,2,⋯,T×M.Since the residualεn obeys zero IC and BCs, so ε0=0. So by mathematical induction, we have(30)εn+1≤1+E2+E3+⋯+En−1Cdt2−α+hp,n=0,1,2,⋯,T×M.Using the condition given in Eq. (28), we have(31)εn+1≤nCdt2−α+hp,n=0,1,2,⋯,T×M.Hence, the scheme is convergent.


## 2.1. Time-Fractional Derivative

Caputo fractional derivative forα∈0,1 is defined by(7)∂αv∂tα=1Γ1−α∫0t∂v∂tt−τ−αdτ,0<α<1,∂v∂t,α=1.Taking derivative att=tn+1, we get(8)∂αvz,tn+1∂tα=1Γ1−α∫0tn+1∂v∂ttn+1−τ−αdτ.This implies(9)∂αv∂tα=1Γ1−α∑k=0n∫tktk+1∂v∂ttn+1−τ−αdτ.Now, we use finite difference scheme to approximate∂n+1v/∂tn+1, as follows:(10)∂v∂t=vk+1−vkdt+odt,where dt is the time step size.Then(11)∂αvz,tn+1∂tα≈1Γ1−α∑k=0n∫tktk+1vk+1−vkdt+odttn+1−τ−αdτ,∂αv∂tα=aα∑k=0nbkvn−k+1−vn−k+odt2−α,where aα=dt−α/Γ2−α and bk=k+11−α−k1−α, k=0,1,⋯,n.Finally, we can write in more precise form as(12)∂αv∂tα=aαvn+1−vn+aα∑k=1nbkvn−k+1−vn−k,n≥1,aαv1−v0,n=0.


## 2.2.θ-Weighted Scheme

Applyingθ-weighted scheme to Eq. (3) and substituting Eq. (12) in Eq. (3), we get(13)aαvn+1+θLvn+1=aαvn−1−θLvn−aα∑k=1nbkvn−k+1−vn−k+ψn+1,n≥1,aαv0−1−θLv0+ψ1,n=0.Now, we employ Kansa’s method and interpolatevz,tn+1 by RBFs. We interpolate the solution at N different points zi∣i=1,2,⋯,N, where zi∣i∈Ω are interior points while z1 and zN are boundary points.(14)vz,tn+1=∑i=1Nλin+1φ∥zi−zj∥,j=1,2,⋯,N=∑i=1Nλin+1φrij,where φrij are the RBFs, . is Euclidian norm, and λi’s are the unknown constants.We can write Eq. (14) in matrix form as(15)vn+1=Aλn+1or(16)λn+1=A−1vn+1.Provided that collocation matrixA 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 λin+1, we can find solution v from Eq. (14).Putting the value from Eq. (15) in Eq. (13), we get the following form(17)aαAλn+1+θLAλn+1=aαAλn−1−θLAλn−aα∑k=1nbkvn−k+1−vn−k+ψn+1,n≥1,aαAλ0−1−θLAλ0+ψ1,n=0.(18)λn+1=M−1ℕλn+M−1Gn+1,n≥0,where(19)M=aαA+θLA,N=aαA−1−θLA,MGn+1=G1n+1+G2n+1,n≥1,ψ1,n=0.Here,(20)G1n+1=g1n+1,0,⋯,g2n+1,G2n+1=−aα∑k=1nbkvn−k+1−vn−k+ψn+1,where g1n+1 and g2n+1 are some known functions given in BCs. Finally, from Eq. (18) and Eq. (16), we get(21)vn+1=AM−1NA−1vn+AM−1Gn+1.We can use this scheme to find the solution at any time leveltn.


## 2.3. Stability and Convergence

The scheme (21) is basically a recurrence relation which gives the values of solution at the time tn+1 by using solution at time tn. The matrix E=AM−1NA−1 is known as amplification matrix, and the elements of 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 vn at time tn.We state a well-known theorem of Fasshauer [41] (see [42] for proof):(22)Dαfz−DαPfz≤Chχ,Ωk−∣α∣CΦxfNϕΩ,provided that hχ,Ω≤h0. Here,(23)CΦz=maxβ,γ∈ℕ0s∣β∣+∣γ∣=2kmaxw,z∈Ω∩Bx,c2hχ,ΩD1βD2γΦw,z.Theorem 2.
[41] Suppose Ω⊆ℝs is open and bounded and satisfies an interior cone condition. Suppose Φ∈C2kΩ×Ω is symmetric and strictly conditionally positive definite of order m on ℝs. Denote the interpolant to f∈NϕΩ on the m−1-unisolvent set χbyPf. Fixα∈ℕ0s with ∣α∣≤k. Then, there exist positive constants h0 and C (independent of z, f, and Φ) such thatApplication of Theorem2 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(24)Dαfz−DαPfz≤Ckhk−∣α∣vNϕΩ,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 Ck on k[43].In this work, MQ-RBF is used, so it is concluded that(25)Dαv^z−Dαvz≤Ckhk−∣α∣vNϕΩ,where v^ and v are the exact and approximate solution, respectively. Now let us assume that the scheme (21) is pth order accurate in space, then(26)vn+1=Evn+AM−1Gn+1+odt2−α+hp,dt,h→0.Let us define the residual byεn=v∧n−vn, then(27)εn+1=Eεn+odt2−α+hp,dt,h→0.By Lax-Richtmyer definition of stability, the scheme (21) is stable if(28)E≤1,when matrix E is normal then E=ρE; otherwise, the inequality ρE≤E 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/hp constant. Therefore, there exist some constant C such that(29)εn+1≤Eεn+Cdt2−α+hp,n=0,1,2,⋯,T×M.Since the residualεn obeys zero IC and BCs, so ε0=0. So by mathematical induction, we have(30)εn+1≤1+E2+E3+⋯+En−1Cdt2−α+hp,n=0,1,2,⋯,T×M.Using the condition given in Eq. (28), we have(31)εn+1≤nCdt2−α+hp,n=0,1,2,⋯,T×M.Hence, the scheme is convergent.


## 3. Formulation of the Method for CaseII

In this section, we take1<α<2 and find ∂αv/∂tα 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α∈1,2 is defined by(32)∂αv∂tα=1Γ2−α∫0t∂2v∂t2t−τ1−αdτ,1<α<2,∂2v∂t2,α=2.Taking derivative att=tn+1, we get(33)∂αvz,tn+1∂tα=1Γ2−α∫0tn+1∂2v∂t2tn+1−τ1−αdτ.This implies(34)∂αv∂tα=1Γ2−α∑k=0n∫tktk+1∂2v∂t2tn+1−τ1−αdτ.Now, we approximate∂2v/∂t2 by finite difference scheme, as follows:(35)∂2vz,τk+1∂τ2=vk+1−2vk+vk−1dt2+odt2.Repeating the same process as we did in Section2.1, we get(36)∂αv∂tα=aαvn+1−2vn+vn−1+aα∑k=1nbkvn−k+1−2vn−k+vn−k−1+odt3−α,n≥0,where aα=dt−α/Γ3−α and bk=k+12−α−k2−α,k=0,1,⋯,n.Now there existv−1 for n=0 and k=n, for which we use the second IC, i.e.,(37)vt0=v1−v−12dt,v−1=v1−2dtvt0.Hence, we get the following value of∂αv/∂tα,(38)∂αv∂tα=aαvn+1−2vn+vn−1+aα∑k=1n−1bkvn−k+1−2vn−k+vn−k−1+2aαbnv1−v0−dtvt0,n≥1,2aαv1−v0−dtvt0,n=0.


### 3.2.θ-Weighted Scheme

Applyingθ-weighted scheme to Eq. (3) and substituting Eq. (38) in Eq. (3), we get the following scheme for n=0 and n≥1, respectively.(39)2aαv1+θLv1=2aαv0+2aαdtvt0−1−θLv0+ψ1,n=0(40)aαvn+1+θLvn+1=2aαvn−aαvn−1−2aαbnv1−v0−dtvt0−aα∑k=1n−1bkvn−k+1−2vn−k+vn−k−1−1−θLvn+ψn+1,n≥1.Proceeding in the same way as in the Section 2.2, we get v1 and vn+1, respectively,(41)v1=AM−1ℕA−1v0+AM−1ℍ,where(42)M=2aαA+θLA,ℕ=2aαA−1−θLA,ℍ=ℍ1+ℍ2,ℍ1=g1,0,⋯,g2,ℍ2=ψ1+2aαdtvt0.Also,(43)vn+1=AM1−1ℕA−1vn+AM1−1ℙA−1vn−1+AM1−1Gn+1,n≥1,where(44)M1=aαA+θLA,ℙ=−aαAλn−1,Gn+1=G1n+1+G2n+1,G1n+1=g1n+1,0,⋯,g2n+1,G2n+1=−2aαbnv1−v0−dtvt0−aα∑k=1n−1bkvn−k+1−2vn−k+vn−k−1+ψn+1.We can use Eq. (41) and Eq. (43) to find the solution at any time level tn for n≥1.


## 3.1. Time-Fractional Derivative

Caputo fractional derivative forα∈1,2 is defined by(32)∂αv∂tα=1Γ2−α∫0t∂2v∂t2t−τ1−αdτ,1<α<2,∂2v∂t2,α=2.Taking derivative att=tn+1, we get(33)∂αvz,tn+1∂tα=1Γ2−α∫0tn+1∂2v∂t2tn+1−τ1−αdτ.This implies(34)∂αv∂tα=1Γ2−α∑k=0n∫tktk+1∂2v∂t2tn+1−τ1−αdτ.Now, we approximate∂2v/∂t2 by finite difference scheme, as follows:(35)∂2vz,τk+1∂τ2=vk+1−2vk+vk−1dt2+odt2.Repeating the same process as we did in Section2.1, we get(36)∂αv∂tα=aαvn+1−2vn+vn−1+aα∑k=1nbkvn−k+1−2vn−k+vn−k−1+odt3−α,n≥0,where aα=dt−α/Γ3−α and bk=k+12−α−k2−α,k=0,1,⋯,n.Now there existv−1 for n=0 and k=n, for which we use the second IC, i.e.,(37)vt0=v1−v−12dt,v−1=v1−2dtvt0.Hence, we get the following value of∂αv/∂tα,(38)∂αv∂tα=aαvn+1−2vn+vn−1+aα∑k=1n−1bkvn−k+1−2vn−k+vn−k−1+2aαbnv1−v0−dtvt0,n≥1,2aαv1−v0−dtvt0,n=0.


## 3.2.θ-Weighted Scheme

Applyingθ-weighted scheme to Eq. (3) and substituting Eq. (38) in Eq. (3), we get the following scheme for n=0 and n≥1, respectively.(39)2aαv1+θLv1=2aαv0+2aαdtvt0−1−θLv0+ψ1,n=0(40)aαvn+1+θLvn+1=2aαvn−aαvn−1−2aαbnv1−v0−dtvt0−aα∑k=1n−1bkvn−k+1−2vn−k+vn−k−1−1−θLvn+ψn+1,n≥1.Proceeding in the same way as in the Section 2.2, we get v1 and vn+1, respectively,(41)v1=AM−1ℕA−1v0+AM−1ℍ,where(42)M=2aαA+θLA,ℕ=2aαA−1−θLA,ℍ=ℍ1+ℍ2,ℍ1=g1,0,⋯,g2,ℍ2=ψ1+2aαdtvt0.Also,(43)vn+1=AM1−1ℕA−1vn+AM1−1ℙA−1vn−1+AM1−1Gn+1,n≥1,where(44)M1=aαA+θLA,ℙ=−aαAλn−1,Gn+1=G1n+1+G2n+1,G1n+1=g1n+1,0,⋯,g2n+1,G2n+1=−2aαbnv1−v0−dtvt0−aα∑k=1n−1bkvn−k+1−2vn−k+vn−k−1+ψn+1.We can use Eq. (41) and Eq. (43) to find the solution at any time level tn for n≥1.


## 4. Numerical Results

This section is devoted to the numerical implementation of the schemes constructed in Section2 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 t and α. We have utilized the following error norms:(45)Eabs=v^i−vi,i=1,2,⋯,N,E∞=maxv^i−vi,ERMS=1N∑i=1NEabs2.Computational orders for time and space are calculated, respectively, by the formula,(46)C1−order=logE∞dt1/E∞dt2logdt1/dt2,C2−order=logE∞dx1/E∞dx2logdx1/dx2.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.(47)∂αv∂tα−∂2v∂z12=ψz1,t,0≤z1≤1,0<α<1,t>0,where(48)ψz1,t=2Γ3−αt2−αsin2πz1+4π2t2sin2πz1,with the exact solution(49)vz1,t=t2sin2πz1.Problem 3.
In the first problem, we take the following time-fractionalPDE [23].IC is given asvz1,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, 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.(50)∂αv∂tα−∂2v∂z12+v=ψz1,t,0≤z1≤1,0<α<1,t>0,where(51)ψz1,t=Γ2+αez1t,with the exact solution(52)vz1,t=ez1t1+α,where IC is considered as vz1,0=0 while the BCs are v0,t=t1+α, v1,t=et1+α. 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.(53)∂αv∂tα−∂2v∂z12−v=ψz1,t,0≤z1≤π,0<α<1,t>0,where(54)ψz1,t=2Γ3−αt2−αsinz1,with the exact solution(55)vz1,t=t2sinz1.Problem 4.
We consider the following problem [44].Problem 5.
In this problem, we take the following time-fractional PDE [44]Figure 1
Comparison ofE∞ of IS and CNS for the Problem 3.Figure 2
(a)Eabs and (b) exact and approximate solution at α=0.2 for the Problem 3.
(a)(b)Figure 3
Approximate solution at different time levels by using (a) IS and (b) CNS for the Problem3.
(a)(b)Table 1
Error norms and computational order for different values ofdt with c=5.1,N=50 for the Problem 3.
θdtE∞ERMSC1-order10.13.2163e-051.7837e-05—0.052.6247e-051.0503e-050.29320.011.4223e-051.1387e-050.38070.0056.7040e-063.0043e-061.08510.0014.0082e-061.7611e-060.3196120.11.0700e-017.4900e-02—0.055.1100e-023.5800e-021.06620.019.8000e-036.9000e-031.02610.0054.9000e-033.4000e-031.00000.0019.8295e-046.8233e-040.9981Table 2
Error norms for different values ofα for the Problem 3.
αISCNSE∞ERMSE∞ERMS0.27.7846e-054.7558e-051.0700e-017.4900e-020.42.2439e-041.4815e-041.0710e-017.5000e-020.65.5165e-043.8107e-041.0750e-017.5300e-020.81.2000e-038.3507e-041.0830e-017.5900e-02Table 3
Comparison with ofE∞ [35] for the Problem 4.
αISCNS[35]0.202.0316e-063.9338e-062.3228e-050.509.9440e-073.6707e-061.9343e-050.752.3590e-061.6847e-061.4279e-050.907.3325e-067.9323e-068.4900e-06Figure 4
(a)Eabs and (b) exact and approximate solution at α=0.2 for the Problem 4.
(a)(b)Figure 5
Approximate solution at different time levels by using (a) IS and (b) CNS for the Problem4.
(a)(b)Table 4
Error norms and computational order for different time steps withc=5.1,N=50 and t=1 for the Problem 4.
θdtE∞ERMSC1−order10.11.2261e-049.3157e-05–0.051.9822e-051.5126e-052.62890.016.4769e-062.8484e-060.69500.0054.0646e-062.2878e-060.67220.0012.0397e-068.3217e-070.4284120.11.1200e-028.5000e-03–0.055.1000e-033.9000e-01.13500.011.1000e-037.6660e-040.95270.0055.3126e-043.8540e-041.05000.0011.0499e-047.7472e-050.9756IC isvz1,0=0 and with BCs v0,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 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.(56)∂αv∂tα−∂2v∂z12=ψz1,t,0≤z1≤1,1<α<2,t>0,where(57)ψz1,t=2t2−αΓ3−αz1z1−1−2t2,with the exact solution(58)vz1,t=t2z1z1−1.Problem 6.
In this problem, we takeLv=−∂2v/∂z12[28] and considerTable 5
Comparison with ofE∞ [44] for the Problem 3.
αISCNS[44]0.203.3143e-075.5885e-062.6640e-050.502.9331e-072.5374e-068.1200e-060.751.0587e-062.9157e-069.5300e-060.905.4445e-063.9059e-064.3240e-05Figure 6
(a)Eabs and (b) exact and approximate solution for α=0.50 for the Problem 5.
(a)(b)Figure 7
Approximate solution at different time levels by using (a) IS and (b) CNS for the Problem3.
(a)(b)Table 6
Error norms and computational order for different time steps forN=50,c=8.1 and t=1 for the Problem 3.
θdtE∞ERMSC1-order10.11.4900e-021.0400e-02—0.055.5000e-033.8000e-033.31070.015.0083e-043.4937e-043.42820.0051.8363e-041.2741e-043.33300.0011.4575e-059.1709e-063.6248120.19.0800e-026.3600e-02—0.054.3100e-023.0200e-022.47530.018.0000e-035.6000e-032.40940.0053.9000e-032.7000e-032.38670.0018.3394e-045.7298e-042.2069ICs are given byvz1,0=0,vtz1,0=0 and with BCs v0,t=0,v1,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 Eabs 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.(59)∂αv∂tα=∂2v∂z12+∂2v∂z22+ψz1,z2,t,Ω=0,1×0,1,1<α<2,t∈0,1,where(60)ψz1,z2,t=sinπz1sinπz2Γ3+α2t2+2t2+α,and the exact solution is(61)v=t2+αsinπz1sinπz2,with the ICs(62)vz1,z2,t=0,vtz1,z2,t=0,z1,z2∈Ω,t=0and the BCs(63)vz1,z2,t=t2+αsinπz1sinπz2,z1,z2∈∂Ω.Problem 7.
In this problem, two-dimensional time fractional diffusion-wave equation [(40)] is consideredTable 7
Comparison ofE∞ with [20] for the Problem 4.
αISCNS[20]1.154.6987e-049.0819e-049.1195e-041.254.7856e-049.3809e-049.0865e-041.54.9736e-049.7102e-048.9931e-04Figure 8
(a)Eabs (b) exact and approximate solution for α=1.15 for the Problem 6.
(a)(b)Figure 9
Approximate solution at different time levels by using (a) IS and (b) CNS for Problem4.
(a)(b)Table 8
Error norms att=1 for different time steps with c=5.9,N=50 for the Problem 4.
θdtE∞ERMSC1−order10.15.2500e-023.8000e-02—0.052.5600e-021.8500e-021.03620.015.0000e-033.6000e-031.01470.0052.5000e-031.8000e-031.00000.0011.2000e-038.9473e-040.4560120.17.7300e-025.5900e-02—0.053.7400e-022.7100e-021.04740.017.3000e-035.3000e-031.01510.0053.7000e-032.7000e-030.98040.0019.6490e-047.0544e-040.8351Table 9
Error norms and computational order for differentdx with c=100,α=1.5, and t=0.1 for the Problem 4.
θdxE∞ERMSC2-order10.12.1000e-031.5000e-03—0.051.7000e-031.2000e-030.30490.014.0345e-043.2477e-040.89370.0051.3119e-041.1574e-041.62070.0014.8815e-053.5377e-050.6143120.12.1000e-031.5000e-03—0.051.7000e-031.2000e-030.30490.013.9522e-043.1891e-040.90650.0051.2686e-041.1204e-041.63940.0015.0791e-053.7090e-050.5687In this problem, we have takenc=5.5, dt=0.001, T=0.1, and N=30. 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.(64)∂αv∂tα=∂2v∂z12+∂2v∂z22,z1,z2∈Ω=0,1×0,1,0<α<1,t∈0,0.5,with the exact solution(65)vz1,z2,t=Eα−12π2tαcosπ2z1cosπ2z2,where Mittag-Leffler function (one parameter) is defined by(66)Eαξ=∑j=0∞ξjΓαj+1,α>0.We consider IC from vz1,z2,0=cosπ/2z1cosπ/2z2, 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].Figure 10
Approximate solution for different values ofα for the Problem 7.
(a)
α=1.2(b)
α=1.5(c)
α=1.7(d)
α=1.9Figure 11
(a) Irregular domain and (b) absolute error for Problem7.
(a)(b)Figure 12
(a) Irregular domain and (b) numerical solution for the Problem7.
(a)(b)Figure 13
(a) Irregular domain and (b) numerical solution for the Problem7.
(a)(b)Figure 14
(a) Irregular domain and (b) absolute error for Problem7.
(a)(b)Figure 15
Approximate solution for different values ofα for the Problem 8.
(a)
α=0.1(b)
α=0.3(c)
α=0.6(d)
α=0.9Figure 16
E∞ with different nonuniform nodes for the Problem 8.
(a)
Chebyshev nodes(b)
E∞ using Chebyshev nodes(c)
Random nodes(d)
E∞ using random nodes(e)
Halton nodes(f)
E∞ using Halton nodes(g)
Scattered data nodes(h)
E∞ 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 for0<α<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.



---

*Source: 8898309-2020-10-24.xml*