# An Efficient Hybrid Numerical Scheme for Nonlinear Multiterm Caputo Time and Riesz Space Fractional-Order Diffusion Equations with Delay

Authors: A. K. Omran; M. A. Zaky; A. S. Hendy; V. G. Pimenov
Journal: Journal of Function Spaces
 (2021)

---

## Abstract

In this paper, we construct and analyze a linearized finite difference/Galerkin–Legendre spectral scheme for the nonlinear multiterm Caputo time fractional-order reaction-diffusion equation with time delay and Riesz space fractional derivatives. The temporal fractional orders in the considered model are taken as0<β0<β1<β2<⋯<βm<1. The problem is first approximated by the L1 difference method on the temporal direction, and then, the Galerkin–Legendre spectral method is applied on the spatial discretization. Armed by an appropriate form of discrete fractional Grönwall inequalities, the stability and convergence of the fully discrete scheme are investigated by discrete energy estimates. We show that the proposed method is stable and has a convergent order of 2−βm in time and an exponential rate of convergence in space. We finally provide some numerical experiments to show the efficacy of the theoretical results.


---

## Body

## 1. Introduction

Fractional-order partial differential equations have evolved into powerful tools for describing a wide range of anomalous behavior and complex systems in natural science and engineering [1–8]. In addition, time delay occurs frequently in realistic world and it has been considered in numerous mathematical models, e.g., automatic control systems with feedback and population dynamics. Moreover, fractional differential equations with delay have been used widely in a variety of scientific and technical disciplines, including the study of natural phenomena, mathematical modelling, and the studies of porous media [9, 10]. Recently, a two-term time-fractional differential equation that contains specific instances of the fractional diffusion-wave problem (see, for example, [11, 12]) has been investigated in the literature.In recent years, multiterm time-fractional differential equations have attracted the attention of many researchers. The ability of these equations to describe complex multirate physical phenomena is a motivating force behind their development (see, e.g., [13–15]). They were proposed to improve the modelling accuracy by accurately depicting the anomalous diffusion process [16], accurately modelling various types of viscoelastic damping [17], and accurately reproducing the unsteady flow of a fractional Maxwell fluid [18] and Oldroyd-B fluid [19]. Daftardar-Gejji and Bhalekar [20] considered the multiterm time-fractional diffusion wave equation with constant coefficients. Through the use of a domain decomposition technique, they were able to derive the linear and nonlinear diffusion-wave equations of the fractional order. Luchko [21] used an appropriate maximum principle and the Fourier technique to study the existence, uniqueness, and a priori estimates for the multiterm time-fractional diffusion equation with variable coefficients. A new analytic technique for solving three types of multiterm time-space fractional advection diffusion equations with nonhomogeneous Dirichlet boundary conditions was proposed by Jiang et al. [22], based on Luchko’s theorem and the equivalent relationship between the Laplacian operator and the Riesz fractional derivative. Ding and Jiang [23] used the technique of spectral representation of the fractional Laplacian operator in order to provide the analytical solutions for the multiterm time-space fractional advection-diffusion equations with mixed boundary conditions. For the solution of initial-boundary value problems of multiterm time fractional diffusion equations, Li et al. [24] examined the well-posedness and long-time asymptotic behaviour of the equations. Zaky [25] constructed a Legendre spectral tau algorithm to deal with the multiterm time-fractional diffusion equations. Hendy [26] presented numerical treatment for solving a class of one-dimensional multiterm time-space fractional advection-diffusion equations with a temporal delay of the functional type. Hendy and De Staelen [27] developed a high-order numerical approximation approach for multiterm time convection diffusion wave equations with a nonlinear fixed time delay. To solve a coupled system of nonlinear multiterm time-space fractional diffusion equations over a nonuniform temporal mesh, Hendy and Zaky [28] developed an effective finite difference/spectral approach. Very recently, Zaky et al. [29] presented a discrete fractional Grönwall inequality that is consistent with the L2−1σ to cope with the analysis of multiterm time-fractional partial differential equations. The key advantage of the proposed discrete Grönwall inequality over earlier efforts was that it can be utilised to provide optimal error estimates for multiterm fractional problems with nonlinear delay. Inspired by these inequalities, we can state and prove the convergence and stability estimates for our proposed fully discrete scheme. The discrete versions of Grönwall inequalities are of high concern in the numerical analysis of the numerical schemes for fractional differential equations [30, 31].Single-term fractional differential equations are often unable to describe some of the changing characteristics of the systems accurately. However, several multiterm fractional differential equations provide us with new tools to solve such problems. The multiterm time-fractional diffusion equation is useful not only for modelling the behaviour of viscoelastic fluids and rheological materials [32] but also for approximating distributed-order differential equations [33]. Hence, studies on the multiterm time-fractional differential equations have become important and useful for different applications. The multiterm time-fractional diffusion equation, whose weight function is taken into the linear combination of the Dirac δ-functions, is an important special case of the time-fractional diffusion equation of distributed order. In this paper, we consider the numerical approximations to the following generalized nonlinear multiterm time-space fractional reaction-diffusion equations with delay:
(1)∑r=0mqr∂βru∂tβr=κ∂αu∂xα+fux,t,ux,t−s~+gx,t,x∈Ω,t∈I,endowed with initial-boundary conditions of the form
(2)ux,t=ψx,t,x∈Ω,t∈−s~,0,ua,t=ub,t=0,t∈I.Here,Ω=a,b⊂ℝ and I=0,T⊂ℝ are space and time domains, respectively. We denote ∂βr/∂tβr as the Caputo fractional derivative with the fractional orders 0<β0<β1<β2<⋯<βm<1. The parameters κ,s are positive constants. The parameters qr are absolutely positive. Also, 1<α<2 is the space fractional order. The left and right Riemann-Liouville fractional derivatives of order αn−1<α<n on the infinite domain [34] are defined as
(3)−∞Dxαux,t=1Γn−α∂n∂xn∫−∞xx−τn−1−αuτ,tdτ,xD∞αux,t=−1nΓn−α∂n∂xn∫x∞τ−xn−1−αuτ,tdτ,where Γx is the gamma function. Thus, the space fractional derivative in the Riesz form on the space interval Ω can be defined as [35]
(4)∂αu∂xα=−cαDaxαux,t+xDbαux,t,cα=12cosπα/2,1<α<2.The Caputo derivative∂β/∂tβ is defined as
(5)∂β∂tβux,t=1Γ1−β∫0tt−r−β∂∂rux,rdr,0<β<1.The main aim of this work is to construct and analyze an efficient linearized numerical scheme for the nonlinear multiterm Riesz space and Caputo time fractional reaction-diffusion problem with fixed delay. A hybrid numerical scheme combines the Galerkin–Legendre spectral schemes, and a uniformL1-type interpolation technique is designed. The main challenges of the considered work are represented in how to numerically approximate the time Caputo fractional derivative, Riesz space fractional derivatives, and the time delay to produce an easy-to-implement and consistent numerical scheme. Overcoming all of these challenges to yield a hybrid linear numerical scheme is a first target. The other target is to analyze the convergence and stability. The theoretical analysis of the constructed fully discrete scheme is successfully estimated using appropriate discrete fractional Grönwall inequalities, and the scheme is proven to be unconditionally stable and convergent.The outline of this paper is as follows. In the next section, we will go over some essential definitions and properties of fractional derivative spaces, fractional Sobolev spaces, and Jacobi polynomials. The steps needed to construct a fully discrete scheme on a uniform mesh are detailed in Section3. Some technical lemmas from the literature are summarized in Section 4. Furthermore, the stability and the convergence analyses of the fully discrete scheme are studied in Section 5. Finally, numerical experiments are performed in Section 6 to illustrate the convergence analysis of the proposed approach.


## 2. Preliminaries

We here give some essential fractional derivative spaces [36] and their required features which will be helpful in the coming analysis. After that, the definition of Jacobi polynomials and their basic properties are recalled. We now fix some notations for the sake of clearness:
(i)
·,·0,Ω is the inner product on the space L2Ω with the L2-norm ·0,Ω(ii)
The maximum norm is defined as·∞(iii)
C0∞Ω is the space of smooth functions with compact support in Ω(iv)
HrΩ and H0rΩ are the usual Sobolev spaces with the norm ·r and seminorm ·r(v)
PNΩ is the space of polynomials defined on the domain Ω with degree at most N(vi)
The approximation spaceWN0 is defined as WN0=PNΩ∩H01Ω(vii)
IN is the Legendre-Gauss-Lobatto interpolation operator IN:CΩ¯⟶WN as
(6)uxk=INuxk∈PN,k=0,1,⋯,N.Definition 1.
Fractional derivative spaces and their related norms and seminorms are defined as follows [36]:
(i)
Left fractional space: letη>0. Define the seminorm uJLηΩ=aDxηu0,Ω and the norm uJLηΩ=uJLη2Ω+u0,Ω21/2, and let JLη (or JL,0η) denote the closure of C∞Ω (or C0∞Ω) with respect to ·JLηΩ(ii)
Right fractional space: letη>0. Define the seminorm uJRηΩ=xDbηu0,Ω and the norm uJRηΩ=uJRηΩ2+u0,Ω21/2, and let JRη (or JR,0η) denote the closure of C∞Ω (or C0∞Ω) with respect to ·JRηΩ(iii)
Symmetric fractional space: letη≠n−1/2,n∈ℕ. Define the seminorm uJsηΩ=Daxηu,Dxbηu0,Ω1/2 and the norm uJsηΩ=uJsη2Ω+u0,Ω21/2, and let Jsη (or Js,0η) denote the closure of C∞Ω (or C0∞Ω) with respect to ·JsηΩ(iv)
Fractional Sobolev space: letη>0. Define the fractional Sobolev space HηΩ as HηΩ=u∈L2ΩωηFu^∈L2ℝ, endowed with the seminorm uHηΩ=ωηFu^0,ℝ and the norm uHηΩ=uHηΩ2+u0,Ω21/2, where Fu^ is the Fourier transformation of u^ and u^ is the extension of zero of u outside Ω. Denote by HηΩ (or H0ηΩ) the closure of C∞Ω (or C0∞Ω) with respect to ·HηΩ.Lemma 2 (see [36]).
The spacesJLη,JRη,Jsη, and Hη are equivalent, with equivalent seminorms and norms if η≠n−1/2,n∈ℕ.Lemma 3 (adjoint property).
Let1<η<2, then for any u∈H0ηΩ and υ∈H0η/2Ω, we get
(7)Daxηu,υ0,Ω=Daxη/2u,Dxbη/2υ0,Ω,Dxxηu,υ0,Ω=Dxxη/2u,Dabη/2υ0,Ω.Spectral methods are characterized by the expansion of the solution in terms of global and, usually, orthogonal polynomials [37–40]. Now, we present the Jacobi polynomials and some of their basic properties. The vital role in the field of spectral methods arose from the nature of Jacobi weights which are related to the singular kernels of time Caputo fractional derivatives of order 0<β<1. Denote Jiμ,νx,μ,ν>−1 as the i-th order Jacobi polynomial of index defined on −1,1. As all classic orthogonal polynomials, Jiμ,νxi=0N satisfies the following three-term-recurrence relation:
(8)J0μ,νx=1,J1μ,νx=122+μ+νx+12μ−ν,Ji+1μ,νx=Aiμ,νx−Biμ,νJiμ,νx−Ciμ,νJi−1μ,νx,if1≤i≤N.The recursion coefficients are given by(9)Aiμ,ν=2i+μ+ν+12i+μ+ν+22i+1i+μ+ν+1,Biμ,ν=2i+μ+ν+1ν2−μ22i+1i+μ+ν+12i+μ+ν,Ciμ,ν=2i+μ+ν+2i+μi+νi+1i+μ+ν+12i+μ+ν.Letωμ,νx=1−xμ1+xν. Then, one has
(10)∫−11Jiμ,νxJjμ,νxωμ,νtdx=γiμ,νδi,j,∀i=0,1,⋯,N,where δi,j is the Kronecker delta function and
(11)γiμ,ν=2μ+ν+1Γi+μ+1Γi+ν+12i+μ+ν+1i!Γi+μ+ν+1,∀i=0,1,⋯,N.In particular, the Legendre polynomial is defined asLix=Ji0,0x.


## 3. The Numerical Scheme

Here, we provide a fully discrete scheme for the problems (1) and (2) based on the L1-type approximation for the Caputo time-fractional derivative and the Legendre-Galerkin spectral method in space. To discretize the time-fractional derivatives, we divide the interval 0,T uniformly with a time step size τ defined by τ=s~/Ns~ such that Ns~ is a positive integer. The uniform partitions given by tn=nτ,∀−Ns~≤n≤M, where M=T/τ. The L1 interpolation scheme for the time-fractional derivative of order 0<βr<1, r=0,1,2,⋯,m, in the Caputo sense at the time tn is defined as
(12)∑r=0mqr∂βru∂tβrt=tn=∑r=0mqr∫0tnu′x,ηω1−βrtn−ηdη=∑r=0mqrΓ1−βr∑i=1nux,ti−ux,ti−1τ∫ti−1titn−η−βrdη+rτn=∑r=0mqrΓ2−βrτβr∑i=1naiβrux,ti−ux,ti−1+rτn,where ωβrt=tβr−1/Γβr,t>0, and ajβr=j+11−βr−j1−βr, for each j≥0. If u∈C20,T;L2Ω, then there exists a constant C>0 such that the truncation error rτn satisfies rτn≤Cτ2−βm, for each n=0,1,⋯,M (see [41]).Definition 4.
Letunn=0M be a sequence of real functions defined on Ω. We define the multiterm discrete time-fractional difference operator ∑r=0mqrDτβr by
(13)∑r=0mqrDτβrun=∑r=0mqrΓ2−βrτβr∑i=1nan−iβrδtui=∑r=0mqrΓ2−βrτβr∑i=0nbn−iβrui,∀n=1,⋯,M.In this expression,δtui=ui−ui−1, and the constants are defined by b0βr=a0βr,bnβr=−an−1βr, and bn−iβr=an−iβr−an−i−1βr, for each i=1,⋯,n−1.In order to provide a semidiscretized form of (1) at each time tn, we approximate the time-fractional term through (13). Taylor approximations are used to approximate the nonlinear source function in a linear style. As a consequence, we obtain the discrete-time system:
(14)∑r=0mqrDτβrun=∂αun∂xα+f2un−1−un−2,un−Ns~+gnx,1≤n≤M,∀x∈Ω,(15)unx=ψx,−Ns~≤n≤0,x∈Ω.We define the following function space to give appropriate base functions such that the boundary conditions are satisfied exactly as clarified in spectral methods for space-fractional differential equations [42, 43].
(16)WN0=spanφnx:n=0,1,⋯,N−2,where for each x^∈−1,1, the function φn is given by
(17)φnx=Lnx^−Ln+2x^=2n+32n+11−x^2Jn1,1x^,and x=1/2b−ax^+a+b∈a,b.We introduce the parameterσr=qr/Γ2−βrτβr. Then, the scheme (14) can be rewritten in the following equivalent form:
(18)∑r=0mσra0βrun−κ∂αun∂xα=∑r=0mσran−1βru0−∑r=0mσr∑i=1n−1bn−iβrui+f2un−1−un−2,un−Ns~+gnx,∀n=1,⋯,M.The fully discreteL1-Galerkin spectral scheme consists of the set of approximations uNn∈WN0, satisfying the system:
(19)∑r=0mσra0βruNn,υ−κ∂α∂xαuNn,υ=∑r=0mσran−1βruN0,υ−∑r=0mσr∑i=1n−1bn−iβruNi,υ+INf2uNn−1−uNn−2,uNn−Ns~,υ+INgnx,υ,∀υ∈WN0,∀n=1,⋯,M,uNn=πN1,0ψtn,x,−Ns~≤n≤0,where πN1,0 is an appropriate projection operator. We expand the approximate solution as
(20)uNn=∑i=0N−2u^inφix.Substituting this expression into (19) and letting υ=φk, for each 0≤k≤N−2, we obtain the following matrix representation of the uniform L1-Galerkin spectral scheme:
(21)∑r=0mσra0βrM¯−κcαS+STUn=Kn−1+Rn−1+Fn.The notations in this expression are given by the system of identities:(22)sij=∫ΩDaxα/2φixDxbα/2φjxdx,S=siji,j=0N−2,mij=∫Ωφixφjxdx,M¯=miji,j=0N−2,hin−1=∫ΩφixINf2uNn−1−uNn−2,uNn−Nsdx,Rn−1=h0n−1,h1n−1,⋯,hN−2n−1⊤,gin=∫ΩφixINgndx,Fn=g0n,g1n,⋯,gN−2n⊤,Un=u^0n,u^1n,⋯,u^N−2n⊤,Kn−1=−∑r=0mσr∑j=0n−1bn−jβrM¯Uj.Lemma 5 (see [42, 43]).
The components of the stiffness matrixS are sij=aij−aij+2−ai+2j+ai+2j+2, for each i,j=0,1,⋯,N−2. Here,
(23)aij=∫ΩDaxα/2Lix^xDbα/2Ljx^dx=b−a21−αΓi+1Γj+1Γi−α/2+1Γj−α/2+1··∑r=0Nϖr−α/2,−α/2Jiα/2,−α/2xr−α/2,−α/2Jj−α/2,α/2xr−α/2,−α/2,and xr−α/2,−α/2,ϖr−α/2,−α/2i=0N is the set of Jacobi–Gauss points and weights with respect to the weight function ω−α/2,−α/2. The mass matrix M¯ is symmetric, with nonzero components:
(24)mij=mji=b−a2j+1+b−a2j+5,∀i=j,−b−a2j+5,∀i=j+2.


## 4. Technical Lemmas

Several lemmas that will be invoked through our analysis appeared in that section. In the sequel,C and Cu will denote generic positive constants independent of τ, N, and n and may be different under different circumstances. We also fix the following notation ℤa,b=ℤ∩a,b, such that ℤ is the set of all integers.Throughout the coming context, we will use the notation(25)Au,w=κcαaDxα/2u,xDbα/2w+xDbα/2u,aDxα/2w.The orthogonal projection operatorπNα/2,0:H0α/2Ω⟶WN0 will be such that
(26)Au−πNα/2,0u,w=0,∀u∈H0α/2Ω,w∈WN0.For convenience of theoretical analysis, we give the following seminorm and norm:(27)uα/2≔Au,u1/2,(28)uα/2≔u2+uα/221/2,which are equivalent to the seminorms and norms of JLα/2Ω,JRα/2Ω,JSα/2Ω, and Hα/2Ω. We recall the following three lemmas from [43].Lemma 6.
Letα and s be arbitrary real numbers satisfying 0<α<1,α<s,α≠1/2. Then, there exists a positive constant C independent of N such that, for any function u∈H0α/2Ω∩HsΩ, the following estimate holds:
(29)u−πNα/2,0uα/2≤CNα/2−sus.Lemma 7.
Suppose thatΩ=a,b, u∈H0α/2Ω. Then, there exist positive constants C1<1 and C2 independent of u, such that
(30)C1uα/2≤uα/2≤uα/2≤C2uHα/2Ω.The following lemma and remark summarize the properties of the interpolation operatorIN.Lemma 8 (see [44]).
Lets≥1. If u∈HsΩ, then there exists a constant C>0 independent of N, such that u−INul≤CNl−sus, for any 0≤l≤1.Remark 9.
A smooth solution of a fractional differential equation does not mean a smooth source term and vice versa. Therefore, the regularity orders of the solution u is not the same as the regularity order r of the source term g, i.e.,
(31)INg−g≤CN−rur,∀g∈HrΩ.Lemma 10 (see [45]).
For any functionut which is absolutely continuous on 0,T, the following inequality is satisfied:
(32)∂β∂tβut,ut≥12∂β∂tβut2.Lemma 11 (see [46]).
The discrete counterpart to the inequality (32) is given as
(33)Dτβuk,uk≥12Dτβuk2,such that Dτβuk is the discrete time-fractional difference operator of the L1 type as defined in (13).Plenty of researchers in recent years are stuck on the study of the continuous fractional Grönwall-type inequalities and their developments. However, the discrete fractional Grönwall-type inequality was far from well investigated, and more recently, the efforts paid in [47–50] tried to fill that gap. In what follows, we present recent discrete fractional-type inequalities. These inequalities play an important role in analyzing stability and convergence of the L1-schemes for the multiterm problems with nonlinear delay.Lemma 12 (discrete fractional Grönwall inequality [29]).
Letϕii=−Ns~∞andξll=0∞ be nonnegative sequences. Let εl, μi,∀l=0,⋯,m,i∈ℤ1,6 and c0 be positive constants independent of τ. The fractional orders are defined as 0<β0≤β1≤⋯≤βm−1≤βm≤1. If ϕi≥0∀i≥0,ϕ0 is known and ϕi=0 if i<0,
(34)∑l=0mεlDτβlϕj≤μ1ϕj+ε0ξj2∀j≤Ns~,∑l=0mεlDτβlϕj≤μ1ϕj+μ2ϕj−1+μ3ϕj−2+μ4ϕj−3+μ5ϕj−Ns~−1+μ6ϕj−Ns+ε0ξj2,∀j>Ns~,
wherec0 and μii=1,⋯,6 are positive constants. Then, there exists a positive constant τ∗≥τ such that
(35)ϕn≤2c0tnβmεmΓ1+βm∑j=1nξj2τ+ϕ0WEβm2μtnβmεm,where Eβz=∑k=0∞zk/Γ1+kβ is the Mittag-Leffler function and
(36)∑l=0mεlΔ1−βlΓ2−βl∑j=1kaj−1βl≔W>0,μ=μ1+μ2a0βm−a1βm+μ3a1βm−a2βm+μ4a2βm−a3βm+μ5aNs~−2βm−aNs~−1βm+μ6aNs~−1βm−aNs~βm.


## 5. Theoretical Analysis

The purpose of this section is to study the efficiency of the fully discrete Galerkin spectral methods for (1) and (2). We start by stability analysis and gives theorem of stability in the first subsection. The second subsection is devoted to the convergence analysis, and the theorem of convergence is given there. For the theoretical analysis requirements, we assume that the function f satisfies the following Lipschitz condition
(37)fu1,v1−fu2,v2≤Lu1−u2+v1−v2,where L is a positive constant.


### 5.1. Stability Analysis

The weak formulation of the scheme is as follows: finduNkk=1M∈PN, such that
(38)∑r=0mqrDτβruNk,υN+AuNk,υN=INf2uNk−1−uNk−2,uNk−Ns~,υN+INgk,υN,∀υN∈PN,with(39)uNk=πN1,0φk,−Ns~≤k≤0.It is a linear iterative scheme which means that we need only to get a solution to a system of linear equations at each time level. The well-posedness of that scheme is satisfied by the well-known Lax-Milgram lemma. Assume thatu~Nkk=1M is the solution of
(40)∑r=0mqrDτβru~Nk,υN+Au~Nk,υN=INf2u~Nk−1−u~Nk−2,u~Nk−Ns~,υN+INg~k,υN,∀υN∈PN,with initial conditions
(41)u~Nk=πN1,0φk,−Ns~≤k≤0.Now, we present the theorem of stability in the following context.Theorem 13.
The fully discrete scheme (38) is unconditionally stable in the sense that for all τ>0, the following holds:
(42)uNk−u~Nk2≤Cmax1≤k≤Mgk−g~k2.Proof.
DenoteηNk=uNk−u~Nk. Subtracting (40) from (38), the following holds:
(43)∑r=0mqrDτβrηNk,υN+AηNk,υN=INf2uNk−1−uNk−2,uNk−Ns~−INf2u~Nk−1−u~Nk−2,u~Nk−Ns~,υN+INgk−INg~k,υN.
According to (37) and using the Hölder inequality and Young’s inequality, we derive that
(44)INf2uNk−1−uNk−2,uNk−Ns~−INf2u~Nk−1−u~Nk−2,u~Nk−Ns~,υN≤CL2ηNk−1−ηNk−2+ηNk−Ns~υN≤ε2CL22ηNk−1−ηNk−22+ε2CL2ηNk−Ns~2+12ευN2≤4εCL2ηNk−12+εL2ηNk−22+ε2CL2ηNk−Ns~2+12ευN2,INgk−INg~k,υN≤ε2Cgk−g~k2+12ευN2.
Then, (43) becomes
(45)∑r=0mqrDτβrηNk,υN+AηNk,υN≤1ευN2+4εCL2ηNk−12+εL2ηNk−22+ε2CL2ηNk−Ns~2+εC2gk−g~k2.
TakingυN=ηNk and using Lemma 11 and (27), we can deduce that
(46)∑r=0mqr2DτβrηNk2+ηα/22≤1εηNk2+4εCL2ηNk−12+εL2ηNk−22+ε2CL2ηNk−Ns~2+εC2gk−g~k2,namely,
(47)∑r=0mqrDτβrηNk2≤2εηNk2+8εCL2ηNk−12+2εL2ηNk−22+εCL2ηNk−Ns~2+εCgk−g~k2.
By means of Lemma12 and since ε>0, there exists a positive constant τ∗=qm/2Γ2−βm2/εβm; when τ<τ∗, we have
(48)ηNk2≤2εCtkβmqmΓ1+βmEβm2μtkβm/qm∑k=1Mgk−g~k2,with μ=2/ε+8CεL2/a0βm−a1βm+2CεL2/a1βm−a2βm+CεL2/aNs−1βm−aNsβm. Thus, the scheme is unconditionally stable.☐


### 5.2. Convergence Analysis

In this section, we investigate the convergence of the fully discrete scheme (38) using error estimation.Theorem 14.
Letukk=−NsM be the exact solution of equation (1) and uNkk=−NsM be the solution of (38). Suppose that D0Ctβru∈L∞0,T;H0α/2Ω∩HsΩ,u∈L∞−s,T;HsΩ; we have
(49)uk−uNkα/2≤CNα/2−s+N−r+τ2−βm,1≤k≤M,where C is independent of N and τ.Proof.
Denoteuk−uNk=eNk=uk−πNα/2,0uk+πNα/2,0uk−uNk=Δe~Nk+e^Nk. The weak formulation of equation (1) is
(50)∑r=0mqr0CDtβruk,υN+Auk,υN=fuk,uk−Ns~,υN+gk,υN.
Subtracting (38) from (50) and owing to the definition of orthogonal projection, the error equation satisfies
(51)∑r=0mqrDτβre^Nk,υN+Ae^N,υN=ΔR1k+R2k+R3k+R4k,where
(52)R1k=INfuk,uk−Ns~−INf2uNk−1−uNk−2,uNk−Ns~,υN,R2k=fuk,uk−Ns~−INfuk,uk−Ns~,υN,R3k=∑r=0mqrDτβrπNα/2,0uk−0CDtβruk,υN,R4k=gk−INgk,υN.
We next estimate the right-hand termsR1k, R2k, R3k, and R4k. For the first term R1k,
(53)R1k=INfuk,uk−Ns~−INf2uk−1−uk−2,uk−Ns~,υN+INf2uk−1−uk−2,uk−Ns~−INf2uNk−1−uNk−2,uNk−Ns~,υN=ΔR11k+R12k.
By applying the Taylor expansion, the following holds:(54)fuk,uk−Ns~=f2uk−1−uk−2,uk−Ns+uk−2uk−1+uk−2f′1ξ,uk−Ns~=f2uk−1−uk−2,uk−Ns~+c~uτ2.Furthermore, by means of the Hölder inequality and Young’s inequality, we have
(55)R11k≤INfuk,uk−Ns~−INf2uk−1−uk−2,uk−Ns~υN≤Cfuk,uk−Ns~−f2uk−1−uk−2,uk−Ns~υN≤ε2c~uτ4+12ευN2.
According to (37), we can deduce that
(56)R12k≤LC2eNk−1−eNk−2+eNk−Ns~υN≤LC2e^Nk−1−e^Nk−2+e^Nk−Ns~+2e~Nk−1−e~Nk−2+e~Nk−Ns~υN≤8ε2CL2e^Nk−12+2ε2CL2e^Nk−22+ε2L2e^Nk−Ns~2+8ε2CL2e~Nk−12+2ε2CL2e~Nk−22+ε2CL2e~Nk−Ns~2+12ευN2.
Moreover, owing to Lemmas6 and 7, the following holds:
(57)e~Nk−12≤CC1Nα−2suk−1s2,e~Nk−22≤CC1Nα−2suk−2s2,e~Nk−Ns~2≤CC1Nα−2suk−Ns~s2.Then, (56) becomes
(58)R12k≤4εCL2e^Nk−12+εCL2e^Nk−22+ε2CL2e^Nk−Ns~2+CC1Nα−2suL∞−s,T;HsΩ2+12ευN2.
Substituting (55) and (58) into (53), we can derive that
(59)R1k≤1ευN2+4εCL2e^Nk−12+εCL2e^Nk−22+ε2CL2e^Nk−Ns~2+CC1Nα−2suL∞−s,T;HsΩ2+ε2c~uτ4.
For the second termR2k, by means of the Hölder inequality, the following holds:
(60)R2k≤ε2CN−2ruL∞−s,T;HsΩ2+12ευN2.
For the third termR3k, the following holds:
(61)R3k=∑r=0mqrDτβrπNα/2,0uk−D0CtβrπNα/2,0uk,υN+∑r=0mqrD0CtβrπNα/2,0uk−0CDtβruk,υN=πNα/2,0∑r=0mqrDτβruk−0CDtβruk,υN−∑r=0mqr0CDtβre~Nk,υN=ΔR31k+R32k.Using (12) and the Hölder inequality, the following holds:
(62)R31k≤ε2∑r=0mqrπNα/2,0Dτβruk−0CDtβruk2+12ευN2≤ε2C∑r=0mqrDτβruk−D0Ctβruk2+12ευN2≤ε2C1,uτ4−2βm+12ευN2.Furthermore, according to Lemma 2, we have
(63)R32k≤ε2CNα−2s∑r=0mqrD0Ctβruks2+12ευN2≤ε2CNα−2s∑r=0mqrD0CtβruL∞−s,T;HsΩ2+12ευN2.
Thus, (61) becomes
(64)R3k≤ε2CNα−2s∑r=0mqr0CDtβruL∞−s,T;HsΩ2+ε2C2,uτ4−2βm+1ευN2.
For the fourth termR4k, the following holds by invoking Remark 9:
(65)R4k≤ε2CNα−2ruL∞−s,T;HrΩ2+12ευN2.
Substituting (59), (60), (64), and (65) into (51), we can infer that
(66)∑r=0mqrDτβre^Nk,υN+Ae^N,υN≤52ευN2+4εCL2e^Nk−12+εCL2e^Nk−22+ε2CL2e^Nk−Ns~2+G~,where
(67)G~=εC~Nα−2s∑r=0mqr0CDtβruL∞0,T;HsΩ2+uL∞−s,T;HsΩ2+εC~N−2ruL∞−s,T;HrΩ2+εC~uτ4−2βm.
TakingυN=e^Nk in (66) and applying Lemma 11, we can conclude that
(68)∑r=0mqr2Dτβre^Nk2+e^Nkα/22≤52εe^Nk2+4εCL2e^Nk−12+εCL2e^Nk−22+ε2CL2e^Nk−Ns~2+G~,namely,
(69)∑r=0mqr2Dτβre^Nk2≤5εe^Nk2+8εCL2e^Nk−12+2εCL2e^Nk−22+εCL2e^Nk−Ns~2+G,with G=2G~. By means of Lemma 12 and since ε>0, there exists a positive constant τ∗=qm/2Γ2−βm10/εβm; when τ<τ∗, we have
(70)e^Nk2≤2GtkβmqmΓ1+βmEβm2μtkβmqm,with μ=10/ε+16εCL2/a0βm−a1βm+4CεL2/a1βm−a2βm+2CεL2/aNs~−1βm−aNs~βm. Thus, the scheme is unconditionally convergent. Finally, by means of the triangle inequality and Lemma 6, we complete the proof of (49).☐


## 5.1. Stability Analysis

The weak formulation of the scheme is as follows: finduNkk=1M∈PN, such that
(38)∑r=0mqrDτβruNk,υN+AuNk,υN=INf2uNk−1−uNk−2,uNk−Ns~,υN+INgk,υN,∀υN∈PN,with(39)uNk=πN1,0φk,−Ns~≤k≤0.It is a linear iterative scheme which means that we need only to get a solution to a system of linear equations at each time level. The well-posedness of that scheme is satisfied by the well-known Lax-Milgram lemma. Assume thatu~Nkk=1M is the solution of
(40)∑r=0mqrDτβru~Nk,υN+Au~Nk,υN=INf2u~Nk−1−u~Nk−2,u~Nk−Ns~,υN+INg~k,υN,∀υN∈PN,with initial conditions
(41)u~Nk=πN1,0φk,−Ns~≤k≤0.Now, we present the theorem of stability in the following context.Theorem 13.
The fully discrete scheme (38) is unconditionally stable in the sense that for all τ>0, the following holds:
(42)uNk−u~Nk2≤Cmax1≤k≤Mgk−g~k2.Proof.
DenoteηNk=uNk−u~Nk. Subtracting (40) from (38), the following holds:
(43)∑r=0mqrDτβrηNk,υN+AηNk,υN=INf2uNk−1−uNk−2,uNk−Ns~−INf2u~Nk−1−u~Nk−2,u~Nk−Ns~,υN+INgk−INg~k,υN.
According to (37) and using the Hölder inequality and Young’s inequality, we derive that
(44)INf2uNk−1−uNk−2,uNk−Ns~−INf2u~Nk−1−u~Nk−2,u~Nk−Ns~,υN≤CL2ηNk−1−ηNk−2+ηNk−Ns~υN≤ε2CL22ηNk−1−ηNk−22+ε2CL2ηNk−Ns~2+12ευN2≤4εCL2ηNk−12+εL2ηNk−22+ε2CL2ηNk−Ns~2+12ευN2,INgk−INg~k,υN≤ε2Cgk−g~k2+12ευN2.
Then, (43) becomes
(45)∑r=0mqrDτβrηNk,υN+AηNk,υN≤1ευN2+4εCL2ηNk−12+εL2ηNk−22+ε2CL2ηNk−Ns~2+εC2gk−g~k2.
TakingυN=ηNk and using Lemma 11 and (27), we can deduce that
(46)∑r=0mqr2DτβrηNk2+ηα/22≤1εηNk2+4εCL2ηNk−12+εL2ηNk−22+ε2CL2ηNk−Ns~2+εC2gk−g~k2,namely,
(47)∑r=0mqrDτβrηNk2≤2εηNk2+8εCL2ηNk−12+2εL2ηNk−22+εCL2ηNk−Ns~2+εCgk−g~k2.
By means of Lemma12 and since ε>0, there exists a positive constant τ∗=qm/2Γ2−βm2/εβm; when τ<τ∗, we have
(48)ηNk2≤2εCtkβmqmΓ1+βmEβm2μtkβm/qm∑k=1Mgk−g~k2,with μ=2/ε+8CεL2/a0βm−a1βm+2CεL2/a1βm−a2βm+CεL2/aNs−1βm−aNsβm. Thus, the scheme is unconditionally stable.☐


## 5.2. Convergence Analysis

In this section, we investigate the convergence of the fully discrete scheme (38) using error estimation.Theorem 14.
Letukk=−NsM be the exact solution of equation (1) and uNkk=−NsM be the solution of (38). Suppose that D0Ctβru∈L∞0,T;H0α/2Ω∩HsΩ,u∈L∞−s,T;HsΩ; we have
(49)uk−uNkα/2≤CNα/2−s+N−r+τ2−βm,1≤k≤M,where C is independent of N and τ.Proof.
Denoteuk−uNk=eNk=uk−πNα/2,0uk+πNα/2,0uk−uNk=Δe~Nk+e^Nk. The weak formulation of equation (1) is
(50)∑r=0mqr0CDtβruk,υN+Auk,υN=fuk,uk−Ns~,υN+gk,υN.
Subtracting (38) from (50) and owing to the definition of orthogonal projection, the error equation satisfies
(51)∑r=0mqrDτβre^Nk,υN+Ae^N,υN=ΔR1k+R2k+R3k+R4k,where
(52)R1k=INfuk,uk−Ns~−INf2uNk−1−uNk−2,uNk−Ns~,υN,R2k=fuk,uk−Ns~−INfuk,uk−Ns~,υN,R3k=∑r=0mqrDτβrπNα/2,0uk−0CDtβruk,υN,R4k=gk−INgk,υN.
We next estimate the right-hand termsR1k, R2k, R3k, and R4k. For the first term R1k,
(53)R1k=INfuk,uk−Ns~−INf2uk−1−uk−2,uk−Ns~,υN+INf2uk−1−uk−2,uk−Ns~−INf2uNk−1−uNk−2,uNk−Ns~,υN=ΔR11k+R12k.
By applying the Taylor expansion, the following holds:(54)fuk,uk−Ns~=f2uk−1−uk−2,uk−Ns+uk−2uk−1+uk−2f′1ξ,uk−Ns~=f2uk−1−uk−2,uk−Ns~+c~uτ2.Furthermore, by means of the Hölder inequality and Young’s inequality, we have
(55)R11k≤INfuk,uk−Ns~−INf2uk−1−uk−2,uk−Ns~υN≤Cfuk,uk−Ns~−f2uk−1−uk−2,uk−Ns~υN≤ε2c~uτ4+12ευN2.
According to (37), we can deduce that
(56)R12k≤LC2eNk−1−eNk−2+eNk−Ns~υN≤LC2e^Nk−1−e^Nk−2+e^Nk−Ns~+2e~Nk−1−e~Nk−2+e~Nk−Ns~υN≤8ε2CL2e^Nk−12+2ε2CL2e^Nk−22+ε2L2e^Nk−Ns~2+8ε2CL2e~Nk−12+2ε2CL2e~Nk−22+ε2CL2e~Nk−Ns~2+12ευN2.
Moreover, owing to Lemmas6 and 7, the following holds:
(57)e~Nk−12≤CC1Nα−2suk−1s2,e~Nk−22≤CC1Nα−2suk−2s2,e~Nk−Ns~2≤CC1Nα−2suk−Ns~s2.Then, (56) becomes
(58)R12k≤4εCL2e^Nk−12+εCL2e^Nk−22+ε2CL2e^Nk−Ns~2+CC1Nα−2suL∞−s,T;HsΩ2+12ευN2.
Substituting (55) and (58) into (53), we can derive that
(59)R1k≤1ευN2+4εCL2e^Nk−12+εCL2e^Nk−22+ε2CL2e^Nk−Ns~2+CC1Nα−2suL∞−s,T;HsΩ2+ε2c~uτ4.
For the second termR2k, by means of the Hölder inequality, the following holds:
(60)R2k≤ε2CN−2ruL∞−s,T;HsΩ2+12ευN2.
For the third termR3k, the following holds:
(61)R3k=∑r=0mqrDτβrπNα/2,0uk−D0CtβrπNα/2,0uk,υN+∑r=0mqrD0CtβrπNα/2,0uk−0CDtβruk,υN=πNα/2,0∑r=0mqrDτβruk−0CDtβruk,υN−∑r=0mqr0CDtβre~Nk,υN=ΔR31k+R32k.Using (12) and the Hölder inequality, the following holds:
(62)R31k≤ε2∑r=0mqrπNα/2,0Dτβruk−0CDtβruk2+12ευN2≤ε2C∑r=0mqrDτβruk−D0Ctβruk2+12ευN2≤ε2C1,uτ4−2βm+12ευN2.Furthermore, according to Lemma 2, we have
(63)R32k≤ε2CNα−2s∑r=0mqrD0Ctβruks2+12ευN2≤ε2CNα−2s∑r=0mqrD0CtβruL∞−s,T;HsΩ2+12ευN2.
Thus, (61) becomes
(64)R3k≤ε2CNα−2s∑r=0mqr0CDtβruL∞−s,T;HsΩ2+ε2C2,uτ4−2βm+1ευN2.
For the fourth termR4k, the following holds by invoking Remark 9:
(65)R4k≤ε2CNα−2ruL∞−s,T;HrΩ2+12ευN2.
Substituting (59), (60), (64), and (65) into (51), we can infer that
(66)∑r=0mqrDτβre^Nk,υN+Ae^N,υN≤52ευN2+4εCL2e^Nk−12+εCL2e^Nk−22+ε2CL2e^Nk−Ns~2+G~,where
(67)G~=εC~Nα−2s∑r=0mqr0CDtβruL∞0,T;HsΩ2+uL∞−s,T;HsΩ2+εC~N−2ruL∞−s,T;HrΩ2+εC~uτ4−2βm.
TakingυN=e^Nk in (66) and applying Lemma 11, we can conclude that
(68)∑r=0mqr2Dτβre^Nk2+e^Nkα/22≤52εe^Nk2+4εCL2e^Nk−12+εCL2e^Nk−22+ε2CL2e^Nk−Ns~2+G~,namely,
(69)∑r=0mqr2Dτβre^Nk2≤5εe^Nk2+8εCL2e^Nk−12+2εCL2e^Nk−22+εCL2e^Nk−Ns~2+G,with G=2G~. By means of Lemma 12 and since ε>0, there exists a positive constant τ∗=qm/2Γ2−βm10/εβm; when τ<τ∗, we have
(70)e^Nk2≤2GtkβmqmΓ1+βmEβm2μtkβmqm,with μ=10/ε+16εCL2/a0βm−a1βm+4CεL2/a1βm−a2βm+2CεL2/aNs~−1βm−aNs~βm. Thus, the scheme is unconditionally convergent. Finally, by means of the triangle inequality and Lemma 6, we complete the proof of (49).☐


## 6. Numerical Experiments

Some numerical experiments are performed here to clarify the convergence orders of the considered scheme in time and space. We also show the impact of time- and space-fractional orders on the behaviour of the dynamics for the solution of nonlinear delay reaction diffusion equations. To examine the temporal and spatial convergence orders separately, the orders of convergence in time and space shall be determined from theL2-error norms defined as
(71)Order=lneN,M1/eN,M2lnM1/M2,where M1≠M2.Example 15.
Consider the nonlinear delay reaction-diffusion problem(72)∑r=1Q∂βru∂tβrx,t=∂αu∂xαx,t−2ux,t+ux,t−1.51+u2x,t−1.5+gx,t,x∈0,1,t∈0,1.We choose the fractional ordersβr=2Q+r−5/3Q. The source function gx,t is given such that problem (72) has the exact solution tβ5+1x21−x2.In Table1, we list the L2-errors and corresponding convergence orders with α=1.1,1.5,1.9 and N=50 for Q=5. We can see that these results confirm the expected theoretical order convergence in time. The convergence orders in space are depicted for different values of α and β at M=500 in Figure 1. All the convergence results are in agreement with the theoretical results.Example 16.
We consider the following fractional problem, where the dynamics of the solution is very interesting and the exact solution is unknown:(73)∑r=18∂βru∂tβrx,t=∂αu∂xαx,t+ux,t1−ux,t1+ux,t−1.5,x∈a,b,t∈0,3,with the initial value ux,0=e−2x2 and βr=r/10.Table 1
TheL2-errors and their convergence orders versus τ and α at N=50 for Example 15.
τα=1.1α=1.5α=1.9ErrorOrderErrorOrderErrorOrder1.5/505.057×10–5—4.692×10–5—4.163×10–5—1.5/1001.927×10–51.3921.989×10–51.2381.783×10–51.2231.5/1501.094×10–51.3951.192×10–51.2621.086×10–51.2241.5/2007.327×10–61.3968.283×10–61.2677.592×10–61.2441.5/2505.365×10–61.3966.226×10–61.2795.732×10–61.2591.5/3004.159×10–61.3974.933×10–61.2774.560×10–61.2552–β5—1.333—1.333—1.333Figure 1
Convergence order in space direction for different values ofα at τ=1.5/300.Figures2 and 3 show the profiles of the numerical solution with the fractional-order parameters α=1.2,1.5,1.8 with N=100 and τ=1.5/500. We can observe that the fractional-order parameter α affects the shape of the solutions. We can say that the fractional-order parameters can be used in physics to modify the shape of waves without changing the nonlinearity and dispersion effects of the fractional nonlinear problems.Figure 2
Evaluations ofu for different fractional-order parameters α.Figure 3
Evaluations ofu at x=0 and t=T for different fractional-order parameters α.


## 7. Conclusion and Remarks

We have constructed and analyzed a novel explicit finite difference/Galerkin–Legendre spectral scheme for the nonlinear multiterm Riesz space and Caputo time fractional reaction-diffusion equation with delay. The problem was first approximated by theL1 difference method on the temporal direction, and then, the Galerkin–Legendre spectral method was applied on the spatial discretization. Using an appropriate form of discrete fractional Grönwall inequality, the stability and the convergence of the fully discrete scheme were investigated. We have proven that the proposed method is stable and has a convergent order 2−βm in time and an exponential rate of convergence in space. High-order difference schemes can be handled to raise the temporal convergence order. This can be done by using the Alikhanov scheme [46], and it can be designed theoretically and numerically easily as shown in our context. Two numerical examples are given to show that the numerical results are consistent with theoretical ones in the case of the smoothness of the solution with respect to time and space.



---

*Source: 5922853-2021-12-06.xml*