    # A second order difference scheme for time fractional diffusion equation with generalized memory kernel

In the current work we build a difference analog of the Caputo fractional derivative with generalized memory kernel (_λL2-1_σ formula). The fundamental features of this difference operator are studied and on its ground some difference schemes generating approximations of the second order in time for the generalized time-fractional diffusion equation with variable coefficients are worked out. We have proved stability and convergence of the given schemes in the grid L_2 - norm with the rate equal to the order of the approximation error. The achieved results are supported by the numerical computations performed for some test problems.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Differential equations with fractional order derivatives represent a powerful mathematical tool for exact and realistic description of physical and chemical processes for which it is needed to take into consideration the background (memory) of the process [1, 2, 3, 4]. The patterns which take memory into consideration in such equations are the memory functions that are the kernels of integrals defining the operators of fractional integro-differentiation. For fractional integro-differentiation operators, the memory functions are namely power functions. The exponent of the power function of memory defines the order of the derivative and is connected with the fractal dimension of the environment in which the described process takes place. For more accurate description of the process in heterogeneous porous media, differential equations with fractional derivatives of distributed order are often used too. Processes of the memory can be described with the help of the memory function which has more complex structure than the power function.

In the rectangle we consider the Dirichlet boundary value problem for time fractional diffusion equation with generalized memory kernel and variable coefficients

 ∂α,λ(t)0tu=Lu+f(x,t),0
 u(0,t)=0,u(1,t)=0,0≤t≤T,u(x,0)=u0(x),0≤x≤1, (2)

where

 Lu=∂∂x(k(x,t)∂u∂x)−q(x,t)u,
 ∂α,λ(t)0tu(x,t)=1Γ(1−α)t∫0λ(t−η)(t−η)α∂u∂η(x,η)dη

is the generalized Caputo fractional derivative of order , with weighting function , where , for all ; , for all .

Diffusion and Fokker-Planck-Smoluchowski equations which have a generalized memory kernel were investigated in . In this work it is demonstrated that the memory kernel appearing in the generalized diffusion equation has diverse potential forms which can describe a broad range of experimental phenomena.

With the help of the energy inequality method, a priori estimates for the solution of both differential and difference problems of the Dirichlet and Robin boundary value problems for the fractional, variable and distributed order diffusion equation with Caputo fractional derivative were derived in [6, 7, 8, 9, 10, 11, 12, 13]. A priori estimates for the difference problems analyzed in  by means of the maximum principle imply the stability and convergence of these difference schemes.

In this work, to construct difference schemes with the order of accuracy in time we have to demand the existence of a sufficiently smooth solution of the original problem. It brings on a significant narrowing of the input data class of the problem for which we apply the proposed method. As it is well known (see for example [15, 16]), in the case of smooth input data for a time-fractional diffusion equation, the solutions are not necessarily smooth in a closed domain, because the derivatives of the function with respect to might possess a singularity at . In such cases, if possible, we present the solution as the sum of two functions: one of which is known but not smooth, whereas the other is smooth but not known, as it is illustrated in work .

In work , we consider a reaction-diffusion problem with a Caputo time derivative of the order . It is shown that the solution of such a problem has in general a weak singularity near the initial time , and we derive sharp pointwise bounds on certain derivatives of this solution. We have given a new analysis of a standard finite difference method for the problem, taking into account this initial singularity.

In , we study an analysis of the L1 scheme for the subdiffusion equation with nonsmooth data. In , error estimates for approximations of distributed order time fractional diffusion equation with nonsmooth data were investigated.

In the current paper, a difference analog of the Caputo fractional derivative with generalized memory kernel (L2-1 formula) is built up. The essential features of this difference operator are investigated and on its ground some difference schemes generating approximations of the second and fourth order in space and the second order in time for the generalized time-fractional diffusion equation with variable coefficients are studied. Stability of the suggested schemes as well as their convergence in the grid - norm with the rate equal to the order of the approximation error are proven. The achieved results are supported by the numerical computations performed for some test problems.

## 2 Stability and convergence of the family of difference schemes

In this section, we consider some families of difference schemes in a general form, set on a non-uniform time grid. A criterion of the stability of the difference schemes in the grid - norm is worked out. The convergence of solutions of the difference schemes to the solution of the corresponding differential problem with the rate equal to the order of the approximation error is proven.

In the rectangle we assign the grid , where , .

The family of difference schemes, approximating problem (1)–(2) on the grid , mainly has the form

 gΔα0tj+1yi=Λy(σj+1)i+φj+1i,i=1,2,…,N−1,j=0,1,…,M−1, (3)
 y(0,t)=0,y(l,t)=0,t∈¯¯¯ωτ,y(x,0)=u0(x),x∈¯¯¯ωh, (4)

where

 gΔα0tj+1yi=j∑s=0(ys+1i−ysi)gj+1s,gj+1s>0, (5)

is a difference analog of the generalized Caputo fractional derivative of the order with weighting function (, , ), is a difference operator which approximates the continuous operator , such that the operator preserves its positive definiteness:

 (−Λy,y)≥ϰ∥y∥20,(y,v)=N−1∑i=1yivih,∥y∥20=(y,y),ϰ>0,

, , at .

###### Lemma 2.1

 If , then for any function defined on the grid the following inequalities hold true

 vj+1gΔα0tj+1v≥12gΔα0tj+1(v2)+12gj+1j(gΔα0tj+1v)2, (6)
 vjgΔα0tj+1v≥12gΔα0tj+1(v2)−12(gj+1j−gj+1j−1)(gΔα0tj+1v)2, (7)

where .

###### Corollary 2.1

 If and , where , , then for any function defined on the grid we have the inequality

 (σj+1vj+1+(1−σj+1)vj)gΔα0tj+1v≥12gΔα0tj+1(v2). (8)
###### Theorem 2.1

 If

 gj+1j>gj+1j−1>…>gj+10≥c2>0,gj+1j2gj+1j−gj+1j−1≤σj+1≤1,

where , , then the difference scheme (3)–(4) is unconditionally stable and its solution satisfies the following a priori estimate:

 ∥yj+1∥20≤∥y0∥20+12ϰc2max0≤j≤M∥φj∥20, (9)

A priori estimate (9) implies the stability of difference scheme (3)–(4).

###### Theorem 2.2

 If the conditions of Theorem (2.1) are fulfilled and difference scheme (3)–(4) has the approximation order , where and are some known positive numbers, then the solution of difference scheme (3)–(4) converges to the solution of differential problem (1)–(2) in the grid - norm with the rate equal to the order of the approximation error .

## 3 A second order numerical differentiation formula for the generalized Caputo fractional derivative

In this section, we construct a difference analog of the Caputo fractional derivative with the approximation order and investigate its essential properties.

Next we consider the uniform grid . Let us find the discrete analog of the at the fixed point , , where , . For all and (, ) the following equalities hold true

 ∂α,λ(t)0tj+σv(t)=1Γ(1−α)tj+σ∫0λ(tj+σ−η)(tj+σ−η)αv′(η)dη
 =1Γ(1−α)⎛⎜ ⎜⎝j∑s=1ts∫ts−1λ(tj+σ−η)(tj+σ−η)αv′(η)dη+tj+σ∫tjλ(tj+σ−η)(tj+σ−η)αv′(η)dη⎞⎟ ⎟⎠
 =1Γ(1−α)j∑s=1ts∫ts−1λ(tj+σ−η)(tj+σ−η)α(Π2,sv(η))′dη
 +1Γ(1−α)j∑s=1ts∫ts−1λ(tj+σ−η)(tj+σ−η)α(v(η)−Π2,sv(η))′dη
 +1Γ(1−α)tj+σ∫tjλ(tj+σ−η)(tj+σ−η)α(Π1,jv(η))′dη
 +1Γ(1−α)tj+σ∫tjλ(tj+σ−η)(tj+σ−η)α(v(η)−Π1,jv(η))′dη
 =1Γ(1−α)j∑s=1vt,s−1ts∫ts−1λ(tj+σ−η)(tj+σ−η)αdη
 +1Γ(1−α)j∑s=1v¯tt,sts∫ts−1λ(tj+σ−η)(η−ts−1/2)(tj+σ−η)αdη+
 +vt,jΓ(1−α)tj+σ∫tjλ(tj+σ−η)(tj+σ−η)αdη+R(1)1j+R(1)jj+σ
 =1Γ(1−α)j∑s=1⎛⎜⎝vt,s−1ts∫ts−1λj−s+σ+1/2−λt,j−s+σ(η−ts−1/2)(tj+σ−η)αdη
 +λj−s+σv¯tt,sts∫ts−1(η−ts−1/2)(tj+σ−η)αdη⎞⎟⎠
 +λσ−1/2vt,jΓ(1−α)tj+σ∫tjdη(tj+σ−η)α+R(1)1j+R(1)jj+σ+R(2)1j+R(2)jj+σ+R(3)1j
 =τ1−αΓ(2−α)j∑s=1(vt,s−1(λj−s+σ+1/2a(α)j−s+1+(λj−s+σ−λj−s+σ+1)b(α)j−s+1)
 +λj−s+σb(α)j−s+1(vt,s−vt,s−1))+τ1−αΓ(2−α)λσ−1/2a(α)0vt,j+Rj+σ1
 =τ1−αΓ(2−α)((λj+σ−1/2a(α)j−λj+σb(α)j)vt,0
 +j−1∑s=1(λj−s+σ−1/2a(α)j−s+λj−s+σb(α)j−s+1−λj−s+σb(α)j−s)vt,s
 +(λσ−1/2a(α)0+λσb(α)1)vt,j)
 =τ1−αΓ(2−α)j∑s=0c(α)j−svt,s+Rj+σ1.

where

 a(α)0=σ1−α,a(α)l=(l+σ)1−α−(l−1+σ)1−α,
 b(α)l=12−α[(l+σ)2−α−(l−1+σ)2−α]−12[(l+σ)1−α+(l−1+σ)1−α],l≥1,
 λs=λ(ts),vt,s=v(ts+1)−v(ts)τ,v¯tt,s=v(ts)−v(ts−1)τ,
 Π1,sv(t)=v(ts+1)t−tsτ+v(ts)ts+1−tτ,
 Π2,sv(t)=v(ts+1)(t−ts−1)(t−ts)2τ2
 −v(ts)(t−ts−1)(t−ts+1)τ2+v(ts−1)(t−ts)(t−ts+1)2τ2,
 Rj+σ1=R(1)1j+R(1)jj+σ+R(2)1j+R(2)jj+σ+R(3)1j,
 R(1)1j=1Γ(1−α)j∑s=1ts∫ts−1λ(tj+σ−η)(tj+σ−η)α(v(η)−Π2,sv(η))′dη,
 R(1)jj+σ=1Γ(1−α)tj+σ∫tjλ(tj+σ−η)(tj+σ−η)α(v(η)−Π1,jv(η))′dη,
 R(2)1j=1Γ(1−α)j∑s=1vt,s−1ts∫ts−1λ(tj+σ−η)−λj−s+σ+1/2+λt,j−s+σ(η−ts−1/2)(tj+σ−η)αdη,
 R(2)jj+σ=vt,jΓ(1−α)tj+σ∫tjλ(tj+σ−η)−λσ−1/2(tj+σ−η)αdη,
 R(3)1j=1Γ(1−α)j∑s=1v¯tt,sts∫ts−1(λ(tj+σ−η)−λj−s+σ)(η−ts−1/2)(tj+σ−η)αdη.

Let us consider the below fractional numerical differentiation formula for the generalized Caputo fractional derivative of order with weighting function ()

 Δα,λ(t)0tj+σv=τ1−αΓ(2−α)j∑s=0c(α)j−svt,s, (10)

where

 c(α)0=λσ−1/2a(α)0,forj=0;and forj≥1,
 c(α)s=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩λσ−1/2a(α)0+λσb(α)1,s=0,λs+σ−1/2a(α)s+λs+σb(α)s+1−λs+σb(α)s,1≤s≤j−1,λj+σ−1/2a(α)j−λj+σb(α)j,s=j. (11)

We call (10) the L2-1 - formula for the generalized Caputo fractional derivative.

###### Lemma 3.1

For any and , it is true that

 ∂α,λ(t)0tj+σv=Δα,λ(t)0tj+σv+O(τ2), (12)

where , and .

Proof. We have .

Estimate the errors , , , and :

 |R(1)1j|=1Γ(1−α)∣∣ ∣ ∣∣j∑s=1ts∫ts−1λ(tj+σ−η)(tj+σ−η)α(v(η)−Π2,sv(η))′dη∣∣ ∣ ∣∣
 ≤1Γ(1−α)j∑s=1∣∣ ∣ ∣∣ts∫ts−1(−λ′(tj+σ−η)(tj+σ−η)α+αλ(tj+σ−η)(tj+σ−η)α+1)(v(η)−Π2,sv(η))dη∣∣ ∣ ∣∣
 ≤Mj+13τ39√3Γ(1−α)j∑s=1ts∫ts−1(mj+11(tj+σ−η)α+αλ(0)(tj+σ−η)α+1)dη
 =Mj+13τ39√3Γ(1−α)tj∫0(mj+11(tj+σ−η)α+αλ(0)(tj+σ−η)α+1)dη
 ≤Mj+13τ39√3Γ(1−α)⎛⎝mj+11t1−αj+σ1−α+λ(0)σατα⎞⎠=O(τ3−α),
 |R(1)jj+σ|=1Γ(1−α)∣∣ ∣ ∣∣tj+σ∫tjλ(tj+σ−η)(tj+σ−η)α(v(η)−Π1,jv(η))′dη∣∣ ∣ ∣∣
 =1Γ(1−α)∣∣ ∣ ∣∣tj+σ∫tjλ(tj+σ−η)(tj+σ−η)α(v′(η)−vt,j)dη∣∣ ∣ ∣∣
 =∣∣ ∣ ∣∣v′′(tj+1)Γ(1−α)tj+σ∫tjλ(tj+σ−η)(η−tj+1/2)(tj+σ−η)αdη+O(τ3−α)∣∣ ∣ ∣∣
 =∣∣ ∣ ∣∣v′′(tj+1)λ(tσ−1/2)Γ(1−α)tj+σ∫tj(η−tj+1/2)(tj+σ−η)αdη+O(τ3−α)∣∣ ∣ ∣∣
 =∣∣ ∣∣v′′(tj+1)λ(tσ−1/2)σ1−ατ2−αΓ(3−α)(σ−1+α/2)+O(τ3−α)∣∣ ∣∣=O(τ3−α),
 |R(2)1j|=1Γ(1−α)∣∣ ∣ ∣∣j∑s=1vt,s−1ts∫ts−1λ(tj+σ−η)−λj−s+σ+1/2+λt,j−s+σ(η−ts−1/2)(tj+σ−η)αdη∣∣ ∣ ∣∣
 ≤Mj+11mj+12τ24Γ(1−α)j∑s=1ts∫ts−1dη(tj+σ−η)α=Mj+11mj+12τ24Γ(1−α)tj∫0dη(tj+σ−η)α
 ≤Mj+11mj+12t1−αj+στ24Γ(1−α)=O(τ2),
 |R(2)jj+σ|=∣∣ ∣ ∣∣vt,jΓ(1−α)tj+σ∫tjλ(tj+σ−η)−λσ−1/2(tj+σ−η)αdη∣∣ ∣ ∣∣
 =∣∣ ∣ ∣∣vt,jΓ(1−α)tj+σ∫tj−λ′(tσ−1/2)(η−tj+1/2)+12λ′′(¯ξ)(η−tj+1/2)2(tj+σ−η)αdη∣∣ ∣ ∣∣
 ≤Mj+11mj+12σ1−α4Γ(2−α)τ3−α=O(τ3−α),
 |R(3)1j|=1Γ(1−α)∣∣ ∣ ∣∣j∑s=1v¯tt,sts∫ts−1(λ(tj+σ−η)−λj−s+σ)(η−ts−1/2)(tj+σ−η)αdη∣∣ ∣ ∣∣
 1Γ(1−α)∣∣ ∣ ∣∣j∑s=1v¯tt,sts∫ts−1−λ′(¯ξ2)(η−ts)(η−ts−1/2)(tj+σ−η)αdη∣∣ ∣ ∣∣
 ≤Mj+12mj+11τ22Γ(1−α)tj∫0dη(tj+σ−η)α=Mj+12mj+11t1−αj+στ