#### 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 as . The problem is first approximated by the 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 in time and an exponential rate of convergence in space. We finally provide some numerical experiments to show the efficacy of the theoretical results.

#### 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 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: endowed with initial-boundary conditions of the form

Here, and are space and time domains, respectively. We denote as the Caputo fractional derivative with the fractional orders . The parameters are positive constants. The parameters are absolutely positive. Also, is the space fractional order. The left and right Riemann-Liouville fractional derivatives of order on the infinite domain [34] are defined as where is the gamma function. Thus, the space fractional derivative in the Riesz form on the space interval can be defined as [35]

The Caputo derivative is defined as

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 uniform -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 Section 3. 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) is the inner product on the space with the -norm (ii)The maximum norm is defined as (iii) is the space of smooth functions with compact support in (iv) and are the usual Sobolev spaces with the norm and seminorm (v) is the space of polynomials defined on the domain with degree at most (vi)The approximation space is defined as (vii) is the Legendre-Gauss-Lobatto interpolation operator as

*Definition 1. *Fractional derivative spaces and their related norms and seminorms are defined as follows [36]:
(i)Left fractional space: let . Define the seminorm and the norm , and let (or ) denote the closure of (or ) with respect to (ii)Right fractional space: let Define the seminorm and the norm , and let (or ) denote the closure of (or ) with respect to (iii)Symmetric fractional space: let Define the seminorm and the norm , and let (or ) denote the closure of (or ) with respect to (iv)Fractional Sobolev space: let . Define the fractional Sobolev space as , endowed with the seminorm and the norm , where is the Fourier transformation of and is the extension of zero of outside . Denote by (or ) the closure of (or ) with respect to .

Lemma 2 (see [36]). *The spaces , and are equivalent, with equivalent seminorms and norms if .*

Lemma 3 (adjoint property). *Let , then for any and , we get
*

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 . Denote as the -th order Jacobi polynomial of index defined on . As all classic orthogonal polynomials, satisfies the following three-term-recurrence relation:

The recursion coefficients are given by

Let . Then, one has where is the Kronecker delta function and

In particular, the Legendre polynomial is defined as .

#### 3. The Numerical Scheme

Here, we provide a fully discrete scheme for the problems (1) and (2) based on the -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 uniformly with a time step size defined by such that is a positive integer. The uniform partitions given by where The interpolation scheme for the time-fractional derivative of order , , in the Caputo sense at the time is defined as where and , for each . If , then there exists a constant such that the truncation error satisfies , for each (see [41]).

*Definition 4. *Let be a sequence of real functions defined on . We define the multiterm discrete time-fractional difference operator by

In this expression, , and the constants are defined by , and , for each .

In order to provide a semidiscretized form of (1) at each time , 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:

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]. where for each , the function is given by and .

We introduce the parameter . Then, the scheme (14) can be rewritten in the following equivalent form:

The fully discrete -Galerkin spectral scheme consists of the set of approximations , satisfying the system: where is an appropriate projection operator. We expand the approximate solution as

Substituting this expression into (19) and letting , for each , we obtain the following matrix representation of the uniform -Galerkin spectral scheme:

The notations in this expression are given by the system of identities:

Lemma 5 (see [42, 43]). *The components of the stiffness matrix are , for each . Here,
and is the set of Jacobi–Gauss points and weights with respect to the weight function . The mass matrix is symmetric, with nonzero components:
*

#### 4. Technical Lemmas

Several lemmas that will be invoked through our analysis appeared in that section. In the sequel, and will denote generic positive constants independent of , , and and may be different under different circumstances. We also fix the following notation , such that is the set of all integers.

Throughout the coming context, we will use the notation

The orthogonal projection operator will be such that

For convenience of theoretical analysis, we give the following seminorm and norm: which are equivalent to the seminorms and norms of , and . We recall the following three lemmas from [43].

Lemma 6. *Let and be arbitrary real numbers satisfying . Then, there exists a positive constant C independent of such that, for any function , the following estimate holds:
*

Lemma 7. *Suppose that , . Then, there exist positive constants and independent of , such that
*

The following lemma and remark summarize the properties of the interpolation operator .

Lemma 8 (see [44]). *Let . If , then there exists a constant independent of , such that , for any .*

*Remark 9. *A smooth solution of a fractional differential equation does not mean a smooth source term and vice versa. Therefore, the regularity order of the solution is not the same as the regularity order of the source term , i.e.,

Lemma 10 (see [45]). *For any function which is absolutely continuous on , the following inequality is satisfied:
*

Lemma 11 (see [46]). *The discrete counterpart to the inequality (32) is given as
such that is the discrete time-fractional difference operator of the 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 -schemes for the multiterm problems with nonlinear delay.

Lemma 12 (discrete fractional Grönwall inequality [29]). *Let be nonnegative sequences. Let , and be positive constants independent of . The fractional orders are defined as . If is known and if ,
**where and are positive constants. Then, there exists a positive constant such that
where is the Mittag-Leffler function and
*

#### 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 satisfies the following Lipschitz condition where is a positive constant.

##### 5.1. Stability Analysis

The weak formulation of the scheme is as follows: find , such that

with

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 that is the solution of with initial conditions

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 , the following holds:
*

*Proof. *Denote . Subtracting (40) from (38), the following holds:
According to (37) and using the Hölder inequality and Young’s inequality, we derive that
Then, (43) becomes
Taking and using Lemma 11 and (27), we can deduce that
namely,
By means of Lemma 12 and since , there exists a positive constant ; when , we have
with . 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. *Let be the exact solution of equation (1) and be the solution of (38). Suppose that ; we have
where is independent of and *

*Proof. *Denote . The weak formulation of equation (1) is
Subtracting (38) from (50) and owing to the definition of orthogonal projection, the error equation satisfies
where
We next estimate the right-hand terms , , , and . For the first term ,
By applying the Taylor expansion, the following holds:
Furthermore, by means of the Hölder inequality and Young’s inequality, we have
According to (37), we can deduce that
Moreover, owing to Lemmas 6 and 7, the following holds:
Then, (56) becomes
Substituting (55) and (58) into (53), we can derive that
For the second term , by means of the Hölder inequality, the following holds:
For the third term , the following holds:
Using (12) and the Hölder inequality, the following holds: