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

Department of Computational Mathematics and Computer Science, Institute of Natural Sciences and Mathematics Ural Federal University, 19 Mira St., Yekaterinburg 620002, Russia Department of Mathematics, Faculty of Science, Al-Azhar University, Assiut 71524, Egypt Department of Mathematics, Nazarbayev University, Nur-Sultan, Kazakhstan Department of Applied Mathematics, Physics Division, National Research Centre, Dokki, Cairo 12622, Egypt Department of Mathematics, Faculty of Science, Benha University, Benha 13511, Egypt Institute of Mathematics and Mechanics, Ural Branch of the Russian Academy of Sciences, 16 Kovalevskoy St., Yekaterinburg 620000, Russia


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][2][3][4][5][6][7][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][14][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 timefractional 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 advectiondiffusion 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 mul-titerm time-space fractional reaction-diffusion equations with delay: endowed with initial-boundary conditions of the form 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 q r 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 where ΓðxÞ 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 ∂ β /∂t β 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 reactiondiffusion problem with fixed delay. A hybrid numerical scheme combines the Galerkin-Legendre spectral schemes, and a uniform L1-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 2 Journal of Function Spaces 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.

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: Spectral methods are characterized by the expansion of the solution in terms of global and, usually, orthogonal polynomials [37][38][39][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 J μ,ν i ðxÞ, μ, ν > −1 as the i-th order Jacobi polynomial of index defined on ½−1, 1. As all classic orthogonal polynomials, fJ μ,ν i ðxÞg N i=0 satisfies the following three-termrecurrence relation:

Journal of Function Spaces
The recursion coefficients are given by Let ω μ,ν ðxÞ = ð1 − xÞ μ ð1 + xÞ ν . Then, one has where δ i,j is the Kronecker delta function and In particular, the Legendre polynomial is defined as

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/N~s such that N~s is a positive integer. The uniform partitions given by t n = nτ, ∀−N~s ≤ n ≤ M, where M = dT/τe: 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 t n is defined as where ω β r ðtÞ = t β r −1 /Γðβ r Þ, t > 0, and a β r j = ðj + 1Þ 1−β r − j 1−β r , for each j ≥ 0. If u ∈ C 2 ð½0, T ; L 2 ðΩÞÞ, then there exists a constant C > 0 such that the truncation error r n τ satisfies kr n τ k ≤ Cτ 2−β m , for each n = 0, 1, ⋯, M (see [41]).
Definition 4. Let fu n g M n=0 be a sequence of real functions defined on Ω. We define the multiterm discrete time- In this expression, δ t u i = u i − u i−1 , and the constants are defined by b In order to provide a semidiscretized form of (1) at each time t n , 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 spacefractional differential equations [42,43].
where for eachx ∈ ½−1, 1, the function φ n is given by We introduce the parameter σ r = q r /Γð2 − β r Þτ β r . Then, the scheme (14) can be rewritten in the following equivalent form: Journal of Function Spaces The fully discrete L1-Galerkin spectral scheme consists of the set of approximations u n N ∈ W 0 N , satisfying the system: where π 1,0 N is an appropriate projection operator. We expand the approximate solution as Substituting this expression into (19) and letting υ = φ k , for each 0 ≤ k ≤ N − 2, we obtain the following matrix repre-sentation of the uniform L 1 -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 and fx −α/2,−α/2 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:

Technical Lemmas
Several lemmas that will be invoked through our analysis appeared in that section. In the sequel, C and C u 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.

Journal of Function Spaces
Throughout the coming context, we will use the notation The orthogonal projection operator π α/2,0 N : H α/2 0 ðΩÞ ⟶ W 0 N 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 J α/2 L ðΩÞ, J α/2 R ðΩÞ, J α/2 S ðΩÞ, 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 ∈ H α/2 0 ðΩÞ ∩ H s ðΩÞ, the following estimate holds: Lemma 7. Suppose that Ω = ða, bÞ, u ∈ H α/2 0 ðΩÞ. Then, there exist positive constants C 1 < 1 and C 2 independent of u, such that The following lemma and remark summarize the properties of the interpolation operator I N .
Lemma 8 (see [44]). Let s ≥ 1. If u ∈ H s ðΩÞ, then there exists a constant C > 0 independent of N, such that ku − I N uk l ≤ CN l−s kuk s , 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 order s of the solution u is not the same as the regularity order r of the source term g, i.e., Lemma 10 (see [45]). For any function uðtÞ which is absolutely continuous on ½0, T, the following inequality is satisfied: Lemma 11 (see [46]). The discrete counterpart to the inequality (32) is given as such that D β τ u k 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][48][49][50] tried to fill that gap. In what follows, we present recent discrete fractionaltype inequalities. These inequalities play an important role in analyzing stability and convergence of the L1-schemes for the multiterm problems with nonlinear delay.

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 where L is a positive constant.

Stability Analysis.
The weak formulation of the scheme is as follows: 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 fũ k N g

M k=1
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 τ > 0, the following holds: Proof. Denote η k N = u k N −ũ k N . 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 υ N = η k N and using Lemma 11 and (27), we can deduce that ; when τ < τ * , we have with μ = 2/ε + ð8CεL 2 /ða Thus, the scheme is unconditionally stable.☐

Convergence Analysis.
In this section, we investigate the convergence of the fully discrete scheme (38) using error estimation.  (38). Suppose that where C is independent of N and τ: 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 R k 1 , R k 2 , R k 3 , and R k 4 . For the first term R k 1 , 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: Journal of Function Spaces Then, (56) becomes Substituting (55) and (58) into (53), we can derive that For the second term R k 2 , by means of the Hölder inequality, the following holds: For the third term R k 3 , the following holds: Using (12) and the Hölder inequality, the following holds: Furthermore, according to Lemma 2, we have Thus, (61) becomes For the fourth term R k 4 , the following holds by invoking Remark 9: Substituting (59), (60), (64), and (65) into (51), we can infer that wherẽ Taking υ N =ê k N in (66) and applying Lemma 11, we can conclude that with μ = 10/ε + ð16εCL 2 /ða Ns ÞÞ. Thus, the scheme is unconditionally convergent. Finally, by means of the triangle inequality and Lemma 6, we complete the proof of (49).☐

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 spacefractional orders on the behaviour of the dynamics for the solution of nonlinear delay reaction diffusion equations. To examine the temporal and spatial convergence orders sepa-rately, the orders of convergence in time and space shall be determined from the L 2 -error norms defined as where Example 15. Consider the nonlinear delay reaction-diffusion problem We choose the fractional orders β r = ð2Q + r − 5Þ/3Q. The source function gðx, tÞ is given such that problem (72) has the exact solution t β 5 +1 x 2 ð1 − xÞ 2 :

Journal of Function Spaces
In Table 1, we list the L 2 -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: with the initial value uðx, 0Þ = e −2x 2 and β r = r/10: Figures 2 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.

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 the L1 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. Highorder 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.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no competing interests.