# Numerical Solution of Fractional Order Anomalous Subdiffusion Problems Using Radial Kernels and Transform

Authors: Muhammad Taufiq; Marjan Uddin
Journal: Journal of Mathematics
 (2021)

---

## Abstract

By coupling of radial kernels and localized Laplace transform, a numerical scheme for the approximation of time fractional anomalous subdiffusion problems is presented. The fractional order operators are well suited to handle by Laplace transform and radial kernels are also built for high dimensions. The numerical computations of inverse Laplace transform are carried out by contour integration technique. The computation can be done in parallel and no time sensitivity is involved in approximating the time fractional operator as contrary to finite differences. The proposed numerical scheme is stable and accurate.


---

## Body

## 1. Introduction

In the last decades, many researchers have studied the fractional calculus [1–3]. Differential equations of fractional order have many applications in the field of science and engineering [4–7]. Analytical solution of many fractional differential equations is not possible or very hard to find, so we need a new numerical technique to find its approximate solution. Various phenomena in viscoelastic materials, economics, chemistry, finance, control theory, hydrology, physics, cosmology, solid mechanics, bioengineering, statistical mechanics, and control theory can be mathematically modeled from fractional calculus [8–17]. In literature, various numerical approaches are available for modeling anomalous diffusive behavior such as Carlo simulations [18]. An introduction of diffusion equations can be found in [19–21].Recently, RBF-based methods were used in solving fractional partial differential equations (FPDEs) [22–24]. These methods have been employed in approximation of partial differential equations with complex domains. An implicit meshless technique based on the radial basis functions for the numerical simulation of the anomalous subdiffusion equation can be found in [25]. The convergence and stability of these mesh-free methods can be found in [26, 27]. These globally defined RBF methods cause ill-condition system matrices [28]. To overcome the problem of ill-conditioning, local RBF techniques were used in [29–31]. Unlike global RBF methods, the RBF method in local setting uses center points in each subdomain area of influence, surrounding each spatial point due to which there is reduction in the computational cost.Recently, Laplace transform is combined with RBF method in [32, 33]. In [34–37], the authors use Laplace transform as tool in spectral method and other mesh-based methods such as finite element methods and finite difference method. To avoid the issues of computational efficiency and instability of the system matrix, we introduce a new technique Laplace transform-based local RBF method in solving the time fractional modified anomalous subdiffusion equations in irregular domain.Here, we consider the following modified anomalous subdiffusion equation of fractional order [38]:(1)∂wx,t∂t=ν1Dt1−α+ν2Dt1−βΔwx,t+fx,t,where x=x,y∈Ω⊂ℝd, d≥1, t>0, subject to the following boundary and initial conditions:(2)ℬwx,t=hx,t,x=x,y∈∂Ω,(3)wx,0=w0,x∈Ω,respectively, where α, β∈0,1, t∈0,T, ν1, ν2 are positive constants, Δ is the Laplace operator, and fx,t is some given function.


## 2. Preliminaries

Here, we introduce some fundamental definitions related to fractional calculus [39, 40].Definition 1.
Letn−1<α<n∈Z+ and α>0, then the Caputo derivative of fractional order is defined as(4)Dtαwt=1Γn−α∫0t1t−zα+1−ndndznwzdz.Definition 2.
Letwt, t≥0, be a given function, then its Laplace transform is defined by(5)w^z=ℒwt=∫0∞e−ztwtdt,provided this integral converges.Lemma 1.
Ifwt∈Cp0,∞, with α∈n−1,n∈Z+, then the Laplace transform of the fractional order Caputo derivative is given by(6)ℒDtαwtz=zαw^−∑i=0n−1zα−i−1wi0.Theorem 1.
the Bromwich inversion theorem [41]). Let w(t) have a continuous derivative and let wt<Keγt, where K and γ are positive constants. Define(7)w^z=∫0∞e−ztwtdt,Rez>γ,then(8)wt=12πi∫ξ−i∞ξ+i∞w^zeztdz.


## 3. Description of the Method

### 3.1. Time Discretization

Here, we apply Laplace transform to models (1)–(3) which gives(9)zI−υ1z1−α+υ2z1−βw^x,z=wx,0−ν1z−α+ν2z−βΔwx,0+f^x,z,x=x,y∈Ω⊂ℝd,ℬw^x,z=h^z,x=x,y∈∂Ω.In more compact form, we have(10)ℒw^x,z=g^x,z,x∈Ω,(11)ℬw^x,z=h^x,z,x∈∂Ω.The transformed problems (10) and (11) will be solved for the solution w^x,z using local RBF method. The solution wx,t of the given models (1)–(3) will be found by using numerical inversion.


### 3.2. Local Radial Basis Functions Method

Here, the linear operatorsℬ and ℒ are discretized by using local RBF [42, 43]. Consider the centers xi,i=1…,N⊂Ω⊂ℝd,d≥1, where Ω is the bounded domain. For each point xi,i=1,2,3,…,N, we can find a subdomain Ωj such that n<N. The unknown function w^x,t can be approximated with RBF in each local subdomain Ωi, i=1,2,…,N, by the following equation:(12)wxi,t≈w^xi,t=∑j=1nλijiϕixi−xj,xj∈Ωj,where λi=λ1i,λ2i,…,λni are the unknown coefficients, and rij=xi−xj is the norm between nodes xi and xj, ϕr, r≥0 is a radial kernel (multiquadric radial basis function), and Ωj⊂Ω is a local domain for around each xi, containing n neighboring nodes around the node xi. So, we have N small size linear systems each of order n×n given by(13)w^1iw^2i⋮w^ni=ϕ11iϕ12i⋯ϕ1niϕ21iϕ22i…ϕ2ni⋮⋮⋱⋮ϕn1iϕn2i…ϕnniλ1iλ2i⋮λni,i=1,2,…,N,which can be denoted by(14)W^i=Siλi,i=1,2,…,N,where ϕjki=ϕixij−xik,xij,xik∈Ωi, and matrix Si is the system matrix.Now, applying the operatorℒ to (12) gives(15)ℒw^xi=∑j=1nλijiℒϕixi−xij,xij∈Ωi.The vector form of (15) is given by(16)ℒw^xi=Gi⋅λi,where Gi is given by(17)Gi=ℒϕixi−xij,xi,xij∈Ωi.From equation (14), the unknown coefficients λi are given by(18)λi=Si−1W^i,and by inserting the values of λi in (16), we have(19)ℒw^xi=GiSi−1W^i=NiW^i,where(20)Ni=GiSi−1.Hence, the discretized form is given by(21)ℒw^≡HW^,where matrix H is called the sparse differentiation matrix of order N×N.


## 3.1. Time Discretization

Here, we apply Laplace transform to models (1)–(3) which gives(9)zI−υ1z1−α+υ2z1−βw^x,z=wx,0−ν1z−α+ν2z−βΔwx,0+f^x,z,x=x,y∈Ω⊂ℝd,ℬw^x,z=h^z,x=x,y∈∂Ω.In more compact form, we have(10)ℒw^x,z=g^x,z,x∈Ω,(11)ℬw^x,z=h^x,z,x∈∂Ω.The transformed problems (10) and (11) will be solved for the solution w^x,z using local RBF method. The solution wx,t of the given models (1)–(3) will be found by using numerical inversion.


## 3.2. Local Radial Basis Functions Method

Here, the linear operatorsℬ and ℒ are discretized by using local RBF [42, 43]. Consider the centers xi,i=1…,N⊂Ω⊂ℝd,d≥1, where Ω is the bounded domain. For each point xi,i=1,2,3,…,N, we can find a subdomain Ωj such that n<N. The unknown function w^x,t can be approximated with RBF in each local subdomain Ωi, i=1,2,…,N, by the following equation:(12)wxi,t≈w^xi,t=∑j=1nλijiϕixi−xj,xj∈Ωj,where λi=λ1i,λ2i,…,λni are the unknown coefficients, and rij=xi−xj is the norm between nodes xi and xj, ϕr, r≥0 is a radial kernel (multiquadric radial basis function), and Ωj⊂Ω is a local domain for around each xi, containing n neighboring nodes around the node xi. So, we have N small size linear systems each of order n×n given by(13)w^1iw^2i⋮w^ni=ϕ11iϕ12i⋯ϕ1niϕ21iϕ22i…ϕ2ni⋮⋮⋱⋮ϕn1iϕn2i…ϕnniλ1iλ2i⋮λni,i=1,2,…,N,which can be denoted by(14)W^i=Siλi,i=1,2,…,N,where ϕjki=ϕixij−xik,xij,xik∈Ωi, and matrix Si is the system matrix.Now, applying the operatorℒ to (12) gives(15)ℒw^xi=∑j=1nλijiℒϕixi−xij,xij∈Ωi.The vector form of (15) is given by(16)ℒw^xi=Gi⋅λi,where Gi is given by(17)Gi=ℒϕixi−xij,xi,xij∈Ωi.From equation (14), the unknown coefficients λi are given by(18)λi=Si−1W^i,and by inserting the values of λi in (16), we have(19)ℒw^xi=GiSi−1W^i=NiW^i,where(20)Ni=GiSi−1.Hence, the discretized form is given by(21)ℒw^≡HW^,where matrix H is called the sparse differentiation matrix of order N×N.


## 4. Numerical Inversion Technique

In this section, the numerical inversion of Laplace transform for approximating the given models (1)–(3) is as follows:(22)wx,t=12πι∫ξ˜−ι∞ξ˜+ι∞w^x,zeztdz=12πι∫Ψeztw^x,zdz,ξ˜>ξ˜0,where Ψ is the suitable path joining ξ˜−ι∞ to ξ˜+ι∞. This Bromwich integral is numerically solved by using the following hyperbolic contour [37]:(23)zη˜=ω+λ˜1−sinσ˜−ιη˜,forη˜∈R,with λ˜>0, ω≥0, 0<σ˜<β˜1/2π, and 1/2π<β˜<π.Integral in (22) gives(24)wx,t=12πι∫−∞∞ezη˜tw^x,zη˜z´η˜dη˜.Next applying trapezoidal rule for approximation of (24), we have(25)wkx,t=k2πι∑j=−MMw^x,zjezjtz´j,zj=zη˜j,η˜j=jk,where k is the step size.


## 5. Application of the Method

In this section, the proposed numerical scheme is applied to multidimensional problems. We solved four test problems and used various domain pointsN∈Ω, stencils points n∈Ωj, and quadrature points M. Three error formulas, the error estimate, Lest=e−cM/logM, L∞, and L2 norms are used. The radial kernel used in our computations is ϕr,ɛ=1+ɛ2r2. The shape parameter ϵ is optimized by the uncertainty rule related to RBFs.Problem 1.
Consider models (1)–(3) to the following form [38]:(26)∂wx,t∂t=D0t1−α+D0t1−β∂2wx,t∂x2+fx,t,fx,t=expx1+αtα−Γ2+α+βΓ1+2α+βt2α+β−Γ2+α+βΓ1+α+2βtα+2β,with the following boundary and initial conditions:(27)w0,t=t1+α+β,w1,t=et1+α+β,0<t≤T,wx,0=0,x∈0,1,respectively, where the actual solution is given by(28)wx,t=expxt1+α+β.
In our numerical scheme, we used the hyperbolic contour (23). The optimal parameter values are taken as(29)λ=θrbNbT,b=cosh−11θτsinσ,ω=2,σ=0.3812,xk=hk,rb=2πr, r=0.3431, h=b/N, τ=t0/T, t0=0.5, t=1, and T=5. This test problem is solved in the domain 0,1. Here, the number of points in domain Ω is denoted by N, the points in local subdomain ∂Ωj are denoted by n, and the number of quadrature points relates to M. The numerical solutions are shown in Table 1 with various values of fractional order α and β and nodal points N. For comparatively smaller values of fractional order α and β, better results in terms of L∞ and L2 error norms are obtained. In the upper part of Table 1, condition number increases, as we increase nodal points N. Error versus various quadrature points M at N=21,n=9, and t=1 and various values of α and β are shown in Figure 1. The error estimate Lest for c=1 is well matched with L∞ and L2 error norms, as shown in Figure 1. Hence, our proposed method is stable and accurate.Table 1
Numerical results using the proposed numerical scheme corresponding to Problem1.
M=50, n=7α=0.2, β=0.1NL∞L2κ118.9603e − 0040.00212.2091e + 021218.5582e − 0040.00289.6329e + 021318.4544e − 0040.00342.2821e + 022418.8689e − 0040.00414.4923e + 022518.6646e − 0040.00457.5954e + 022719.2962e − 0040.00571.6860e + 023M=50, n=7,N=21L∞L2κα,β=0.2,0.60.01410.04621.0025e + 021α,β=0.2,0.40.00720.02363.6776e + 020α,β=0.2,0.18.5582e − 0040.00289.6329e + 021α,β=0.6,0.30.00790.02583.1036e + 020α,β=0.4,0.30.00630.02073.1036e + 020α,β=0.1,0.30.00350.01169.2023e + 021Figure 1
Error versus quadrature pointsM at t=1 and various values of α,β corresponding to Problem 1.
(a)(b)(c)Problem 2.
Consider models (1)–(3) corresponding to the form [38](30)∂wx,t∂t=12∂1−αwx,t∂t1−α+∂1−βwx,t∂t1−β∂2wx,t∂x2+fx,t,fx,t=expx1+αtα−Γ2+αΓ1+2αt2α+1+βtβ−Γ2+βΓ1+2βt2β,initial and boundary conditions given by(31)wx,0=0,x∈0,1,w0,t=t1+α+t1+β,w1,t=exp1t1+α+t1+β,0<t≤1.
The actual solution is(32)wx,t=expxt1+α+t1+β.
The same domain and same parameter values as used in Problem1 are incorporated. The numerical results are shown in Table 2 with the same as well as with various values of fractional order α and β and nodal points N. For comparatively identical values of fractional order α and β, better results in terms of L∞ and L2 error norms are obtained. In the upper part of Table 2, condition number of the system matrix is fixed for 11≤N≤71. Error versus various quadrature points M at N=41,n=9, and t=1 and various values of α and β are depicted in Figure 2. The error estimate Lest for c=0.7 is well agreed with L∞ and L2error norms, as shown in Figure 1. The results obtained by our proposed numerical scheme are comparatively identical with the results in Table 2 [38].Table 2
Numerical results using the proposed numerical scheme corresponding to Problem2.
M=80, n=9α=0.5, β=0.5NL∞L2κ115.4148e − 0051.2209e − 0045.0986e + 033211.3432e − 0042.7509e − 0045.0985e + 033311.1372e − 0043.2761e − 0045.0985e + 033418.6091e − 0053.1340e − 0045.0985e + 033511.2453e − 0045.3029e − 0045.0985e + 033714.5906e − 0052.1446e − 0045.0985e + 033M=80, n=9,N=41L∞L2κα,β=0.2,0.60.00720.03315.1040e + 033α,β=0.2,0.40.00190.00895.1040e + 033α,β=0.2,0.15.1121e − 0040.00232.6797e + 035α,β=0.6,0.30.00390.01795.0985e + 033α,β=0.4,0.34.3443e − 0040.00195.0985e + 033α,β=0.1,0.30.00210.00962.1983e + 034Figure 2
Error versus quadrature pointsM at N=41, n=9, and t=1 and various values of α,β corresponding to Problem 2.
(a)(b)(c)Problem 3.
Next, we consider models (1)–(3) corresponding to the form [44](33)∂wx,y,t∂t=D0t1−α+D0t1−β∂2wx,y,t∂x2+∂2wx,y,t∂y2+fx,t,where(34)fx,y,t=2tsin2πxsin2πy1+8π.2Γ2+αtα+8π.2Γ2+βtβ,initial and boundary conditions given by(35)wx,y,0=0,x,y∈Ω,w0,t=0,w1,t=0,t>0.
The exact solution is(36)wx,y,t=t2sin2πxsin2πy.
This problem is solved over the domainΩ=0,1. In Table 3, for various nodal points N and stencils points n=11,15 and with various values of α and β, the L∞ error norm is well matched with L2 error norm. The condition number is increasing steadily as we decrease both the values of α and β at the same time.Table 3
Numerical results using the proposed numerical scheme corresponding to Problem 3.
M=50, n=11α=0.5, β=0.5NL∞L2κ157.3893e − 0040.00503.1023e + 020200.00250.02323.1024e + 020260.00290.03443.0989e + 020418.6091e − 0053.1340e − 0045.0985e + 033M=50, n=15,N=20L∞L2κα,β=0.6,0.52.2399e − 0040.00173.1000e + 020α,β=0.5,0.32.2409e − 0040.00175.0520e + 021α,β=0.3,0.22.2419e − 0040.00171.1756e + 024Problem 4.
Finally, we consider models (1)–(3) corresponding to the form [38](37)∂wx,y,t∂t=D0t1−α+D0t1−β∂2wx,y,t∂x2+∂2wx,y,t∂y2+fx,t,0<t≤1,0<x,y<1,fx,y,t=exp−x−0.52γ−y−0.52γf1x,y,t+f2x,y,t,where(38)f1x,y,t=1+α+βtα+β+2Γ2+α+βΓ1+2α+βt2α+β2γ+2Γ2+α+βΓ1+α+2βtα+2β2γ,f2x,y,t=−4Γ2+α+βΓ1+2α+βt2α+β+Γ2+α+βΓ1+α+2βtα+2βx−0.52γ2+y−0.52γ2.
The exact solution is(39)wx,y,t=t1+α+βexp−x−0.52γ−y−0.52γ.
Here, the problem is solved over the domainΩ=0,1×0,1. In the upper section of Table 4, the ℓ∞ and ℓ2 error norms are decreasing with α=0.5, β=0.3, n=9, and γ=0.2 and for nodal points 11≤N≤25. In the lower section of Table 4, for same values of α and β, the ℓ∞ and ℓ2 error norms are decreasing steadily at N=20, n=11, and M=50. The results are comparatively identical with the results of the paper [38]. Figure 3 shows the error with varying quadrature points M and various values of α and β at N=20, M=50, and γ=0.2. The error L∞ is well matched with estimate Lest for c=0.5 and L2 error norm, as shown in Figure 3. The present method is stable and accurate in multidimensional fractional order partial differential equations.Table 4
Numerical results using the proposed numerical scheme corresponding to Problem4.
M=50, n=9α=0.5, β=0.3, γ=0.2NL∞L2κ110.00690.02491.3362e + 021150.00220.01022.7478e + 021214.4497e − 0040.00576.7578e + 021257.8109e − 0040.01121.3745e + 022M=50,N=20,n=11L∞L2κα,β=0.2,0.26.2477e − 0040.00609.6487e + 025α,β=0.5,0.56.2141e − 0040.00603.1021e + 020α,β=0.7,0.76.1740e − 0040.00593.0989e + 020α,β=0.9,0.96.1134e − 0040.00583.0989e + 020Figure 3
Error versus quadrature pointsM at N=20, n=11, and t=1 and various values of α,β corresponding to Problem 4.
(a)(b)


## 6. Conclusion

In this work, a numerical scheme is constructed which is based on Laplace transform and radial basis functions in the local setting. The proposed numerical scheme efficiently approximated time fractional anomalous subdiffusion equation. The supremacy of this method particularly for fractional order equations is its nonsensitive nature in time as contrary to finite difference approximation for fractional order operators. Since the fractional order derivative is of integral convolution type and suited to handle by Laplace transform, the spatial operators in multidimensions can be approximated by RBF in the local setting which generates small size differentiation matrices in local subdomains and these are assembled as a single sparse matrix in the global domain. So, large amount of data can be manipulated very easily and accurately.



---

*Source: 9965734-2021-07-09.xml*