# A Numerical Method for Time-Fractional Reaction-Diffusion and Integro Reaction-Diffusion Equation Based on Quasi-Wavelet

Authors: Sachin Kumar; Jinde Cao; Xiaodi Li
Journal: Complexity
 (2020)

---

## Abstract

In this research work, we focused on finding the numerical solution of time-fractional reaction-diffusion and another class of integro-differential equation known as the integro reaction-diffusion equation. For this, we developed a numerical scheme with the help of quasi-wavelets. The fractional term in the time direction is approximated by using the Crank–Nicolson scheme. The spatial term and the integral term present in integro reaction-diffusion are discretized and approximated with the help of quasi-wavelets. We study this model with Dirichlet boundary conditions. The discretization of these initial and boundary conditions is done with a different approach by the quasi-wavelet-based numerical method. The validity of this proposed method is tested by taking some numerical examples having an exact analytical solution. The accuracy of this method can be seen by error tables which we have drawn between the exact solution and the approximate solution. The effectiveness and validity can be seen by the graphs of the exact and numerical solutions. We conclude that this method has the desired accuracy and has a distinctive local property.


---

## Body

## 1. Introduction

Differential and integral calculus is a vital and broad branch of mathematics [1, 2]. First of all, the integer-order differential and integer calculus came into existence. About a few decades ago, the concept of fractional calculus came into existence. Day by day the new concept, new definitions, and new terminology are emerging. The concepts and ideas of fractional calculus were evolved by Liouville and Abel. The more basic fundamentals and concepts are available in the literature [3, 4]. This branch of mathematics allows us to generalize the order of derivatives from the integer order to the real and arbitrary order. There are so many types available in the literature of fractional calculus. Some of them are Caputo, Hadamard, Riemann–Liouville, and Grunwald–Letnikov. As time passes, in the development of this branch, the variable-order derivative concepts are introduced in which orders of derivatives are functions of variables. The most prevailing derivatives are power law kernel derivatives. But, in recent few years, the derivatives with nonsingular kernel are introduced. Derivatives with the Mittag–Leffler kernel are known as Atangana–Baleanu derivatives, and derivatives with the exponential kernel are known as the Caputo–Fabrizio derivative. The fractional differential equations are in trends to study. There are so many physical phenomena in which fractional differential equations depict and represent the behavior and dynamics of the phenomena in a more accurate way than integer-order differential equations. To find out the solution of the fractional differential equation and fractional partial differential equations is a tough task as it is not possible always to find out the exact solution. The researchers started to develop the methods and schemes for finding the numerical solution of fractional differential equations. Many numerical methods are developed and available in the literature. Some of them are the fractional differential transform method [5], predictor-corrector method [6], Adomian decomposition method [7], homotopy perturbation method [8], and generalized block pulse operational matrix method [9]. The operational matrix method is a spectral method which is easy to be applicable for many types of fractional differential equations whether they are of integer, fractional, and variable type. Its results are so accurate, and it is easy to apply. Some of them are Gegenbauer wavelet [9], Haar wavelets [10], Chebyshev wavelets [11], and Legendre wavelets [12]. The stability of discrete-time delayed switched systems and time delayed systems is given in the literature [13, 14].FDEs, FPDEs, and fractional integro-differential equations are used to model many physical phenomena [15–17], which have arisen in physics, chemistry, biology, electrical engineering, nuclear reactor, viscoelastic, fluid dynamics, population dynamics, and medical science. To show the memory effect of the system, in the branch of nuclear reactor and thermo elasticity, we use the fractional differential equations. But, while modeling, the effect of history is being ignored. This can be taken into account by adding a term of integration in these fractional PDEs. The integro-differential equations have so many applications in control theory, financial market, fluid dynamics, aerospace systems, and viscoelastic mechanism. The most occurring physical process in nature is the diffusion process. In this process of diffusion, the molecules or energy or any physical quantity is transferred to the low concentration area from the high concentration area. If another process known as reaction is added to the process of diffusion, then it is called the reaction-diffusion process. In this process, molecules react with the fluids. This phenomenon of diffusion can be depicted by the following PDE known as the reaction-diffusion partial differential equation:(1)∂Ψζ,η∂η=D∇2Ψζ,η+PΨ,ζ,η,where PΨ,η represents the reaction term and the first term on the right-hand side of the equation represents the diffusion process with D as the diffusion coefficient.We divided our article into many sections as follows: Section1 contains the preliminary definitions which will be used throughout the article. This section contains the theory and methods of the quasi-wavelet to find out an approximation of an unknown function. The taken model of integro reaction-diffusion and reaction-diffusion is given in Section 3. How to deal with discretization of spatial and time-fractional terms with Dirichlet boundary conditions is also given in Section 3. The numerical results which depict the effectiveness and accuracy of our proposed method are present in Section 4. The conclusion of this paper is contained in Section 5.


## 2. Preliminary Definition about Fractional Derivative and Quasi-Wavelets

In this section, we described the definitions of fractional derivatives which will be used throughout the article. There are two types of definitions available in the literature of the fractional derivative, i.e., Riemann–Liouville and Caputo definitions. We have taken Caputo’s one because it is more reliable in modeling.


### 2.1. Integration and Differentiation of Riemann–Liouville Type

The R-L-type definition of fractional integration is given as follows [18]:(2)Iϱζθ=1Γϱ∫0θθ−Ωϱ−1ζΩdΩ,θ>0,ϱ∈R+.The definition of fractional differentiation in the R-L sense is given by the following definition:(3)Dlϱζθ=ddθmIm−ϱζθ,ϱ>0,m−1<ϱ<m.


### 2.2. Definition of Fractional Derivative in Caputo Sense

The following equation represents the definition of fractional differentiation in Caputo sense:(4)Dcϱζθ=dqζθdθq,ϱ=q∈N,1Γϱ∫0θθ−ηq−ϱ−1ζqηdη,q−1<ϱ<q,where the parameter q represents the integer. A resemblance of the Caputo definition with the integer-order derivative is as follows:(5)DcϱM=0,with M being a constant.(6)Dcϱθσ=0,σ∈N∪0andσ<ϱ,Γ1+σΓ1−ϱ+σθ−ϱ+σ,σ∈N∪0andσ≥ϱorσ∉Nandσ>ϱ,where ⌊ϱ⌋ represents the floor function. All these types of fractional operators follow the linear property:(7)DcϱC1h1θ+C2h2θ=C1Dcϱh1θ+C2Dcϱh2θ,with C1 and C2 being constants. A useful relation between the R-L definition and Caputo is given as follows:(8)IϑDcϑςη=ςη−∑k=0l−1ςk0+ηkk!,l−1<ϑ≤l.


### 2.3. Approximation of Functions with the Help of Quasi-Wavelets

The method based upon quasi-wavelets nowadays is developing in the form of the local spectral collocation method. It is used for seeking the numerical solution of the fractional-order partial differential equation and integral equations. The singular convolution used in this method is a significant transformation. We define this transformation as follows:(9)ϖθ=R∗pθ=∫−∞∞Rθ−ζpζdζ,where pθ is known as the test function and R is a singular kernel. We know that the construction of the family of wavelets is done with the help of a function known as the mother wavelet φ. We use the operations of dilation and transformation to form the family of wavelets:(10)φϑ,γζ=ϑ−1/2φζ−γϑ.We have used two parameters here. The parameterϑ is the dilation parameter, and γ is the translation parameter. The kernel known as Shannon’s delta sequence is used in the construction of quasi-wavelets. It is defined as(11)δθη=1π∫0πcosηydy=sinθηπη,where limθ⟶θ0δθη=δη and δ is a well-known Dirac delta function. Shannon’s delta sequence kernel produces the basis for the Paley–Wiener reproducing kernel Hilbert space Bθ2 [19] for a given θ>0 which is a subspace of L2R. We can produce uniquely a function ςζ∈Bθ2 as follows:(12)ςζ=∫−∞∞ςζδθζ−ηdt=∫−∞∞ςζsinθζ−ηπζ−ηdη,∀ς∈Bθ2.The sampling scaling function can be written in another form as follows:(13)δθ,p=δθζ−ζp=sinζ−ζpθζ−ζpπ.We have denoted the sampling points around the spaceζ by ζp which are sampling points centered around ζ. The discrete form of function ςζ which belongs to the space Bθ2 is presented by the following infinite series expansion:(14)ςζ=∑p=−∞∞ςζpδθζ−ζp.The above approximation can be rewritten with the help of Shannon sampling theorem at the Nyquist frequencyθ in the band-limited signal in Bθ2:(15)ςζ=∑p=−∞∞ςζpδθζ−ζp=∑p=−∞∞ςζpsinπζ−ζp/Δπζ−ζp/Δ,where the value of parameter θ is π/Δ and the grid size in spatial direction is denoted by Δ. We improved the localized asymptotic behavior of Dirichlet’s delta sequence kernel by using the regularizer Rσζ. In the literature, many regularizers are introduced, and we use here the following Gaussian type regularizer:(16)Rσζ=exp−ζ22σ2,σ>0,where σ represents the width parameter of the Gaussian envelope. The parameters Δ and σ are related by the expression σ=r×Δ, with r as the computation parameter. The sampling scaling function is modified as follows by using the Gaussian regularizer:(17)δΔ,σζ=sinπζ/Δπζ/Δexp−ζ22σ2,where(18)limσ⟶∞δΔ,σζ=sinπζ/Δπζ/Δ,where the Gaussian-regularized sampling scaling function follows the property of the quasi-scaling function. So, we can approximate an unknown function ςζ belonging to the space textbfBθ2 with the help of the quasi-scaling function:(19)ςζ=∑p=−∞∞ςζpδθζ−ζp=∑p=−∞∞ςζpδθζ−ζpRθζ−ζp.Here, the infinite sampling points have been used, but while commutating in practical manner, we cannot use these infinite points. By seeing this problem, the domain of computation is restricted to the finite sampling points. For this purpose, we select2W+1. The simplified form of the above equation is as follows:(20)ςζ=∑p=−WWςζpδΔ,σζ−ζp.Now, for the approximation of integer-order derivatives of ordern of the function fζ,(21)ςnζ=∑p=−WWςζpδΔ,σnζ−ζp,n=1,2,….The sampling points vary fromp=−W to p=+W, i.e., there are 2W+1 sampling points. We give some of the formula of δΔ,σ,δΔ,σ1, and δΔ,σ2 [20] for the calculation purpose:(22)δΔ,σζ=exp−ζ2/2σ2sinζπ/Δπζ/Δ,ζ≠01,ζ=0,δΔ,σ1ζ=−sinπζ/Δπζ2/Δ−Δsinπζ/Δπσ2/Δ+cosπζ/Δζexp−ζ22σ2,ζ≠0,0,ζ=0,δΔ,σ2ζ=2Δsinζπ/Δπζ3−2cosπζ/Δζ2+Δζsinζπ/Δπσ4+Δsinπζ/Δπσ2ζ−2cosπζ/Δσ2−πsinζπ/ΔζΔexp−ζ22σ2,ζ≠0,0,ζ=0.


## 2.1. Integration and Differentiation of Riemann–Liouville Type

The R-L-type definition of fractional integration is given as follows [18]:(2)Iϱζθ=1Γϱ∫0θθ−Ωϱ−1ζΩdΩ,θ>0,ϱ∈R+.The definition of fractional differentiation in the R-L sense is given by the following definition:(3)Dlϱζθ=ddθmIm−ϱζθ,ϱ>0,m−1<ϱ<m.


## 2.2. Definition of Fractional Derivative in Caputo Sense

The following equation represents the definition of fractional differentiation in Caputo sense:(4)Dcϱζθ=dqζθdθq,ϱ=q∈N,1Γϱ∫0θθ−ηq−ϱ−1ζqηdη,q−1<ϱ<q,where the parameter q represents the integer. A resemblance of the Caputo definition with the integer-order derivative is as follows:(5)DcϱM=0,with M being a constant.(6)Dcϱθσ=0,σ∈N∪0andσ<ϱ,Γ1+σΓ1−ϱ+σθ−ϱ+σ,σ∈N∪0andσ≥ϱorσ∉Nandσ>ϱ,where ⌊ϱ⌋ represents the floor function. All these types of fractional operators follow the linear property:(7)DcϱC1h1θ+C2h2θ=C1Dcϱh1θ+C2Dcϱh2θ,with C1 and C2 being constants. A useful relation between the R-L definition and Caputo is given as follows:(8)IϑDcϑςη=ςη−∑k=0l−1ςk0+ηkk!,l−1<ϑ≤l.


## 2.3. Approximation of Functions with the Help of Quasi-Wavelets

The method based upon quasi-wavelets nowadays is developing in the form of the local spectral collocation method. It is used for seeking the numerical solution of the fractional-order partial differential equation and integral equations. The singular convolution used in this method is a significant transformation. We define this transformation as follows:(9)ϖθ=R∗pθ=∫−∞∞Rθ−ζpζdζ,where pθ is known as the test function and R is a singular kernel. We know that the construction of the family of wavelets is done with the help of a function known as the mother wavelet φ. We use the operations of dilation and transformation to form the family of wavelets:(10)φϑ,γζ=ϑ−1/2φζ−γϑ.We have used two parameters here. The parameterϑ is the dilation parameter, and γ is the translation parameter. The kernel known as Shannon’s delta sequence is used in the construction of quasi-wavelets. It is defined as(11)δθη=1π∫0πcosηydy=sinθηπη,where limθ⟶θ0δθη=δη and δ is a well-known Dirac delta function. Shannon’s delta sequence kernel produces the basis for the Paley–Wiener reproducing kernel Hilbert space Bθ2 [19] for a given θ>0 which is a subspace of L2R. We can produce uniquely a function ςζ∈Bθ2 as follows:(12)ςζ=∫−∞∞ςζδθζ−ηdt=∫−∞∞ςζsinθζ−ηπζ−ηdη,∀ς∈Bθ2.The sampling scaling function can be written in another form as follows:(13)δθ,p=δθζ−ζp=sinζ−ζpθζ−ζpπ.We have denoted the sampling points around the spaceζ by ζp which are sampling points centered around ζ. The discrete form of function ςζ which belongs to the space Bθ2 is presented by the following infinite series expansion:(14)ςζ=∑p=−∞∞ςζpδθζ−ζp.The above approximation can be rewritten with the help of Shannon sampling theorem at the Nyquist frequencyθ in the band-limited signal in Bθ2:(15)ςζ=∑p=−∞∞ςζpδθζ−ζp=∑p=−∞∞ςζpsinπζ−ζp/Δπζ−ζp/Δ,where the value of parameter θ is π/Δ and the grid size in spatial direction is denoted by Δ. We improved the localized asymptotic behavior of Dirichlet’s delta sequence kernel by using the regularizer Rσζ. In the literature, many regularizers are introduced, and we use here the following Gaussian type regularizer:(16)Rσζ=exp−ζ22σ2,σ>0,where σ represents the width parameter of the Gaussian envelope. The parameters Δ and σ are related by the expression σ=r×Δ, with r as the computation parameter. The sampling scaling function is modified as follows by using the Gaussian regularizer:(17)δΔ,σζ=sinπζ/Δπζ/Δexp−ζ22σ2,where(18)limσ⟶∞δΔ,σζ=sinπζ/Δπζ/Δ,where the Gaussian-regularized sampling scaling function follows the property of the quasi-scaling function. So, we can approximate an unknown function ςζ belonging to the space textbfBθ2 with the help of the quasi-scaling function:(19)ςζ=∑p=−∞∞ςζpδθζ−ζp=∑p=−∞∞ςζpδθζ−ζpRθζ−ζp.Here, the infinite sampling points have been used, but while commutating in practical manner, we cannot use these infinite points. By seeing this problem, the domain of computation is restricted to the finite sampling points. For this purpose, we select2W+1. The simplified form of the above equation is as follows:(20)ςζ=∑p=−WWςζpδΔ,σζ−ζp.Now, for the approximation of integer-order derivatives of ordern of the function fζ,(21)ςnζ=∑p=−WWςζpδΔ,σnζ−ζp,n=1,2,….The sampling points vary fromp=−W to p=+W, i.e., there are 2W+1 sampling points. We give some of the formula of δΔ,σ,δΔ,σ1, and δΔ,σ2 [20] for the calculation purpose:(22)δΔ,σζ=exp−ζ2/2σ2sinζπ/Δπζ/Δ,ζ≠01,ζ=0,δΔ,σ1ζ=−sinπζ/Δπζ2/Δ−Δsinπζ/Δπσ2/Δ+cosπζ/Δζexp−ζ22σ2,ζ≠0,0,ζ=0,δΔ,σ2ζ=2Δsinζπ/Δπζ3−2cosπζ/Δζ2+Δζsinζπ/Δπσ4+Δsinπζ/Δπσ2ζ−2cosπζ/Δσ2−πsinζπ/ΔζΔexp−ζ22σ2,ζ≠0,0,ζ=0.


## 3. Mathematical Model and Used Proposed Method

In this section, we propose the integrated mathematical model of the reaction-diffusion and integro reaction-diffusion equation. In this section, we develop the Crank–Nicolson scheme with the combination of the quasi-wavelet-based numerical method to find the numerical solution of the taken model:(23)∂ϑς∂ηϑ=∂2ς∂ζ2+aςη,ζ1−ςη,ζ+bςη,ζ∫0ηκη−s,ζςs,ζds+fζ,η.We have taken the following Dirichlet boundary conditions and initial conditions:(24)ς0,η=g1η,ς1,η=g2η,ςζ,0=g3ζ,where the term fζ,η denotes the forced function and κζ,η is known as a kernel.


### 3.1. Using Crank–Nicolson Scheme for the Time Discretization

To discretize the time derivatives present in model, we will use the Crank–Nicolson scheme. We divided the time interval into discrete points asηn=n×Δη where n=0,1,…,M and time length is noted by Δη. The notations ςn and fn represent the values of ςζ,η and fζ,η at time ηn=n×Δη. The term ∂ϑςζ,η/∂ηϑ from the model is approximated as at the grid point ζ,ηn+1/2:(25)∂ϑςζ,ηn+1/2∂ηϑ=1Γ1−ϑ∫0ηn+1/2ηn+1/2−s−ϑ×∂ςζ,s∂sds=1Γ1−ϑ∫0ηnn+12Δη−s−ϑ×∂ςxi,s∂sds+∫ηnηηn+1/2ςn+1−ςnΔη+OΔηn+12Δη−s−ϑds=1Γ1−ϑ∑m=1n∫m−1ΔηmΔηn+12Δη−s−ϑ×ςm+1−ςmΔη+s−ηm−1/2uttcm,ζds+1Γ1−ϑ∫nΔηn+1/2Δηςn+1−ςnΔη+OΔηn+12Δη−s−ϑds.In aberration form, we can rewrite as(26)∂ϑςζ,ηn+1/2∂ηϑ=Δη−ϑΓ2−ϑ∑m=1nςm+1−ςmn+32−m1−ϑ+n+12−m1−ϑ+Δη−ϑΓ2−ϑςn+1−ςn21−ϑ+1Γ1−ϑ∑m=1n∫m−1ΔηmΔηn+12Δη−s−ϑ×s−ηm−1/2uttcm,ζds+Δη−ϑΓ2−ϑ121−ϑOΔη2−ϑ.Hence, we obtain(27)∂ϑςxi,ηn+1/2∂ηϑ=Δη−ϑΓ2−ϑ∑m=1nuim+1−uimn+32−m1−ϑ+n+12−m1−ϑ+uin+1−uin21−ϑ+R1+R2,where(28)R1=1Γ1−ϑ∑m=1n∫m−1ΔηmΔηn+12Δη−s−ϑ×s−ηm−1/2uttcm,xids,R2=2ϑ−1Γ2−ϑOΔη2−ϑ.The bound ofR1 is as follows after simplifying:(29)R1≤max0≤ζ≤1uttcm,ζΓ1−ϑΔη2−ϑ.The trapezoidal technique is implemented to deal with the integral term∫0ηκη−s,ζςs,ζds:(30)∫0ηn+1/2κζ,η−sςζ,sds=12∫0ηnκζ,η−sςζ,sds+∫0ηn+1κζ,η−sςζ,sds.Taking the exponential form of the kernel asκη−s,ζ=e−η−s, then(31)∫0ηnκζ,η−sςζ,sds=∑j=0n−1∫ηjηj+1exp−ηn−sςζ,sds=∑j=0n−1∫ηjηj+1exp−ηn−sςζ,ηj+1s−ηjΔη+ςζ,ηjηj+1−sΔηds=e−ηnΔη∑j=0n−1ςj+1eηj+1Δη−eηj+1+eηj+ςj−eηjΔη+eηj+1−eηj.Similarly,(32)∫0ηn+1κζ,η−sςζ,sds=e−ηn+1Δηςn+1eηn+1Δη−eηn+1+eηn+ςn−eηnΔη+eηn+1−eηn+e−ηn+1Δη∑j=0n−1ςj+1eηj+1Δη−eηj+1+eηj+ςj−eηjΔη+eηj+1−eηj.The temporal semidiscrete form is obtained with the help of equations (32)–(34):(33)Δη−ϑΓ2−ϑ∑m=1nςm+1−ςmn+32−m1−ϑ+n+12−m1−ϑ+ςn+1−ςn21−ϑ=bςnζe−ηn+e−ηn+12Δη∑j=0n−1ςj+1eηj+1Δη−eηj+1+eηj+ςj−eηjΔη+eηj+1−eηj+bςnζe−ηn+12Δηςn+1eηn+1Δη−eηn+1+eηn+ςn−eηnΔη+eηn+1−eηn+∂2ςnζ∂ζ2+aςnζ1−ςnζ.


### 3.2. Use of Quasi-Wavelet-Based Approximation for the Discretization of Spatial Derivatives

By using the approximation by quasi-wavelets given in Section2, we discretize the unknown function and its spatial derivatives. The spatial step is chosen as Δζ=1/N. We denote the value of the unknown function ςζ,η at spatial point ζi and time η=ti by ςin with i=0,1,…,N and n=0,1,…,M. By using the approximation from equation (21), we have(34)ςnxi=∑p=i−Wi+Wςζk,ηnδΔ,σnζi−ζk,n=0,1,2,…,i=0,1,…,N−1.With the help of equation (33), we can rewrite equation (34) as follows:(35)Δη−ϑΓ2−ϑ∑m=1nςm+1−ςmn+32−m1−ϑ+n+12−m1−ϑ+ςn+1−ςn21−ϑ=b∑k=i−Wi+Wςζk,ηnδΔ,σζi−ζke−tηn+e−ηn+12Δη∑j=0n−1ςj+1eηj+1Δη−eηj+1+eηj+ςj−eηjΔη+eηj+1−eηj+b∑p=i−Wi+Wςζk,ηnδΔ,σζi−ζke−ηn+12Δηςn+1eηn+1Δη−eηn+1+eηn+ςn−eηnΔη+eηn+1−eηn+∑p=i−Wi+Wςζk,ηnδΔ,σ2ζi−ζk+a∑p=i−Wi+Wςζk,ηnδΔ,σζi−ζk1−∑p=i−Wi+Wςζk,ηnδΔ,σζi−ζk.Consideringxp−xi=xk, we have the following:(36)Δη−ϑΓ2−ϑ∑m=1nςim+1−ςimn+32−m1−ϑ+n+12−m1−ϑ+ςin+1−ςin21−ϑ=b∑k=−W+Wuk+inδΔ,σ−kΔζ×e−ηn+e−ηn+12Δη∑j=0n−1ςij+1eηj+1Δη−eηj+1+eηj+ςij−eηjΔη+eηj+1−eηj+b∑k=−W+Wuk+inδΔ,σ−kΔζ×e−ηn+12Δηςin+1eηn+1Δη−eηn+1+eηn+ςin−eηnΔη+eηn+1−eηn+∑k=−W+Wuk+inδΔ,σ2ζi−ζk+a∑k=−Wi+Wuk+inδΔ,σ−kΔζ1−∑k=−W+Wuk+inδΔ,σ−kΔζ.This is the full discrete form of the taken model in case of the kernelκζ,η=e−η. For other kernels κζ,η=1 and κζ,η=e−ζ2t, the discrete form is obtained as follows:(37)Δη−ϑΓ2−ϑ∑m=1nςim+1−ςimn+32−m1−ϑ+n+12−m1−ϑ+ςin+1−ςin21−ϑ=bΔη4∑j=0n−1uij+ςij+1+Δη4∑j=0nuij+ςij+1+∑k=−W+Wuk+inδΔ,σ2ζi−ζk+a∑k=−Wi+Wuk+inδΔ,σ−kΔζ1−∑k=−W+Wuk+inδΔ,σ−kΔζ,Δη−ϑΓ2−ϑ∑m=1nςim+1−ςimn+32−m1−ϑ+n+12−m1−ϑ+ςin+1−ςin21−ϑ=b∑k=−W+Wuk+inδΔ,σ−kΔζ×e−ηn+1xi22Δη∑j=0n−1ςij+1eηj+1xi2xi2Δη−eηj+1xi2xi4+eηjxi2xi4+uij−eηj+1xi2xi2Δη+eηj+1xi2xi4−eηjxi2xi4+∑k=−W+Wuk+inδΔ,σ2ζi−ζk+a∑k=−Wi+Wuk+inδΔ,σ−kΔζ1−∑k=−W+Wuk+inδΔ,σ−kΔζ.The discretization of initial and boundary conditions is performed by adopting a technique because the value ofςζk is not defined outside the domain 0,1:(38)ςi,n=ςin=0,i>M,i<0,n=0,1,…,N,ui0=f3ζi,i=0,1,…,M.


### 3.3. Theorem (see [21])

Considering the functionfζ which is a member of the space L∞∩L2Ω∩CsΩ and band limited to B, s∈Z+, σ=Δr, W∈N, W≥rs/2. Then, we have the following required result:(39)fs−∑k=−WWδσ,Δsζ−ζkfζk≤β×exp−γ22r2,where(40)γ=minr2π−BΔ,W,β=2BfLsΩ+2rfL∞Ω×eπs+1!rγπΔs.


## 3.1. Using Crank–Nicolson Scheme for the Time Discretization

To discretize the time derivatives present in model, we will use the Crank–Nicolson scheme. We divided the time interval into discrete points asηn=n×Δη where n=0,1,…,M and time length is noted by Δη. The notations ςn and fn represent the values of ςζ,η and fζ,η at time ηn=n×Δη. The term ∂ϑςζ,η/∂ηϑ from the model is approximated as at the grid point ζ,ηn+1/2:(25)∂ϑςζ,ηn+1/2∂ηϑ=1Γ1−ϑ∫0ηn+1/2ηn+1/2−s−ϑ×∂ςζ,s∂sds=1Γ1−ϑ∫0ηnn+12Δη−s−ϑ×∂ςxi,s∂sds+∫ηnηηn+1/2ςn+1−ςnΔη+OΔηn+12Δη−s−ϑds=1Γ1−ϑ∑m=1n∫m−1ΔηmΔηn+12Δη−s−ϑ×ςm+1−ςmΔη+s−ηm−1/2uttcm,ζds+1Γ1−ϑ∫nΔηn+1/2Δηςn+1−ςnΔη+OΔηn+12Δη−s−ϑds.In aberration form, we can rewrite as(26)∂ϑςζ,ηn+1/2∂ηϑ=Δη−ϑΓ2−ϑ∑m=1nςm+1−ςmn+32−m1−ϑ+n+12−m1−ϑ+Δη−ϑΓ2−ϑςn+1−ςn21−ϑ+1Γ1−ϑ∑m=1n∫m−1ΔηmΔηn+12Δη−s−ϑ×s−ηm−1/2uttcm,ζds+Δη−ϑΓ2−ϑ121−ϑOΔη2−ϑ.Hence, we obtain(27)∂ϑςxi,ηn+1/2∂ηϑ=Δη−ϑΓ2−ϑ∑m=1nuim+1−uimn+32−m1−ϑ+n+12−m1−ϑ+uin+1−uin21−ϑ+R1+R2,where(28)R1=1Γ1−ϑ∑m=1n∫m−1ΔηmΔηn+12Δη−s−ϑ×s−ηm−1/2uttcm,xids,R2=2ϑ−1Γ2−ϑOΔη2−ϑ.The bound ofR1 is as follows after simplifying:(29)R1≤max0≤ζ≤1uttcm,ζΓ1−ϑΔη2−ϑ.The trapezoidal technique is implemented to deal with the integral term∫0ηκη−s,ζςs,ζds:(30)∫0ηn+1/2κζ,η−sςζ,sds=12∫0ηnκζ,η−sςζ,sds+∫0ηn+1κζ,η−sςζ,sds.Taking the exponential form of the kernel asκη−s,ζ=e−η−s, then(31)∫0ηnκζ,η−sςζ,sds=∑j=0n−1∫ηjηj+1exp−ηn−sςζ,sds=∑j=0n−1∫ηjηj+1exp−ηn−sςζ,ηj+1s−ηjΔη+ςζ,ηjηj+1−sΔηds=e−ηnΔη∑j=0n−1ςj+1eηj+1Δη−eηj+1+eηj+ςj−eηjΔη+eηj+1−eηj.Similarly,(32)∫0ηn+1κζ,η−sςζ,sds=e−ηn+1Δηςn+1eηn+1Δη−eηn+1+eηn+ςn−eηnΔη+eηn+1−eηn+e−ηn+1Δη∑j=0n−1ςj+1eηj+1Δη−eηj+1+eηj+ςj−eηjΔη+eηj+1−eηj.The temporal semidiscrete form is obtained with the help of equations (32)–(34):(33)Δη−ϑΓ2−ϑ∑m=1nςm+1−ςmn+32−m1−ϑ+n+12−m1−ϑ+ςn+1−ςn21−ϑ=bςnζe−ηn+e−ηn+12Δη∑j=0n−1ςj+1eηj+1Δη−eηj+1+eηj+ςj−eηjΔη+eηj+1−eηj+bςnζe−ηn+12Δηςn+1eηn+1Δη−eηn+1+eηn+ςn−eηnΔη+eηn+1−eηn+∂2ςnζ∂ζ2+aςnζ1−ςnζ.


## 3.2. Use of Quasi-Wavelet-Based Approximation for the Discretization of Spatial Derivatives

By using the approximation by quasi-wavelets given in Section2, we discretize the unknown function and its spatial derivatives. The spatial step is chosen as Δζ=1/N. We denote the value of the unknown function ςζ,η at spatial point ζi and time η=ti by ςin with i=0,1,…,N and n=0,1,…,M. By using the approximation from equation (21), we have(34)ςnxi=∑p=i−Wi+Wςζk,ηnδΔ,σnζi−ζk,n=0,1,2,…,i=0,1,…,N−1.With the help of equation (33), we can rewrite equation (34) as follows:(35)Δη−ϑΓ2−ϑ∑m=1nςm+1−ςmn+32−m1−ϑ+n+12−m1−ϑ+ςn+1−ςn21−ϑ=b∑k=i−Wi+Wςζk,ηnδΔ,σζi−ζke−tηn+e−ηn+12Δη∑j=0n−1ςj+1eηj+1Δη−eηj+1+eηj+ςj−eηjΔη+eηj+1−eηj+b∑p=i−Wi+Wςζk,ηnδΔ,σζi−ζke−ηn+12Δηςn+1eηn+1Δη−eηn+1+eηn+ςn−eηnΔη+eηn+1−eηn+∑p=i−Wi+Wςζk,ηnδΔ,σ2ζi−ζk+a∑p=i−Wi+Wςζk,ηnδΔ,σζi−ζk1−∑p=i−Wi+Wςζk,ηnδΔ,σζi−ζk.Consideringxp−xi=xk, we have the following:(36)Δη−ϑΓ2−ϑ∑m=1nςim+1−ςimn+32−m1−ϑ+n+12−m1−ϑ+ςin+1−ςin21−ϑ=b∑k=−W+Wuk+inδΔ,σ−kΔζ×e−ηn+e−ηn+12Δη∑j=0n−1ςij+1eηj+1Δη−eηj+1+eηj+ςij−eηjΔη+eηj+1−eηj+b∑k=−W+Wuk+inδΔ,σ−kΔζ×e−ηn+12Δηςin+1eηn+1Δη−eηn+1+eηn+ςin−eηnΔη+eηn+1−eηn+∑k=−W+Wuk+inδΔ,σ2ζi−ζk+a∑k=−Wi+Wuk+inδΔ,σ−kΔζ1−∑k=−W+Wuk+inδΔ,σ−kΔζ.This is the full discrete form of the taken model in case of the kernelκζ,η=e−η. For other kernels κζ,η=1 and κζ,η=e−ζ2t, the discrete form is obtained as follows:(37)Δη−ϑΓ2−ϑ∑m=1nςim+1−ςimn+32−m1−ϑ+n+12−m1−ϑ+ςin+1−ςin21−ϑ=bΔη4∑j=0n−1uij+ςij+1+Δη4∑j=0nuij+ςij+1+∑k=−W+Wuk+inδΔ,σ2ζi−ζk+a∑k=−Wi+Wuk+inδΔ,σ−kΔζ1−∑k=−W+Wuk+inδΔ,σ−kΔζ,Δη−ϑΓ2−ϑ∑m=1nςim+1−ςimn+32−m1−ϑ+n+12−m1−ϑ+ςin+1−ςin21−ϑ=b∑k=−W+Wuk+inδΔ,σ−kΔζ×e−ηn+1xi22Δη∑j=0n−1ςij+1eηj+1xi2xi2Δη−eηj+1xi2xi4+eηjxi2xi4+uij−eηj+1xi2xi2Δη+eηj+1xi2xi4−eηjxi2xi4+∑k=−W+Wuk+inδΔ,σ2ζi−ζk+a∑k=−Wi+Wuk+inδΔ,σ−kΔζ1−∑k=−W+Wuk+inδΔ,σ−kΔζ.The discretization of initial and boundary conditions is performed by adopting a technique because the value ofςζk is not defined outside the domain 0,1:(38)ςi,n=ςin=0,i>M,i<0,n=0,1,…,N,ui0=f3ζi,i=0,1,…,M.


## 3.3. Theorem (see [21])

Considering the functionfζ which is a member of the space L∞∩L2Ω∩CsΩ and band limited to B, s∈Z+, σ=Δr, W∈N, W≥rs/2. Then, we have the following required result:(39)fs−∑k=−WWδσ,Δsζ−ζkfζk≤β×exp−γ22r2,where(40)γ=minr2π−BΔ,W,β=2BfLsΩ+2rfL∞Ω×eπs+1!rγπΔs.


## 4. Numerical Results and Validation

In this section, we take some numerical examples of the reaction-diffusion equation and the integro reaction-diffusion equation. We show the validity and feasibility by showing the absolute error tables. The formula of finding errorsEr between exact and numerical solutions is as follows:(41)Er=ςζi,ηj−ς˜ζi,ηj,where ς˜ζi,ηj denotes the numerical solution at point ζi,ηj. We do all numerical computation in Wolfram Mathematica version-11.3.Example 1.
We take the following particular case of our model with parametersa=1, ϑ=0.9, and b=0:(42)∂ϑςη,ζ∂ηϑ=∂2ςη,ζ∂ζ2+ςη,ζ1−ςη,ζ+fζ,η.
The initial and boundary conditions are taken as follows:(43)ςζ,0=sin2πζ,ς0,η=0,ς1,η=0.
We determine the force functionfζ,η by putting the value of the exact solution ςζ,η=sin2πζ/η+1 in the taken particular case of the model. The graph is plotted for the exact and approximate solutions with parameters’ value M=30,Δη=0.0001,andN=10 which can be seen in Figure 1. The variation of absolute errors between exact and numerical solutions is depicted in Table 1.Figure 1
Graph of exact and numerical solutions with parametric valueςζ,η for M=30, ϑ=0.9, W=30, Δη=0.00001, and r=3.2.Table 1
Representation of absolute errors for differentMandΔη.
M↓Δη=0.0001Δη=0.000001106.50×10−23.6×10−3154.5×10−21.1×10−3207.60×10−31.15×10−4309.10×10−42.2×10−5Example 2.
Taking the following integro diffusion equation as a particular case of our model with parametric valueκζ,η=1,a=0, ϑ=1, and b=1:(44)∂ςη,ζ∂η=∂2ςη,ζ∂ζ2−∫0ηςs,ζds+fζ,η.
The boundary and initial conditions are taken as follows:(45)ς0,ζ=121−ζ2,ςη,0=coshηsinh2η+2,ςη,1=0.
The exact solution of the above problem is taken asςζ,η=1−ζ2coshη/sinh2η+2 while force function fζ,η is chosen with the help of exact solutions and the given problem. The graph is plotted for the numerical and exact solutions with parameters’ value N=10,M=30,andΔη=0.0001 which can be seen in Figure 2. The variation of absolute errors between exact and numerical solutions is depicted in Table 2.Figure 2
Graph of approximate and exact and numerical solutions with parametric valueςζ,η for M=30, ϑ=1, W=30, Δη=0.00001, and r=3.2.Table 2
Representation of absolute errors for differentMandΔη.
M↓Δη=0.0001Δη=0.000001109.40×10−43.4×10−5157.1×10−42.4×10−5204.6×10−41.9×10−5304.0×10−41.3×10−5Example 3.
If we takea=1, ϑ=1, b=1, and kernel κη,ζ=e−ζ2t, then we get(46)∂ςη,ζ∂η=∂2ςη,ζ∂ζ2+ςη,ζ1−ςη,ζ+ςη,ζ∫0ηe−ζ2η,ζςη,ζds+fη,ζ.
With the following initial and boundary conditions:(47)ςζ,0=0,ς0,η=sinη,ς1,η=sinη2.
The functionςζ,η=sinη/1+ζ2 is an exact solution of the problem. The graph is plotted for the numerical and exact solutions with parameters’ value N=10,M=30,andΔη=0.0001 which can be seen in Figure 3. The variation of absolute errors between exact and numerical solutions is depicted in Table 3.Figure 3
Graph of approximate and exact solutions with parametric valueςζ,η for ϑ=1, W=30, Δη=0.00001, M=30, and r=3.2.Table 3
Representation of absolute errors for differentMandΔη.
M↓Δη=0.0001Δη=0.000001101.97×10−77.4×10−13151.90×10−74.8×10−13201.84×10−73.5×10−13301.7×10−72.2×10−13Example 4.
If we takeκη,ζ=e−η, a=1, ϑ=0.9, and b=1, then our model (42) is reduced to(48)∂0.9ςη,ζ∂η0.9=∂2ςη,ζ∂ζ2+ςη,ζ1−ςη,ζ+ςη,ζ∫0ηe−η−sςs,ζds+fζ,η.
The following equations are chosen as initial and boundary conditions:(49)ςζ,0=1−ζ4sinζ,ς0,η=sinη,ς1,η=0,with the exact solution ςζ,η=sinη+ζ1−ζ4.
The graph is plotted for the exact and approximate solutions with parameters’ valueN=10,M=30,andΔη=0.0001 which can be seen in Figure 4. The variation of absolute errors between exact and numerical solutions is depicted in Table 4. The graphs and error tables clearly show that our numerical results are in feasible and desirable accuracy.Figure 4
Graph of approximate and exact solutions with parametric valueςζ,η for M=30, ϑ=0.9, W=30, Δη=0.00001, and r=3.2.Table 4
Representation of absolute errors for differentM and Δη.
M↓Δη=0.0001Δη=0.000001109.80×10−44.9×10−5155.5×10−43.3×10−5202.90×10−42.5×10−5301.50×10−41.7×10−5Table 5
Comparison of absolute errors for our proposed method and the method given in [22].
ζ,η↓Our method withΔη=10−7Method given in literature(0.1,1)6,9×10−82.51×10−7(0.2,1)5.9×10−81.07×10−6(0.3,1)3.5×10−82.13×10−6(0.4,1)6.5×10−83.21×10−6(0.5,1)5.7×10−84.13×10−6(0.6,1)4.9×10−84.74×10−6(0.7,1)8.2×10−74.89×10−6(0.8,1)1.5×10−84.37×10−6(0.9,1)4.40×10−72.8×10−6Example 5.
Let us take another example of the Burger–Fisher equation:(50)∂ςζ,η∂η=∂2ςζ,η∂ζ2+ς∂ςζ,η∂ζ+ς1−ς,0<ζ<1,η>0,with following initial and boundary conditions:(51)ςζ,0=12+12tanhζ4,ς0,η=12+12tanh5η8,ς1,η=12+12tanh141+5η2,where the exact solution of the above problem is 1/2+1/2tanh1/4ζ+5η/2. Now, our obtained numerical results are compared with the results given in [22].


## 5. Conclusion

In this research work, we developed a numerical scheme for time-fractional partial differential equations and integro fractional partial differential equations. The discretization of the time-fractional derivative is done by using the Crank–Nicolson scheme, and the discretization of the unknown function and its spatial derivatives is performed with the help of the quasi-wavelet-based numerical method. For different kernels, we show the validity and effectiveness of our method. Our method has a good accuracy for both fractional PDEs and fractional integro PDEs as seen from the error tables. The same can be seen by the plotted graph for the existing exact solution and numerical solution.



---

*Source: 3291723-2020-09-26.xml*