Journal of Computational and Applied Mathematics 420 (2023) 114832 Contents lists available at ScienceDirect Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam An L1 type difference/Galerkin spectral scheme for variable-order time-fractional nonlinear diffusion–reaction equationswith fixed delay M.A. Zaky a,∗, K. Van Bockstal b, T.R. Taha c, D. Suragan d, A.S. Hendy e,f,∗∗ a Department of Applied Mathematics, National Research Centre, Dokki, 33 El-Bohouth St., Giza 12622, Egypt b Research Group NaM2 , Department of Electronics and Information systems, Ghent University, Krijgslaan 281, 9000 Ghent, Belgium c Department of Computer Science, University of Georgia, Athens, GA 30602-7404, USA d Department of Mathematics, Nazarbayev University, Nur-Sultan, Kazakhstan e Department of Computational Mathematics and Computer Science, Institute of Natural Sciences and Mathematics, Ural Federal University, 19 Mira St., Yekaterinburg 620002, Russia f Department of Mathematics, Faculty of Science, Benha University, Benha 13511, Egypt a r t i c l e i n f o Article history: Received 30 December 2021 Received in revised form 24 May 2022 Keywords: Variable order diffusion Time delay L1 difference scheme Galerkin spectral method Convergence and stability estimates a b s t r a c t A linearized spectral Galerkin/finite difference approach is developed for variable fractional-order nonlinear diffusion–reaction equations with a fixed time delay. The temporal discretization for the variable-order fractional derivative is performed by the L1-approximation. An appropriate basis function in terms of Legendre polynomials is used to construct the Galerkin spectral method for the spatial discretization of the second-order spatial operator. The main advantage of the proposed approach is that the implementation of the iterative process is avoided for the nonlinear term in the variable fractional-order problem. Convergence and stability estimates for the con- structed scheme are proved theoretically by discrete energy estimates. Some numerical experiments are finally provided to demonstrate the efficiency and accuracy of the theoretical findings. © 2022 Elsevier B.V. All rights reserved. 1. Introduction Fractional theory of calculus of variable-order has recently attracted considerable attention as a powerful mathematical tool for modeling a wide range of discontinuities and nonlinear phenomena [1]. Because many physical systems involve dynamics with memory effects whose behavior changes over time, even while transiting from one fractional- order to another, interest in fractional operators moved increasingly to their variable-order counterparts. The powerful practical implications of these variable-order objects, of course, come at the price of a more complicated mathematical characterization. In contrast to integer-order and fixed fractional-order operators, the variable-order of the fractional operators can be considered as a function of internal or external system variables such as, for example, time, space, state of stress, system energy, temperature, or even a combination of the various variables. Such an extension from ∗ Corresponding author. ∗∗ Corresponding author at: Department of Computational Mathematics and Computer Science, Institute of Natural Sciences and Mathematics, Ural Federal University, 19 Mira St., Yekaterinburg 620002, Russia. E-mail addresses: ma.zaky@yahoo.com (M.A. Zaky), karel.vanbockstal@ugent.be (K. Van Bocksta), trtaha@uga.edu (T.R. Taha), durvudkhan.suragan@nu.edu.kz (D. Suragan), ahmed.hendy@fsc.bu.edu.eg (A.S. Hendy). https://doi.org/10.1016/j.cam.2022.114832 0377-0427/© 2022 Elsevier B.V. All rights reserved. https://doi.org/10.1016/j.cam.2022.114832 http://www.elsevier.com/locate/cam http://www.elsevier.com/locate/cam http://crossmark.crossref.org/dialog/?doi=10.1016/j.cam.2022.114832&domain=pdf mailto:ma.zaky@yahoo.com mailto:karel.vanbockstal@ugent.be mailto:trtaha@uga.edu mailto:durvudkhan.suragan@nu.edu.kz mailto:ahmed.hendy@fsc.bu.edu.eg https://doi.org/10.1016/j.cam.2022.114832 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 o f v p c t t u a S h M d f t o p d a i t o d c f m t K r d t a f f f f a o e w w s w fixed-order to variable-order operators also allows the formulation of mathematical models in which the order of the underlying governing equations can be modified using either the system’s instantaneous state or its history. As a result, the corresponding model can evolve seamlessly to capture widely dissimilar dynamics without changing the structure of the governing equations which describe the system’s response [2]. This emphasizes the evolutionary nature of the variable- order fractional calculus formalism, which indeed can play a critical role in the simulation of nonlinear dynamical models. Recognizing this untapped potential and unique capability of variable-order operators, the scientific community has been intensively investigating applications of variable-order fractional calculus to the modeling of physical and engineering systems [3–9]. It is worth noting that viscoelasticity is certainly the playground of the most widespread applications of variable- rder fractional calculus since its appearance [10–12]. Although some examples of fractional operators with constant ractional-order successfully mimic experimental test findings, these operators may not always be suited for characterizing iscoelastic behavior. Indeed, nonlinear and complex phenomena occur in several materials, such as biological tissues, olymers, and rubbers. These phenomena imply significant variations in the material’s mechanical characteristics that annot be described by fractional viscoelastic models of constant-orders [2]. Similar problems arise to characterize he mechanical behavior of viscoelastic materials under variable environmental conditions, for instance, changes in emperature or viscoelastic aging effects. As a result, variable-order fractional operators have recently been proposed and tilized in a wide range of applications. For more comprehensive literature review on variable-order fractional operators, ssociated differential equations and latest applications in natural sciences, we refer the reader to Patnaik et al. [13], amko [14], Sun et al. [1] and Ortigueira et al. [15]. Time-fractional diffusion equations serve as effective tools for modeling and simulation of various processes with ereditary or memory features, such as anomalous diffusion transport, and hence attract considerable interest [16–19]. oreover, the theory, application, and implementation of several numerical approaches for the solution of time-fractional ifferential equations have been proposed in [20–25]. However, the singularity of the solutions of constant-order time- ractional diffusion equations at the beginning of time seems to be physically irrelevant to the subdiffusive transport of he model. The fundamental reason why this phenomenon occurs lies in the incompatibility between the global nature f the power law decaying tails and the locality of the classical initial condition at the beginning of time [26,27]. It was ointed out in [28] that the regularity of the solution to the variable-order Riemann–Liouville fractional diffusion equation epends on the behavior of the variable order and its derivatives at time t = 0, in addition to the standard smoothness ssumptions. It was demonstrated that the solution to the problem exhibits full regularity like its integer-order analogue f the variable order has an integer limit at t = 0, or it has a singularity at t = 0 like in the case of the constant-order ime-fractional diffusion equations if the variable order has a non-integer value at time t = 0. Accordingly, the variable- rder fractional operators can eliminate the nonphysical singularity of the solutions to constant-order time-fractional iffusion equations and produce solutions with full regularity. To the best of our knowledge, up to now, there are no losed form solutions to variable-order fractional diffusion equations. The main reason is that, due to the impact of variable ractional order, analytic techniques such as the Laplace transform cannot be utilized to solve variable-order fractional odels. Recently, there are several investigations on numerical schemes to approximate the solutions of variable-order ime-fractional diffusion equations in the literature. We mention here the works by Tavares et al. [29], Chen et al. [30], arniadakis et al. [31–33], Zheng et al. [34–36], Pang and Sun [37], Wei et al. [38,39] and Liu et al. [40]. Nevertheless, igorous numerical and mathematical analysis of variable-order fractional differential equations remains wide-open. Most studies in the literature have been absorbed in the linear variable-order fractional diffusion equations without elay [27,41,42]. Although the generalized nonlinear ones are much more useful in real applications, the existing works on he numerical solutions of such problems are quite sparse. Fractional differential equations with delay have recently played significant role in modeling of many real areas of sciences. A consideration of a method of backward differentiation ormula type for solving them is done in [43]. Most recent studies investigated the finite difference/spectral solutions of the ractional diffusion–reaction equations only in implicit schemes [23,44,45]. In this work, we develop a linearized explicit inite difference/spectral Galerkin approach for nonlinear variable fractional-order diffusion–reaction equations with a ixed time delay and a drift term, and discuss its stability and the convergence rate. The main advantage of the proposed pproach is that the implementation of the iterative process is avoided for the nonlinear term in the variable fractional- rder problem. More specifically, we consider the following nonlinear variable order time-fractional reaction–diffusion quations with delay ∂β(t)u ∂tβ(t) + λ ∂u ∂t = κ ∂2u ∂x2 + f (u(x, t), u(x, t − s)) + g(x, t), 0 < β(t) ≤ β̄ < 1, x ∈ Ω, t ∈ I, (1.1) ith the initial–boundary conditions of the form{ u(x, t) = φ(x, t), x ∈ Ω, t ∈ [−s, 0], u(a, t) = u(b, t) = 0, t ∈ I, (1.2) here I = (0, T ] ⊂ R and Ω = (a, b) ⊂ R are the time and the space domains, respectively. The parameters κ, λ, and are positive constants. Throughout the paper, we assume that g is a continuous function and f is Lipschitz continuous ith the Lipschitz constant L, that is, |f (u , v ) − f (u , v )| ≤ L |u − u | + |v − v | , ∀u , v ∈ R. 1 1 2 2 ( 1 2 1 2 ) i i 2 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 d T i a d t a t m s o c H p d i m t e 2 − w e D In addition, we suppose that φ is a smooth function which vanishes at the endpoints of Ω . The variable-order fractional integral operator 0I β(τ ) τ , and the variable-order fractional Caputo operator ∂β(τ ) ∂τβ(τ ) are efined as [46,47] 0Iβ(τ )τ u(τ ) := 1 Γ (β(τ )) ∫ τ 0 u(r) (τ − r)1−β(τ ) dr, ∂β(τ ) ∂τ β(τ ) u(τ ) := 0I1−β(τ ) τ u′(τ ) = 1 Γ (1 − β(τ )) ∫ τ 0 u′(r) (τ − r)β(τ ) dr. (1.3) he well-posedness of the linear hidden-memory variable–order Caputo time-fractional diffusion equations was discussed n [48,49], and the space-dependent variable-order fractional Caputo case was discussed in [42]. The well-posedness nd the regularity of the solutions to Riemann–Liouville variable-order fractional nonlinear differential equations were iscussed in [50,51]. As a special class of the problem under consideration, an initial boundary value problem for he fractional delayed semilinear diffusion equation with the fractional Laplacian was discussed in [52]. The existence nd uniqueness of mild solutions for the abstract time–space evolution equation with delay are investigated by using he semigroup theory of operators and the monotone iterative technique under some quasi-monotone conditions. The aximum principle for multi-term space–time variable-order Riesz–Caputo fractional diffusion equations with nonlinear ource function was discussed in [53]. This principle was employed to show the uniqueness of solutions of the variable- rder Riesz–Caputo fractional diffusion equations and the continuous dependence of solutions on initial–boundary value onditions. The existence of a solution to nonlinear equations with delay of type (1.1) is a fundamental research question. owever, we leave this as an open problem for future work. To the best of our knowledge, we can construct a maximum rinciple for the problem under consideration by invoking the methodology in [53] and by noticing the existence of a elay in (1.1)–(1.2). This can also be done by invoking the techniques of Barrier analysis in [52]. In the meantime, we nvestigate the problem’s numerical solution and its convergence and stability analysis. The rest of the paper is arranged as follows. In Section 2, we construct the fully discrete finite difference/Galerkin ethod for the above model equation. In Section 3, we prove that the fully discrete scheme is unconditionally stable, and he numerical approximation is convergent. In Section 4, we perform two numerical examples to show the accuracy and fficiency of the proposed scheme. Finally, we give the concluding remarks in Section 5. . The linearized numerical scheme Let τ = s/Ns be the temporal step size, where Ns is a positive integer. Denote M = ⌈ T τ ⌉ . Define tk = kτ , for each Ns ≤ n ≤ M and uk = u(·, tk). The L1-approximation is defined as [54] ∂β(t)f (x, t) ∂tβ(t) |t=tk = ∫ tk 0 f ′(x, z)ω1−β(tk)(tk − z)dz = 1 Γ (1 − β(tk)) k∑ q=1 f (x, tq) − f (x, tq−1) τ ∫ tq tq−1 (tk − z)−β(tk)dz + rkτ = 1 τ β(tk)Γ (2 − β(tk)) k∑ q=1 akk−q ( f (x, tq) − f (x, tq−1) ) + rkτ , (2.4) here ωβ(t)(t) = tβ(t)−1 Γ (β(t)) , t > 0, akq = (q+ 1)1−β(tk) − q1−β(tk), for each q ≥ 0. Let f ∈ C2([0, T ]; L2(Ω)). Then the truncation rror rkτ fulfills ∥rkτ ∥ ≤ Cτ 2−β̄ , for each k = 0, 1, . . . ,M . efinition 1. Define [54] Dβk τ uk = τ 1−βk Γ (2 − βk) k∑ q=1 akk−qδtu q = τ−βk Γ (2 − βk) k∑ q=0 bkk−qu q, ∀k = 1, . . . ,M. (2.5) In this expression, we define δtuq = uq−uq−1 τ and bk0 = ak0, bkk = −akk−1, bkk−q = akk−q − akk−q−1, for q = 1, . . . , k − 1. The coefficients akq satisfy 1 = ak0 > ak1 > ak2 > · · · > akk > 0, (2.6) and (1 − βk)(q + 1)−βk ≤ ak ≤ (1 − βk)q−βk . q 3 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 T t w t w L s w According to (2.6) and [55], the following property is satisfied:( Dβk τ uk, uk ) ≥ 1 2 Dβk τ (uk)2. (2.7) he nonlinear source term is discretized and linearized by Taylor’s expansion. As a result, we get the following ime-discrete system λ δtum + Dβm τ um = κ ∂2um ∂x2 + f (−um−2 + 2um−1, um−Ns ) + gm, 1 ≤ m ≤ M, ∀x ∈ Ω; (2.8) um(x) = φ(x), −Ns ≤ m ≤ 0, x ∈ Ω. We here introduce some fundamental properties of Jacobi polynomials. Let N be a positive integer. Denote Pγ ,ς q (y), γ ith ς > −1 as the qth order Jacobi polynomial of indices γ , ς defined on (−1, 1). They satisfy the following hree-term-recurrence relation⎧⎪⎪⎨⎪⎪⎩ Pγ ,ς 0 (y) = 1, Pγ ,ς 1 (y) = 1 2 (2 + γ + ς )y + 1 2 (γ − ς ), Pγ ,ς q+1(y) = ( Aγ ,ς q y − Bγ ,ς q ) Pγ ,ς q (y) − Cγ ,ς q Pγ ,ς q−1(y), if 1 ≤ q ≤ N, here⎧⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎩ Aγ ,ς q = (2q + γ + ς + 1)(2q + γ + ς + 2) 2(q + 1)(q + γ + ς + 1) , Bγ ,ς q = (2q + γ + ς + 1)(ς2 − γ 2) 2(q + 1)(q + γ + ς + 1)(2q + γ + ς ) , Cγ ,ς q = (2q + γ + ς + 2)(q + γ )(q + ς ) (q + 1)(q + γ + ς + 1)(2q + γ + ς ) . et ωγ ,ς (y) = (1 − y)γ (1 + y)ς . Then, one has∫ 1 −1 Pγ ,ς q (y)Pγ ,ς j (y)ωγ ,ς (t)dy = hγ ,ς q δq,j, ∀q = 0, 1, 2, . . . , where δq,j is the Kronecker delta and hγ ,ς q = 2(γ+ς+1)Γ (q + γ + 1)Γ (q + ς + 1) (2q + γ + ς + 1)q!Γ (q + γ + ς + 1) , ∀q = 0, 1, 2, . . . . For the special case γ = ς = 0, we get the Legendre polynomials Lr (t) = P0,0 r (t). As a result, we introduce the following pace to provide adequate base functions satisfying the boundary conditions [56,57] W0 N = span {ϕs(x) : s = 0, 1, . . . ,N − 2} , here, for each x ∈ [a, b], the function ϕs is given by ϕs(x) = Ls(x̂) − Ls+2(x̂) = 2s + 3 2(s + 1) (1 − x̂2)P1,1 s (x̂), where x̂ := 2x−b−a b−a ∈ [−1, 1]. Next, we introduce the coefficients dk = λ τ + τ−βk Γ ( 2 − βk ) , b̂ki = λ τ δi,k−1 − τ−βk Γ ( 2 − βk )bkk−i. Hence, the scheme (2.8) can be rewritten in the following equivalent form: dmum − κ ∂2um ∂x2 = m−1∑ i=0 b̂mi u i + f (−um−2 − 2um−1, um−Ns ) + gm, ∀m = 1, . . . ,M. The fully scheme includes the set of approximations uk N ∈ W0 N , satisfying the system⎧⎪⎪⎪⎨⎪⎪⎪⎩ dk ( uk N , υ ) − κ ( ∂2 ∂x2 uk N , υ ) = ∑k−1 i=0 b̂ki ( ui N , υ ) + ( IN f (2uk−1 N − uk−2 N , uk−Ns N ), υ ) + ( INgk, υ ) , ∀υ ∈ W0 N , ∀k = 1, . . . ,M, k 1,0 (2.9) uN = πN φ(tk, x), −Ns ≤ k ≤ 0, 4 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 n L d where π 1,0 N is an appropriate projection operator and IN is the Lobatto–Gauss type Legendre interpolation operator. The umerical solution can be expanded as uk N = N−2∑ m=0 ûk mϕm. emma 2.1 (Lemma 4.2, [58]). The stiffness matrix S and the mass matrix M̄ which invoked later in (2.10) are respectively a iagonal matrix and symmetric matrix with the nonzero elements. There elements can be specified as sii = 4i + 6, i = 0, 1, . . . , mij = mji = ⎧⎨⎩ b−a 2j+1 + b−a 2j+5 , i = j, − b−a 2j+5 , i = j + 2. Substituting this formula into (2.9) and taking υ = ϕk, we get the following matrix form representation of the proposed scheme( dk M̄ + κS ) Uk = K k−1 + Rk−1 + Gk, (2.10) where⎧⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ sij = ∫ Ω ϕ′ i (x)ϕ ′ j (x)dx, S = ( sij )N−2 i,j=0 , mij = ∫ Ω ϕi(x)ϕj(x)dx, M̄ = ( mij )N−2 i,j=0 , gk−1 i = ∫ Ω ϕi(x) ( INgk) (x)dx, Gk = (gk−1 0 , gk−1 1 , . . . , gk−1 N−2) ⊤, hk−1 i = ∫ Ω ϕi(x) ( IN f (2uk−1 N − uk−2 N , uk−Ns N ) ) (x)dx, Rk−1 = (hk−1 0 , hk−1 1 , . . . , hk−1 N−2) ⊤, Uk = (ûk 0, û k 1, . . . , û k N−2) ⊤, K k−1 = − ∑k−1 j=0 b̂kk−jM̄U j. 3. Theoretical analysis of the fully discrete scheme The goal of this section is to investigate the effectiveness of the proposed approach for approximating the solution to problem (1.1)–(1.2). In the first subsection, we discuss the stability analysis and provide the stability theorem. In the second subsection, we discuss the convergence analysis. We begin with introducing some technical lemmas that will be of great importance in the following context. 3.1. Technical lemmas This section introduces several lemmas that will be used in our analysis [see, Section 3 [58]]. In the sequel, C and Cv will be used to denote generic positive constants that are independent of n, N , and τ and may vary depending on the conditions. We shall also utilize the following notation B(v, f ) = κ (∂xv, ∂xf ) . (3.11) The orthogonal projection operator π 1,0 N : H1 0 (Ω) → W0 N satisfies B(v − π 1,0 N v, f ) = 0, v ∈ H1 0 (Ω), ∀ f ∈ W0 N . We define the following semi-norms and norms for theoretical purposes [59]: |v|B := B(v, v)1/2, (3.12) ∥v∥B := (∥v∥ 2 + |v| 2 B) 1/2 . We recall the following three lemmas from [56]. 5 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 e L R s s L t w L L P H L Lemma 3.1. Let L be an arbitrary real number satisfying L > 1. Then, for any function w ∈ H1 0 (Ω) ∩ HL(Ω), the following stimate holds true |w − π 1,0 N w|B ≤ CN1−L ∥w∥HL . The following lemma and remark summarize the properties of the interpolation operator IN , see [44,56,58]. emma 3.2. Let L ≥ 1. If w ∈ HL(Ω) then ∥w − INw∥Hq ≤ CNq−L ∥w∥HL , ∀0 ≤ q ≤ 1. emark 3.1. A smooth solution to a partial differential equation with fractional order does not always imply a smooth ource term, and vice versa. As a result, the regularity order L of the solution u differs from the regularity order r of the ource term g , i.e. ∥g − INg∥ ≤ CN−r ∥g∥Hr , ∀g ∈ Hr (Ω), ∥INg∥ ≤ C∥g∥. The following discrete Grönwall inequality will be needed in our analysis later. emma 3.3. Assume that {uk |k = 0, 1, . . . ,M} be non-negative sequence, and it satisfies uk ≤ A + Bτ k∑ s=1 us, k = 0, 1, . . . ,M, hen, when τ ≤ 1 2B , it holds that uk ≤ A exp(2Bkτ ), k = 0, 1, . . . ,M, here A and B are non-negative constants. The following lemma is a direct consequence of Young inequality and inner product properties. emma 3.4. Suppose that {uk N} M k=1 ∈ W0 N , then it holds that( δtuk N , uk N ) ≥ 1 2τ ( ∥uk N∥ 2 − ∥uk−1 N ∥ 2) . emma 3.5. Suppose that {uk N} M k=1 ∈ W0 N , then it holds that( Dβk τ uk N , uk N ) + (δtuk N , uk N ) ≥ 1 2τ ( ∥uk N∥ 2 − ∥uk−1 N ∥ 2) + τ−βk 2Γ (2 − βk) [ ∥uk N∥ 2 − k−1∑ i=1 (akk−i−1 − akk−i)∥u i N∥ 2 − akk−1∥u 0 N∥ 2 ] . roof. A direct consequence of Lemma 3.4 and Eq. (2.7) gives that( Dβk τ uk N , uk N ) + (δtuk N , uk N ) ≥ 1 2τ ( ∥uk N∥ 2 − ∥uk−1 N ∥ 2) + 1 2 Dβk τ ∥uk N∥ 2. ence, the result follows from Eq. (2.5). □ emma 3.6. Suppose that {uk N} M k=1 ∈ W0 N , then it holds that k∑ p=1 ( ∥up N∥ 2 − ∥up−1 N ∥ 2 ) + k∑ p=1 τ−βp Γ (2 − βp) [ ∥up N∥ 2 − p−1∑ i=1 (app−i−1 − app−i)∥u i N∥ 2 − app−1∥u 0 N∥ 2 ] ≥ ( ∥uk N∥ 2 − ∥u0 N∥ 2) + 1 Tβk Γ (1 − βk) k∑ i=1 ∥ui N∥ 2 + ∥u0 N∥ 2 k∑ p=1 τ−βp Γ (2 − βp) app−1. 6 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 3 I t P Proof. Invoking (2.20) in [60], we similarly have that I := k∑ p=1 ( ∥up N∥ 2 − ∥up−1 N ∥ 2 ) + k∑ p=1 τ−βp Γ (2 − βp) [ ∥up N∥ 2 − p−1∑ i=1 (app−i−1 − app−i)∥u i N∥ 2 − app−1∥u 0 N∥ 2 ] = k∑ p=1 ( ∥up N∥ 2 − ∥up−1 N ∥ 2 ) + k∑ p=1 p∑ i=1 τ−βp Γ (2 − βp) app−i ( ∥ui N∥ 2 − ∥ui−1 N ∥ 2) = ( ∥uk N∥ 2 − ∥u0 N∥ 2) + k∑ i=1 τ−βk Γ (2 − βk) akk−i∥u i N∥ 2 + ∥u0 N∥ 2 k∑ p=1 τ−βp Γ (2 − βp) app−1. As ak0 > ak1 > ak2 > · · · > akk−1 ≥ (1 − βk)(k)−βk , we see that I ≥ ( ∥uk N∥ 2 − ∥u0 N∥ 2) + 1 Tβk Γ (1 − βk) k∑ i=1 ∥ui N∥ 2 + ∥u0 N∥ 2 k∑ p=1 τ−βp Γ (2 − βp) app−1. □ .2. Stability analysis The weak formulation of the proposed method is as follows: find {um N } M m=1 ∈ W0 N such that( Dβm τ um N , υN ) + λ(δtum N , υN ) + B ( um N , υN ) = ( IN f (2um−1 N − um−2 N , um−Ns N ), υN ) + ( IN gm, υN ) , ∀υN ∈ W0 N , (3.13) with um N = π 1,0 N ϕm, −Ns ≤ m ≤ 0. t is a linear iterative approach, which means that we just need to find a solution to a system of linear equations at each ime level. The well-known Lax–Milgram lemma proves that the scheme is well-posed. Now, assume that { ũk N }M k=1 is the solution to( Dβm τ ũm N , υN ) + λ(δt ũm N , υN ) + B ( ũm N , υN ) = ( IN f (2ũm−1 N − ũm−2 N , ũm−Ns N ), υN ) + ( IN g̃m, υN ) , ∀υN ∈ W0 N , (3.14) with initial conditions ũm N = π 1,0 N ϕm, −Ns ≤ m ≤ 0. In the following theorem, we provide the stability with respect to the source function g . Theorem 3.1. Assume that λ > 0, 1 ≤ m ≤ M. Then, the fully discrete scheme (3.13) is unconditional stable, which means that for all τ > 0, it holds thatum N − ũm N 2 ≤ CC̃T λ m∑ s=1 ∥g s − g̃ s ∥ 2 , where C̃ = exp ( 4 λ ( 1 + 7C2L2 ) T ) is a positive constant, in which C is a generic positive constant independent on τ and N. roof. Denote ηk N = uk N − ũk N . Subtracting (3.14) from (3.13), it holds that( Dβk τ ηk N , υN ) + λ(δtηk N , υN ) + B ( ηk N , υN ) = ( IN f ( 2uk−1 N − uk−2 N , uk−Ns N ) − IN f ( 2ũk−1 N − ũk−2 N , ũk−Ns N ) , υN ) + ( IN gk − IN g̃k, υN ) . (3.15) 7 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 f A 3 Using the ϵ-Young inequality and the Minkowski inequality side by side to the Lipschitz condition on f , we derive the ollowing by the aid of Remark 3.1,( IN f (2uk−1 N − uk−2 N , uk−Ns N ) − IN f (2ũk−1 N − ũk−2 N , ũk−Ns N ), υN ) ≤ ∥IN ( f (2uk−1 N − uk−2 N , uk−Ns N ) − f (2ũk−1 N − ũk−2 N , ũk−Ns N ) ) ∥∥υN∥ ≤ ϵ 2 ∥IN ( f (2uk−1 N − uk−2 N , uk−Ns N ) − f (2ũk−1 N − ũk−2 N , ũk−Ns N ) ) ∥ 2 + 1 2ϵ ∥υN∥ 2 ≤ ϵ 2 C2L2 ( ∥2ηk−1 N − ηk−2 N ∥ + ∥η k−Ns N ∥ )2 + 1 2ϵ ∥υN∥ 2 ≤ ϵC2L2∥2ηk−1 N − ηk−2 N ∥ 2 + ϵC2L2∥ηk−Ns N ∥ 2 + 1 2ϵ ∥υN∥ 2 and ( IN gk − IN g̃k, υN ) ≤ ϵ 2 C ∥gk − g̃k ∥ 2 + 1 2ϵ ∥υN∥ 2. Then, (3.15) becomes (Dβk τ ηk N , υN ) + λ(δtηk N , υN ) + B(ηk N , υN ) ≤ 1 ϵ ∥υN∥ 2 + 4ϵC2L2∥ηk−1 N ∥ 2 + 2ϵC2L2∥ηk−2 N ∥ 2 + ϵC2L2∥ηk−Ns N ∥ 2 + ϵC 2 ∥gk − g̃k ∥ 2, Invoking Lemma 3.4, Lemma 3.5 and taking υN = ηk N yield ( ∥ηk N∥ 2 − ∥ηk−1 N ∥ 2) + τ 1−βk λΓ (2 − βk) [ ∥ηk N∥ 2 − k−1∑ i=1 (akk−i−1 − akk−i)∥η i N∥ 2 − akk−1∥η 0 N∥ 2 ] + 2τ λ |ηk N | 2 B ≤ 2τ λϵ ∥ηk N∥ 2 + 8ϵτC2L2 λ ∥ηk−1 N ∥ 2 + 4ϵτC2L2 λ ∥ηk−2 N ∥ 2 + 2ϵτC2L2 λ ∥η k−Ns N ∥ 2 + τϵC λ ∥gk − g̃k ∥ 2. (3.16) Eliminating the positive term |ηk N | 2 B, replacing k by p in (3.16) and summing the result up for p from 1 to k give k∑ p=1 ( ∥η p N∥ 2 − ∥η p−1 N ∥ 2 ) + k∑ p=1 τ 1−βp λΓ (2 − βp) [ ∥η p N∥ 2 − p−1∑ i=1 (app−i−1 − app−i)∥η i N∥ 2 − app−1∥η 0 N∥ 2 ] ≤ ( 2τ λϵ + 14ϵτC2L2 λ ) k∑ p=1 ∥η p N∥ 2 + τϵC λ k∑ p=1 ∥gp − g̃p ∥ 2. Let ϵ = 1. By the aid of Lemma 3.6, we get that( ∥ηk N∥ 2 − ∥η0 N∥ 2) + 1 Tβk Γ (1 − βk) k∑ i=1 ∥ηi N∥ 2 + ∥η0 N∥ 2 k∑ p=1 τ−βp Γ (2 − βp) app−1 ≤ ( 2τ λ + 14τC2L2 λ ) k∑ p=1 ∥η p N∥ 2 + τC λ k∑ p=1 ∥gp − g̃p ∥ 2. (3.17) As η0 N = 0 and due to the positivity of the second term of the l.h.s of (3.17), we finally get that ∥ηk N∥ 2 ≤ 2 λ ( 1 + 7C2L2 ) τ k∑ p=1 ∥η p N∥ 2 + CT λ k∑ p=1 ∥gp − g̃p ∥ 2. An application of the Grönwall inequality (see Lemma 3.3) yields for τ ≤ 4 λ ( 1 + 7C2L2 ) that ∥ηk N∥ 2 ≤ CC̃T λ k∑ p=1 ∥gp − g̃p ∥ 2, C̃ = exp ( 4 λ ( 1 + 7C2L2 ) T ) . s a result, the scheme is unconditionally stable. □ .3. Convergence analysis In this subsection, we use error estimates to analyze the convergence of the fully discrete scheme (3.13). 8 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 b t w P g S w W a A S Theorem 3.2. Assume that λ > 0. Let um be the exact solution to (1.1) at time tm for m = −Ns, . . . ,M, and let {um N } M m=−Ns e the solution to (3.13). Assume that u ∈ C2 ( [0, T ]; L2(Ω) ) ∩ C1 ( [0, T ];HL(Ω) ) and g ∈ C ((0, T );Hr (Ω)). Then, it holds hat ∥um − um N ∥ ≤ C(N1−L + N−r + τ ) , 1 ≤ m ≤ M, here C is a suitable positive constant independent of N and τ . roof. Let um − um N = emN = (um − π 1,0 N um) + (π1,0 N um − um N ) ∆ = ẽmN + êmN . The weak formulation of Eq. (1.1) at time tm is iven by (C0D βm t um, υN ) + λ(∂tum, υN ) + B ( um, υN ) = (f (um, um−Ns ), υN ) + (gm, υN ). (3.18) ubtracting (3.13) from (3.18), then the error equation satisfies (Dβm τ êmN , υN ) + λ(δt êmN , υN ) + B(êN , υN ) ∆ = Υ m 1 + Υ m 2 + Υ m 3 + Υ m 4 , (3.19) here Υ m 1 = ( IN f (um, um−Ns ) − IN f (2um−1 N − um−2 N , um−Ns N ), υN ) , Υ m 2 = ( f (um, um−Ns ) − IN f (um, um−Ns ), υN ) , Υ m 3 = ( (Dβm τ + λδt )π 1,0 N um − (C0D βm t + λ∂t )um, υN ) , Υ m 4 = ( gm − INgm, υN ) . e next estimate the right-hand terms Υ m 1 , Υ m 2 , Υ m 3 and Υ m 4 . For the first term Υ m 1 , we have that Υ m 1 = (IN f (um, um−Ns ) − IN f (2um−1 − um−2, um−Ns ), υN ) + (IN f (2um−1 − um−2, um−Ns ) − IN f (2um−1 N − um−2 N , um−Ns N ), υN ) (3.20) ∆ = Υ m 11 + Υ m 12. Applying Taylor expansion, it holds that f (um, um−Ns ) = f (2um−1 − um−2, um−Ns ) + (um − 2um−1 + um−2) f ′ 1 (ξ, um−Ns ) = f (2um−1 − um−2, um−Ns ) + c̃uτ 2. Moreover, by means of Hölder’s inequality and Young’s inequality, we have that Υ m 11 ≤ ∥IN f (um, um−Ns ) − IN f (−um−2 + 2um−1, um−n)∥ ∥υN∥ ≤ C∥f (um, um−Ns ) − f (−um−2 + 2um−1, um−Ns )∥ ∥υN∥ ≤ ϵ 2 c̃uτ 4 + 1 2ϵ ∥υN∥ 2, (3.21) nd Υ m 12 ≤ LC(∥2em−1 N − em−2 N ∥ + ∥em−Ns N ∥)∥υN∥ ≤ LC(∥2êm−1 N − êm−2 N ∥ + ∥êm−Ns N ∥ + ∥2ẽm−1 N − ẽm−2 N ∥ + ∥ẽm−Ns N ∥)∥υN∥ ≤ 8ϵ 2 CL2∥êm−1 N ∥ 2 + 2ϵ 2 CL2∥êm−2 N ∥ 2 + ϵ 2 L2∥êm−Ns N ∥ 2 + 8ϵ 2 CL2∥ẽm−1 N ∥ 2 (3.22) + 2ϵ 2 CL2∥ẽm−2 N ∥ 2 + ϵ 2 CL2∥ẽm−Ns N ∥ 2 + 1 2ϵ ∥υN∥ 2. dditionally, owing to Lemma 3.1, we notice that ∥ẽm−1 N ∥ 2 ≤ C C1 N2−2L ∥um−1 ∥ 2 L, ∥ẽm−2 N ∥ 2 ≤ C C1 N2−2L ∥um−2 ∥ 2 L, ∥ẽm−Ns N ∥ 2 ≤ C C1 N2−2L ∥um−Ns∥ 2 L. Then, (3.22) becomes Υ m 12 ≤ 4ϵCL2∥êm−1 N ∥ 2 + ϵCL2∥êm−2 N ∥ 2 + ϵ 2 CL2∥êm−Ns N ∥ 2 + C C1 N2−2L ∥u∥2 L + 1 2ϵ ∥υN∥ 2. (3.23) ubstituting (3.21) and (3.23) into (3.20), we can derive that Υ m 1 ≤ 1 ∥υN∥ 2 + 4ϵCL2∥êm−1 N ∥ 2 + ϵCL2∥êm−2 N ∥ 2 + ϵ CL2∥êm−Ns N ∥ 2 + C N2−2L ∥u∥2 L + ϵ c̃uτ 2. (3.24) ϵ 2 C1 2 9 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 F U F F S w T w For the second term Υ m 2 , by means of Hölder’s inequality, we see that Υ m 2 ≤ ϵ 2 CN−2r ∥u∥2 r + 1 2ϵ ∥υN∥ 2. (3.25) or the third term Υ m 3 , we have that Υ m 3 = ((Dβm τ + λδt )π 1,0 N um − (C0D βm t + λ∂t )π 1,0 N um, υN ) + ((C0D βm t + λ∂t )π 1,0 N um − (C0D βm t + λ∂t )um, υN ) = (π1,0 N ((Dβm τ + λδt )um − (C0D βm t + λ∂t )um), υN ) − ((C0D βm t + λ∂t )ẽmN , υN ) (3.26) ∆ = Υ m 31 + Υ m 32. sing (2.4) and Hölder’s inequality, we obtain that Υ m 31 ≤ ϵ 2 ∥π 1,0 N ((Dβm τ + λδt )um − (C0D βm t + λ∂t )um)∥2 + 1 2ϵ ∥υN∥ 2 ≤ ϵ 2 C∥((Dβm τ + λδt )um − (C0D βm t + λ∂t )um)∥2 + 1 2ϵ ∥υN∥ 2 ≤ ϵ 2 C1,uτ 2 + 1 2ϵ ∥υN∥ 2. urthermore, according to Lemma 3.1, we have that Υ m 32 ≤ ϵ 2 CN2−2L ∥ C 0D βm t um ∥ 2 L + 1 2ϵ ∥υN∥ 2 ≤ ϵ 2 CN2−2L ∥(C0D βm t + λ∂t )u∥ 2 L + 1 2ϵ ∥υN∥ 2. Thus (3.26) becomes Υ m 3 ≤ ϵ 2 CN2−2L ∥(C0D βm t + λ∂t )u∥2 L + ϵ 2 C2,uτ 2 + 1 ϵ ∥υN∥ 2. (3.27) or the fourth term Υ m 4 , by invoking Remark 3.1, it holds that Υ m 4 ≤ ϵ 2 CN−2r ∥u∥2 r + 1 2ϵ ∥υN∥ 2. (3.28) ubstituting (3.24), (3.25), (3.27) and (3.28) into (3.19), we can infer that (Dβm τ êmN , υN ) + λ(δt êmN , υN ) + B(êN , υN ) ≤ 5 2ϵ ∥υN∥ 2 + 4ϵCL2∥êm−1 N ∥ 2 + ϵCL2∥êm−2 N ∥ 2 + ϵ 2 CL2∥êm−Ns N ∥ 2 + G̃, (3.29) here G̃ = ϵC̃N2−2L ( ∥(C0D βm t + λ∂t )u∥2 L + ∥u∥2 L ) + ϵC̃N−2r ∥u∥2 r + ϵC̃uτ 2. aking υN = êmN in (3.29) and applying (3.12), we can conclude that (Dβm τ êmN , êmN ) + λ(δt êmN , êmN ) + |êmN | 2 ≤ 5 2ϵ ∥êmN ∥ 2 + 4ϵCL2∥êm−1 N ∥ 2 + ϵCL2∥êm−2 N ∥ 2 + ϵ 2 CL2∥êm−Ns N ∥ 2 + G̃. Now, we can use the same methodology as in Theorem 3.1 to complete the proof. □ 4. Numerical verification In this section, two numerical experiments are considered to demonstrate the effectiveness and convergence orders of the scheme. Example 1. We consider the following nonlinear delay reaction–diffusion problem ∂β(t)u ∂tβ(t) (x, t) + ∂u ∂t (x, t) = ∂2u ∂x2 (x, t) − 2u(x, t) + u(x, t − 0.1) 1 + u2(x, t − 0.1) + g(x, t), x ∈ (0, 1), t ∈ (0, 1], u(0, t) = u(1, t) = 0, t ∈ (0, 1), u(x, t) = (1 + t)σ sin(π x), x ∈ Ω, t ∈ [−0.1, 0], β(t) = 2 + sin(t) 4 , (4.30) here g(x, t) is a given function such that problem (4.30) has the exact solution u(x, t) = (1 + t)σ sin(π x). 10 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 b o f o F Table 1 The rate of convergence and the errors versus τ and σ with N = 50 for Example 1. τ σ = 2 σ = 1.5 σ = 1.1 Error Order Error Order Error Order 1/100 6.585 × 10−4 – 1.882 × 10−4 – 2.509 × 10−5 – 1/200 3.209 × 10−4 1.037 2.537 × 10−5 1.025 1.237 × 10−5 1.020 1/400 1.570 × 10−4 1.031 4.566 × 10−5 1.019 6.122 × 10−6 1.015 1/800 7.715 × 10−5 1.025 2.260 × 10−5 1.014 3.035 × 10−6 1.012 1/1600 3.802 × 10−5 1.021 1.121 × 10−5 1.010 1.508 × 10−6 1.009 Fig. 1. The rate of Convergence in the spatial direction with different σ at τ = 1/1600. Fig. 2. The absolute error function with σ = 1.1, 1.5, 2, N = 50 and τ = 1/2000. According to the assumptions on the solution u in Theorem 3.2, We can testify the numerical experiments at σ = 2 ut additionally we also run some numerical experiments for different values of σ not exceeds 2 for which the conditions f Theorem 3.2 cannot be satisfied. So, we devote Table 1 to show the L2-errors and corresponding convergence orders or σ = 2, 1.5, 1.1 with N = 50. These findings support the theoretically obtained convergence order in time that was btained in Theorem 3.2. Fig. 1 shows the spatial convergence orders for various values of σ at τ = 1/1600. Additionally, ig. 2 shows the space–time error functions with σ = 1.1, 1.5, 2, N = 50 and τ = 1/2000. Example 2. We investigate the following variable-order fractional equation, where the dynamics of the solution are intriguing and exact solution is unknown ∂β(t)u (x, t) + λ ∂u (x, t) = ∂2u (x, t) + u(x, t)(1 − u(x, t))(1 + u(x, t − 1.5)), x ∈ (−15, 15), t ∈ (0, 1], (4.31) ∂tβ(t) ∂t ∂x2 11 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 w a c c 5 m t Fig. 3. Behavior of the solutions to model (4.31) with λ = 0, N = 100 and τ = 0.003. Fig. 4. Behavior of the solutions to model (4.31) with λ = 1, N = 100 and τ = 0.003. ith the initial value u(x, 0) = e−2x2 . The variable orders β(t) is given by β(t) = β(T ) + (β(0) − β(T )) ( 1 − t T − sin ( 2π ( 1 − t T )) 2π ) . For N = 100 and τ = 0.003, the behavior of the numerical solution to problem (4.31) is displayed in Fig. 3 with λ = 0 nd in Fig. 4 with λ = 1 for the three cases: (I) β(0) = 0.0 and β(1) = 0.9. (II) β(0) = 0.5 and β(1) = 0.9. (III) β(0) = 0.8 and β(1) = 0.9. For λ = 0, we see that the solution for case (I) has initial weak singularities near t = 0, whereas the solution for ases (II) and (II) is smooth near t = 0. For λ = 1, these observations do not appear, and the solution is smooth for all onsidered cases. . Conclusions and remarks A numerical framework based on a combination of L1-difference approximation and Galerkin Legendre spectral ethod is introduced to solve variable order subdiffusion equation with delay. The resulted scheme is linear despite he nonlinearity of the problem under consideration. Standard techniques based on mathematical induction and discrete 12 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 A V A R energy estimates are handled to prove the unconditional estimates of convergence and stability. For the sake of clearness, some numerical experiments are constructed to show the scheme’s efficiency by examining the convergence order in both temporal and spatial directions, the absolute error for different time-fractional orders, and the numerical solution’s physical behavior. Data availability The data used to support the findings of this study are included within the article. cknowledgments The first and the fourth authors were supported by the Nazarbayev University, Kazakhstan Program 091019CRP2120. K. an Bockstal is supported by a postdoctoral fellowship of the Research Foundation - Flanders, Belgium (106016/12P2919N). . S. Hendy wishes to acknowledge the support of the RSF, Russia grant, project 22-21-00075. eferences [1] H. Sun, A. Chang, Y. Zhang, W. Chen, A review on variable-order fractional differential equations: mathematical foundations, physical models, numerical methods and applications, Fract. Calc. Appl. Anal. 22 (1) (2019) 27–59. [2] A. Burlon, G. Alotta, M. Di Paola, G. Failla, An original perspective on variable-order fractional operators for viscoelastic materials, Meccanica 56 (4) (2021) 769–784. [3] K. Van Bockstal, M. Zaky, A. Hendy, On the existence and uniqueness of solutions to a nonlinear variable order time-fractional reaction–diffusion equation with delay, Commun. Nonlinear Sci. Numer. Simul. 115 (2022) 106755. [4] N. Kadkhoda, H. Jafari, R. Ganji, A numerical solution of variable order diffusion and wave equations, Int. J. Nonlinear Anal. Appl. 12 (1) (2021) 27–36. [5] R. Ganji, H. Jafari, A numerical approach for multi-variable orders differential equations using Jacobi polynomials, Int. J. Appl. Comput. Math. 5 (2) (2019) 1–9. [6] V. Hosseini, M. Koushki, W.-N. Zou, The meshless approach for solving 2D variable-order time-fractional advection–diffusion equation arising in anomalous transport, Eng. Comput. (2021) 1–19. [7] V. Hosseini, F. Yousefi, W.-N. Zou, The numerical solution of high dimensional variable-order time fractional diffusion equation via the singular boundary method, J. Adv. Res. 32 (2021) 73–84. [8] A. Babaei, H. Jafari, S. Banihashemi, Numerical solution of variable order fractional nonlinear quadratic integro-differential equations based on the sixth-kind Chebyshev collocation method, J. Comput. Appl. Math. 377 (2020) 112908. [9] R. Ganji, H. Jafari, S. Nemati, A new approach for solving integro-differential equations of variable order, J. Comput. Appl. Math. 379 (2020) 112946. [10] L. Ramirez, C. Coimbra, A variable order constitutive relation for viscoelasticity, Ann. Phys. 16 (7–8) (2007) 543–552. [11] G. Scarpi, Sulla possibilita di un modello reologico intermedio di tipo evolutivo, Atti Accad. Naz. Lincei. Cl. Sci. Fis. Mat. Natur. Rend. 52 (6) (1972) 912–917. [12] R. Garrappa, A. Giusti, F. Mainardi, Variable-order fractional calculus: A change of perspective, Commun. Nonlinear Sci. Numer. Simul. 102 (2021) 105904. [13] S. Patnaik, J. Hollkamp, F. Semperlotti, Applications of variable-order fractional operators: a review, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 476 (2234) (2020) 20190498. [14] S. Samko, Fractional integration and differentiation of variable order: an overview, Nonlinear Dynam. 71 (4) (2013) 653–662. [15] M. Ortigueira, D. Valério, J. Machado, Variable order fractional systems, Commun. Nonlinear Sci. Numer. Simul. 71 (2019) 231–243. [16] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339 (1) (2000) 1–77. [17] Y. Luchko, Initial-boundary-value problems for the one-dimensional time-fractional diffusion equation, Fract. Calc. Appl. Anal. 15 (1) (2012) 141–160. [18] N. Kopteva, Error analysis of the L1 method on graded and uniform meshes for a fractional-derivative problem in two and three dimensions, Math. Comp. 88 (319) (2019) 2135–2155. [19] D. Benson, R. Schumer, M. Meerschaert, S. Wheatcraft, Fractional dispersion, Lévy motion, and the MADE tracer tests, Transp. Porous Media 42 (1) (2001) 211–240. [20] M. Abbaszadeh, M. Dehghan, Meshless upwind local radial basis function-finite difference technique to simulate the time-fractional distributed-order advection–diffusion equation, Eng. Comput. 2021 (37) (2012) 873–889. [21] M. Abbaszadeh, M. Dehghan, A finite-difference procedure to solve weakly singular integro partial differential equation with space-time fractional derivatives, Eng. Comput. 37 (2021) (2021) 2173–2182. [22] A. Hendy, M. Zaky, D. Suragan, Discrete fractional stochastic Grönwall inequalities arising in the numerical analysis of multi-term fractional order stochastic differential equations, Math. Comput. Simulation 193 (2022) 269–279. [23] A. Hendy, M. Zaky, R. Hafez, R. De Staelen, The impact of memory effect on space fractional strong quantum couplers with tunable decay behavior and its numerical simulation, Sci. Rep. 11 (1) (2021) 1–15. [24] A. Bhrawy, M. Zaky, An improved collocation method for multi-dimensional space–time variable-order fractional Schrödinger equations, Appl. Numer. Math. 111 (2017) 197–218. [25] A. Bhrawy, M. Zaky, Highly accurate numerical schemes for multi-dimensional space variable-order fractional Schrödinger equations, Comput. Math. Appl. 73 (6) (2017) 1100–1117. [26] B. Baeumer, M. Kovács, M.M. Meerschaert, H. Sankaranarayanan, Boundary conditions for fractional diffusion, J. Comput. Appl. Math. 336 (2018) 408–424. [27] J. Jia, H. Wang, X. Zheng, A preconditioned fast finite element approximation to variable-order time-fractional diffusion equations in multiple space dimensions, Appl. Numer. Math. 163 (2021) 15–29. [28] H. Wang, X. Zheng, Wellposedness and regularity of the variable-order time-fractional diffusion equations, J. Math. Anal. Appl. 475 (2) (2019) 1778–1802. [29] D. Tavares, R. Almeida, D. Torres, Caputo derivatives of fractional variable order: numerical approximations, Commun. Nonlinear Sci. Numer. Simul. 35 (2016) 69–87. 13 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb1 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb1 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb1 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb2 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb2 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb2 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb3 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb3 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb3 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb4 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb4 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb4 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb5 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb5 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb5 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb6 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb6 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb6 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb7 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb7 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb7 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb8 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb8 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb8 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb9 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb9 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb9 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb10 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb11 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb11 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb11 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb12 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb12 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb12 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb13 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb13 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb13 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb14 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb15 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb16 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb17 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb17 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb17 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb18 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb18 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb18 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb19 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb19 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb19 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb20 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb20 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb20 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb21 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb21 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb21 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb22 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb22 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb22 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb23 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb23 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb23 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb24 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb24 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb24 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb25 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb25 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb25 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb26 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb26 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb26 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb27 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb27 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb27 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb28 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb28 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb28 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb29 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb29 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb29 M.A. Zaky, K. Van Bocksta, T.R. Taha et al. Journal of Computational and Applied Mathematics 420 (2023) 114832 [30] S. Chen, F. Liu, K. Burrage, Numerical simulation of a new two-dimensional variable-order fractional percolation equation in non-homogeneous porous media, Comput. Math. Appl. 68 (12) (2014) 2133–2141. [31] F. Zeng, Z. Zhang, G. Karniadakis, A generalized spectral collocation method with tunable accuracy for variable-order fractional differential equations, SIAM J. Sci. Comput. 37 (6) (2015) A2710–A2732. [32] M. Zayernouri, G. Karniadakis, Fractional spectral collocation methods for linear and nonlinear variable order FPDEs, J. Comput. Phys. 293 (2015) 312–338. [33] T. Zhao, Z. Mao, G. Karniadakis, Multi-domain spectral collocation method for variable-order nonlinear fractional differential equations, Comput. Methods Appl. Mech. Engrg. 348 (2019) 377–395. [34] X. Zheng, H. Wang, Optimal-order error estimates of finite element approximations to variable-order time-fractional diffusion equations without regularity assumptions of the true solutions, IMA J. Numer. Anal. 41 (2) (2021) 1522–1545. [35] X. Zheng, H. Wang, A hidden-memory variable-order time-fractional optimal control model: Analysis and approximation, SIAM J. Control Optim. 59 (3) (2021) 1851–1880. [36] J. Jia, H. Wang, X. Zheng, A fast collocation approximation to a two-sided variable-order space-fractional diffusion equation and its analysis, J. Comput. Appl. Math. 388 (2021) 113234. [37] H.-K. Pang, H.-W. Sun, A fast algorithm for the variable-order spatial fractional advection-diffusion equation, J. Sci. Comput. 87 (1) (2021) 1–28. [38] L. Wei, Y. Yang, Optimal order finite difference/local discontinuous Galerkin method for variable-order time-fractional diffusion equation, J. Comput. Appl. Math. 383 (2021) 113129. [39] L. Wei, S. Zhai, X. Zhang, Error estimate of a fully discrete local discontinuous Galerkin method for variable-order time-fractional diffusion equations, Commun. Appl. Math. Comput. (2020) 1–15. [40] H. Liu, A. Cheng, H. Wang, A parareal finite volume method for variable-order time-fractional diffusion equations, J. Sci. Comput. 85 (1) (2020) 1–27. [41] X.-M. Gu, H.-W. Sun, Y.-L. Zhao, X. Zheng, An implicit difference scheme for time-fractional diffusion equations with a time-invariant type variable order, Appl. Math. Lett. 120 (2021) 107270. [42] K. Van Bockstal, Existence of a unique weak solution to a non-autonomous time-fractional diffusion equation with space-dependent variable order, Adv. Difference Equ. 2021 (1) (2021) 1–43. [43] V. Pimenov, A. Hendy, Numerical studies for fractional functional differential equations with delay based on BDF-type shifted Chebyshev approximations, in: Abstract and Applied Analysis, Vol. 2015, Hindawi, 2015. [44] M. Zaky, A. Hendy, J. Macías-Díaz, Semi-implicit Galerkin–Legendre spectral schemes for nonlinear time-space fractional diffusion–reaction equations with smooth and nonsmooth solutions, J. Sci. Comput. 82 (1) (2020) 1–27. [45] A. Hendy, M. Zaky, Combined Galerkin spectral/finite difference method over graded meshes for the generalized nonlinear fractional Schrödinger equation, Nonlinear Dynam. 103 (3) (2021) 2493–2507. [46] C. Lorenzo, T. Hartley, Variable order and distributed order fractional operators, Nonlinear Dynam. 29 (1) (2002) 57–98. [47] P. Zhuang, F. Liu, V. Anh, I. Turner, Numerical methods for the variable-order fractional advection-diffusion equation with a nonlinear source term, SIAM J. Numer. Anal. 47 (3) (2009) 1760–1781. [48] X. Zheng, H. Wang, An error estimate of a numerical approximation to a hidden-memory variable-order space-time fractional diffusion equation, SIAM J. Numer. Anal. 58 (5) (2020) 2492–2514. [49] X. Zheng, H. Wang, A time-fractional diffusion equation with space-time dependent hidden-memory variable order: analysis and approximation, BIT Numer. Math. (2021) 1–29. [50] H. Wang, X. Zheng, Analysis and numerical solution of a nonlinear variable-order fractional differential equation, Adv. Comput. Math. 45 (5) (2019) 2647–2675. [51] X. Zheng, Z. Zhang, H. Wang, Analysis of a nonlinear variable-order fractional stochastic differential equation, Appl. Math. Lett. 107 (2020) 106461. [52] Q. Li, G. Wang, M. Wei, Monotone iterative technique for time-space fractional diffusion equations involving delay, Nonlinear Anal. Model. Control 26 (2) (2021) 241–258. [53] Z. Liu, S. Zeng, Y. Bai, Maximum principles for multi-term space-time variable-order fractional diffusion equations and their applications, Fract. Calc. Appl. Anal. 19 (1) (2016) 188–211. [54] S. Shen, F. Liu, J. Chen, I. Turner, V. Anh, Numerical techniques for the variable order time fractional diffusion equation, Appl. Math. Comput. 218 (22) (2012) 10861–10870. [55] A. Alikhanov, Boundary value problems for the diffusion equation of the variable order in differential and difference settings, Appl. Math. Comput. 219 (8) (2012) 3938–3946. [56] F. Zeng, F. Liu, C. Li, K.n. Burrage, I. Turner, V. Anh, A Crank–Nicolson ADI spectral method for a two-dimensional Riesz space fractional nonlinear reaction-diffusion equation, SIAM J. Numer. Anal. 52 (6) (2014) 2599–2622. [57] J. Shen, Efficient spectral-Galerkin method I. Direct solvers of second-and fourth-order equations using Legendre polynomials, SIAM J. Sci. Comput. 15 (6) (1994) 1489–1505. [58] J. Shen, T. Tang, L.-L. Wang, Spectral Methods: Algorithms, Analysis and Applications, Vol. 41, Springer Science & Business Media, 2011. [59] Y. Lin, C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys. 225 (2) (2007) 1533–1552. [60] M. Zaky, A. Hendy, A. Alikhanov, V. Pimenov, Numerical analysis of multi-term time-fractional nonlinear subdiffusion equations with time delay: What could possibly go wrong? Commun. Nonlinear Sci. Numer. Simul. 96 (2021) 105672. 14 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb30 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb30 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb30 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb31 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb31 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb31 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb32 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb32 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb32 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb33 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb33 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb33 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb34 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb34 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb34 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb35 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb35 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb35 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb36 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb36 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb36 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb37 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb38 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb38 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb38 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb39 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb39 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb39 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb40 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb40 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb40 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb41 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb41 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb41 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb42 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb42 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb42 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb43 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb43 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb43 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb44 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb44 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb44 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb45 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb45 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb45 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb46 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb47 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb47 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb47 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb48 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb48 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb48 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb49 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb49 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb49 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb50 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb50 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb50 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb51 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb51 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb51 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb52 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb52 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb52 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb53 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb53 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb53 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb54 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb54 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb54 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb55 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb55 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb55 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb56 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb56 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb56 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb57 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb57 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb57 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb58 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb59 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb60 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb60 http://refhub.elsevier.com/S0377-0427(22)00430-7/sb60 An L1 type difference/Galerkin spectral scheme for variable-order time-fractional nonlinear diffusion–reaction equations with fixed delay Introduction The linearized numerical scheme Theoretical analysis of the fully discrete scheme Technical lemmas Stability analysis Convergence analysis Numerical verification Conclusions and remarks Data availability Acknowledgments References