Compact difference scheme for semi-linear delayed diffusion
      wave equation with fractional order in time

                           Ahmed S. Hendy1,2                           Vladimir G. Pimenov1,3
                        ahmed.hendy@fsc.bu.edu.eg                       v.g.pimenov@urfu.ru
      1 – Institute of Natural Sciences and Mathematics, Ural Federal University (Yekaterinburg, Russia)
                       2 – Department of Mathematics, Benha University (Benha, Egypt)
                 3 – Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg, Russia)




                                                        Abstract
                       A compact difference scheme for a class of non-linear fractional
                       diffusion-wave equations with fixed time delay is considered. Analysis
                       of the constructed difference scheme is done in L∞ -norm by means of
                       the discrete energy method. A numerical test example is introduced to
                       illustrate the accuracy and efficiency of the proposed method.




1    Introduction
An important class of fractional differential equations which has been studied widely in recent years is the
time fractional diffusion-wave equation (FDWE). The time FDWE is obtained from the classical diffusion-wave
equation by replacing the second-order time derivative term by a fractional derivative of order 1 < α ≤ 2.
The fractional diffusion equation was introduced in physics by Nigmatullin [16] to describe diffusion in media
with fractal geometry, which is a special type of porous media. Gorenflo et al. [10] presented the scale-invariant
solutions for the time-fractional diffusion-wave equation in terms of the generalized Wright function. Agrawal [2]
extended this formulation to a diffusion-wave equation that contains a fourth-order space derivative term in a
bounded space domain. Recently, simulations of the approximation solutions of time-fractional wave, forced wave
(shear wave) and damped wave equations are given in [1]. A novel fractional diffusion-wave equation with non-
local regularization for noise removal was presented in [24]. The existence and uniqueness of solutions for Dirichlet
initial-boundary value problem associated to the semi linear fractional wave equation was recently studied in [14].
As a numerical approach to solve FDWE, Sun and his co-authors proposed a high order difference methods for the
fractional diffusion-wave equation [8]. Also, in [22], the efforts of the authors were devoted to the application of
fractional multi-step method to obtain a numerical solution of time fractional diffusion-wave equation. Compact
finite difference schemes for the modified anomalous fractional sub-diffusion equation and fractional diffusion-
wave equation were studied in [21]. A numerical solution for a general class of diffusion problem was considered
in [3], where the standard time derivative is replaced by a fractional one. An efficient numerical method was
constructed in [26] to solve this moving boundary problem.
   Time delay occurs in many realistic applications which are modeled mathematically, e.g. [9, 6, 13, 4, 15, 5].
Reaction-diffusion equations with time delay effect have been proposed as models for the population ecology, the
cell biology and the control theory in recent years [19]. The existence of mild solutions for initial value problem
for nonlinear time fractional non-autonomous evolution equations with delay in Banach space E was studied
in [7]. Numerically, a linearized quasi-compact difference scheme was proposed for semi-linear space-fractional

Copyright c by the paper’s authors. Copying permitted for private and academic purposes.
In: A.A. Makhnev, S.F. Pravdin (eds.): Proceedings of the International Youth School-conference «SoProMat-2017», Yekaterinburg,
Russia, 06-Feb-2017, published at http://ceur-ws.org




                                                             325
diffusion equations with a fixed time delay [11]. A linearized compact finite difference scheme was presented for
the semi-linear fractional delay convection-reaction–diffusion equation in [23]. The authors of the manuscript at
hand, recently proposed a difference scheme for a class of non-linear delay distributed order fractional diffusion
equations in [17]. As an extension to this contribution and depending on Sun ’s work [20], we seek to derive a
compact linear difference scheme to solve numerically FDWE effected with a non-linear delayed source function,
more specific we consider

                       ∂ α u(x, t)    ∂ 2 u(x, t)
                             α
                                   =K             + f (x, t, u(x, t), u(x, t − s)), t > 0, 0 ≤ x ≤ L,                  (1a)
                           ∂t             ∂x2
with the following initial and boundary conditions

                                                                               ∂u(x, 0)               ∂ r̃(x, t)
                   u(x, t) = r̃(x, t),    0 ≤ x ≤ L,          t ∈ [−s, 0],              = ψ̃(x) = lim            ,     (1b)
                                                                                 ∂t              t→−0     ∂t
                   u(0, t) = φ0 (t),     u(L, t) = φL (t),         t > 0,                                              (1c)

where s > 0 is the delay parameter, K is a positive constant. The fractional derivative of order 1 < α ≤ 2 is
defined in Caputo sense.
   In order to transform (1) to a system with zero Dirichlet boundary conditions, we define h(x, t) := φ0 (t) +
x
L (φL (t) − φ0 (t)) and introduce the new function v(x, t) = u(x, t) − h(x, t). Hence, we have

                        ∂ α v(x, t)     ∂ 2 v(x, t)
                                    = K             + f (x, t, v(x, t), v(x, t − s)), t > 0, 0 ≤ x ≤ L,                (2a)
                            ∂tα             ∂x2
with the following initial and boundary conditions

                                                                               ∂v(x, 0)              ∂r(x, t)
                   v(x, t) = r(x, t),     0 ≤ x ≤ L,          t ∈ [−s, 0],              = ψ(x) = lim          ,        (2b)
                                                                                 ∂t             t→−0   ∂t
                   v(0, t) = v(L, t) = 0,      t > 0.                                                                  (2c)

   We need to overcome two degrees of complexity; how to employ a suitable approximation of the time fractional
derivative on the one hand, and how to approximate the non linear delay source function linearly on the other,
in order to obtain a numerical solution for (2). Throughout this work and by the aid of [23], we suppose that the
function f (x, t, µ, v) and the solution u(x, t) of (2) are sufficiently smooth in the following sense:
    • Let m be an integer
                        Smsatisfying ms ≤ T < (m + 1)s, define Ir = (rs, (r + 1)s), for r = −1, 0, . . . , m − 1,
      Im = (ms, T ), I = q=−1 Iq and assume that u(x, t) ∈ C (6,3) ([0, L] × [0, T ]),
    • The partial derivatives fµ (x, t, µ, v) and fv (x, t, µ, v) are continuous in the 0 -neighborhood of the solution.
      Define

                                  c1 =         sup            |fµ (x, t, u(x, t) + 1 , u(x, t − s) + 2 )| ,          (3a)
                                         0<x<L, 0<t≤T
                                         |1 |≤0 ,|2 |≤0

                                  c2 =         max            |fv (x, t, u(x, t) + 1 , u(x, t − s) + 2 )| .          (3b)
                                         0<x<L, 0<t≤T
                                         |1 |≤0 ,|2 |≤0


The structure of this paper is arranged as: a derivation of the linear difference scheme is done in the following
section. Next, in the third section, the solvability, convergence and stability for the difference scheme are carried
out. In the fourth section, numerical examples are given to illustrate the accuracy of the presented scheme and
to support our theoretical results.

2     Construction of the difference scheme
A numerical solution based on the Crank-Nicholson method is derived. Before we continue, some further notations
                                                          L
are fixed. Take two positive integers M and n, let h = M    , τ = ns and denote xi = i h for i = 0, . . . , M ; tk = k τ and
tk−1/2 = k − 2 τ = 2 (tk +tk−1 ), for k = −n, . . . , N , where N = Tτ . Using the points xi in space and tk in time
               1      1
                 
                                                      

we cover the space-time domain by Ωhτ = Ωh × Ωτ , where Ωh = {xi | 0 ≤ i ≤ M } and Ωτ = {tk | −n ≤ k ≤ N }.




                                                                     326
          
Let W = w : Ωhτ → R | w(xi , tk ) = Wik ; i = 0, 1, . . . , M ; k = −n, −n + 1, . . . , N        be a grid function space on
                            k−1/2
                                   := 21 Wik + Wik−1 .
                                                           

Ωhτ . For w ∈ W, we define Wi
  Lemma 1. Let q(x) ∈ C 6 ([xi−1 , xi+1 ]), then

               1 00                                          1                                   h4 (6)
                 (q (xi−1 ) + 10q 00 (xi ) + q 00 (xi+1 )) − 2 (q(xi−1 ) − 2q(xi ) + q(xi+1 )) =     q (ωi ),
              12                                            h                                    240
where ωi ∈ (xi−1 , xi+1 ).[25]
  In [20], an approximation for the time Caputo fractional derivative at tk−1/2 with 1 < α < 2 was given:
                                                                                                  
                                                       k−1
              ∂ α v(xi , tk−1/2 )  1           k−1/2
                                                       X                        j−1/2
                                  = bα 0 δt Vi      −     (bα        α
                                                             k−j−1 − bk−j )δt Vi      − bα
                                                                                         k−1 ψ(xi )
                                                                                                    + rk ,                   (4a)
                     ∂tα           µ̄                  j=1


where ψ(x) is defined in (2b),

                                      τ 2−α
                                  bα         (k + 1)2−α − k 2−α , µ̄ = τ Γ(2 − α),
                                                               

                                   k =                                                                                       (4b)
                                      2−α

                                                 23−α                             ∂ 3 v 3−α
                                                                       
                                 1       2−α                          1
                     |rk | ≤                  +       − (1 + 21−α ) +      max          τ   ,                                 (4c)
                             Γ(3 − α)     12     3−α                  12 0≤t≤tk ∂t3

and for any function v : [0, L] × [−s, +∞) → R one denotes v(xi , tj ) = Vij for i ∈ N, j ∈ Z and defines

                              k−1/2       1                           k−1/2       1
                                            Vik + Vik−1 ,                           Vik − Vik−1 ,
                                                       
                                       

                            Vi        =                          δt Vi        =                                              (4d)
                                          2                                       τ

                               k          1                               1
                                            Vik − Vi−1
                                                    k
                                                              δx2 Vik =         k
                                                                                   − 2Vik + Vi−1
                                                                                              k
                                                       
                                         

                          δx Vi−1/2 =                    ,                  2
                                                                              Vi+1                 .                          (4e)
                                          h                               h
We are now in a position to apply and combine the above, that is (4), to (2a) at the points (xi , tk−1/2 ), and
arrive at
                                                                     
                       k−1
   1
   bα        k−1/2
                       X                      j−1/2
                           bα        α
                                                    − bα
                                        

        0 δt Vi      −      k−j−1 − bk−j δt Vi         k−1 ψ(xi )
                                                                  + rk 
   µ̄                  j=1

                                                   ∂ 2 v(xi , tk−1/2 )
                                              =K                       + f (xi , tk−1/2 , v(xi , tk−1/2 ), v(xi , tk−1/2 − s)), (5)
                                                          ∂x2
such that i = 0, . . . , M, k = 1, . . . , N.
   Lemma 2. For g = (g0 , g1 , . . . , gM ), let the linear operator A be defined as

                                             1
                                   Agi =       (gi−1 + 10gi + gi+1 ),       1 ≤ i ≤ M − 1.
                                            12
Then, we obtain
                                                                    
                         k−1
     1           k−1/2
                         X                      j−1/2
  A  bα                    bα          α
                                                       − bα
                                            

          0 δt Vi      −      k−j−1 − bk−j δt Vi           k−1 ψ(xi )
                                                                      
     µ̄                  j=1
                                                                                          
                                   k−1/2                   3           1   1         1         k−1/2
                          = Kδx2 Vi      + Af xi , tk−1/2 , Vik−1 + Vik−2 , Vik−n−1 + Vik−n + Ri     , (6)
                                                           2           2   2         2

where
                                 k−1/2
                                          ≤ C̄ τ 3−α + h4 ,
                                                         

                             Ri                                 1 ≤ i ≤ M − 1,        1 ≤ k ≤ N.                               (7)




                                                              327
  Proof. By using the following Taylor expansions

                           ∂ 2 v(xi , tk−1/2 )   1 ∂ 2 v(xi , tk ) ∂ 2 v(xi , tk−1 )
                                                                                    
                                                                                       + O τ2 ,
                                                                                             

                                       2
                                               =          2
                                                                  +          2
                                   ∂x            2      ∂x               ∂x
                                                 3 k−1 1 k−2
                                                                    + O τ2 ,
                                                                              

                               v(xi , tk−1/2 ) = Vi     − Vi
                                                 2         2
                                                 1 k−n 1 k−n−1
                                                                       + O τ2 ,
                                                                                 

                          v(xi , tk−1/2 − s) = Vi       + Vi
                                                 2         2

in (5) we obtain

                                                                          
                         k−1
   1
   bα        k−1/2
                         X                         j−1/2
        0 δt Vi      −       (bα          α
                               k−j−1 − bk−j )δt Vi       − bα
                                                            k−1 ψ(xi )
                                                                       + rk  =
   µ̄                    j=1

        K  ∂ 2 v(xi , tk ) ∂ 2 v(xi , tk−1 ) 
                                                                                          
                                                                 3 k−1 1 k−2 1 k−n 1 k−n−1
                           +                    + f xi , tk−1/2 , Vi    − Vi , Vi  + Vi      + O(τ 2 ), (8)
        2        ∂x2              ∂x2                            2        2      2  2

where we used the continuity of the derivatives of f in its third and fourth component when letting τ → 0.
According to Lemma 1, we have

                                 ∂ 2 v(xi , tk )             h4 ∂ 6 v k 

                                                
                                                       2 k
                         A                         = δ  V
                                                       x i +          θ , tk ,    θik ∈ (xi−1 , xi+1 ),
                                      ∂x2                    240 ∂x6 i

so applying A to (8) we arrive at, using (4e),

                                                                        
                         k−1
     1           k−1/2
                         X                      j−1/2
  A  bα                    bα        α
                                                      − bα           + rk  + O h4 + O τ 2
                                          
                                        
        

          0 δt Vi      −      k−j−1 − bk−j δt Vi         k−1 ψ(xi )
     µ̄                  j=1
                                                                                                 
                                 K                                     3        1     1       1
                              = δx2 Vik + Vik−1 + Af xi , tk−1/2 , Vik−1 − Vik−2 , Vik−n + Vik−n−1 ,
                                                  

                                 2                                     2        2     2       2

                                                 k−1/2             
       

as v(x, t) ∈ C (6,3) ([0, L] × [0, T ]). Define Ri     = Ark + O h4 + O τ 2 , then from (4c), the estimate (7) is
achieved and the proof is complete.                                                                             
                                                                       k−1/2
   The final form of our difference scheme is obtained by neglecting Ri       and replacing Vik with vik in (6)

                                                               
                      k−1
     1    α   k−1/2
                      X
                            α        α        j−1/2    α
  A  b0 δt vi     −     (bk−j−1 − bk−j )δt vi     − bk−1 ψ(xi )
     µ̄               j=1
                                                                                             
                                         2 k−1/2                    3 k−1 1 k−2 1 k−n−1 1 k−n
                                   = Kδx vi       + Af xi , tk−1/2 , vi − vi , vi      + vi     ,               (9a)
                                                                    2     2     2       2

such that 1 ≤ i ≤ M − 1,         1 ≤ k ≤ N, and supplying appropriate initial and boundary conditions

                                        v0k = φ0 (tk ),       k
                                                             vM = φL (tk ),   1 ≤ k ≤ N,                        (9b)
                                        vik = r(xi , tk ),    0 ≤ i ≤ M,      −n ≤ k ≤ 0.                       (9c)

Recall that v(xi , tk ) = Vik and vik is the solution of the difference scheme, hoping to have v(xi , tk ) ≈ vik , we
discuss in the next section ki := |Vik − vik |.
   We now prove that our difference scheme admits a unique solution. Next, we show that the obtained solution
solves (2).
   Theorem 1. (Solvability). The difference scheme (9) is uniquely solvable.




                                                                  328
    Proof. We can arrange the system (9) as follows

         1 bα                10 bα                     1 bα
                                                                   
             0   K   k            0      K                  0      K
               − 2 vi+1 +             + 2 vik +                − 2 vi−1     k
        12 µ̄τ  2h           12 µ̄τ     h             12 µ̄τ      2h
                                                                                                  
                                1                      K        k−1        10                     K k−1
                         =        µ̄τ (2bα
                                         0 − b α
                                               1 ) −          v      +        µ̄τ (2b α
                                                                                      0 − b α
                                                                                            1 ) +      v
                               12                     2h2 i+1              12                     h2 i
                                                                                 
                                                   1                          K
                                             +        µ̄τ (2bα0  − b α
                                                                     1 ) −          v k−1
                                                   12                       2h2 i−1
                                                                                                            
                                  k−2
                                  X                              k−1
                                                                 X
                                                            j                                j−1
                       + µ̄τ A  (bα                α
                                        k−j−1 − bk−j )vi −            (bα             α
                                                                         k−j−1 − bk−j )vi        + τ bα
                                                                                                      k−1 ψi
                                                                                                             
                                       j=1                             j=1
                                                                                                                         
                                                                                          3       1       1         1
                                                                 + Af        xi , tk+1/2 , vik−1 − vik−2 , vik+1−n + vik−n ,   (10)
                                                                                          2       2       2         2

or written in a more concise form
                                                    Av k = ϕk (v k−1 , v k−1 , . . . , v −n ).
                                                                                                             P
The tridiagonal coefficient matrix A = (aij ) is strictly diagonally dominant because aii >                         |aij |;
                                                                                                             j6=i

                                              10 bα0  K                         1 bα0   K
                                      aii =          + 2,         ai+1,i =            − 2 = ai−1,i .
                                              12 µ̄τ  h                        12 µ̄τ  2h
Therefore, the coefficient matrix is nonsingular and the theorem is readily proved by strong induction.                          

3       Convergence and stability for the difference scheme
Now, we introduce the uniqueness, stability and convergence theorems in L∞ -norm using the discrete energy
method for the proposed difference scheme.
The spatial domain [0, L] is covered by Ωh = {xi | 0 ≤ i ≤ M, } and let Vh = {v | v = (v0 , . . . , vM ), v0 = vM =
0} be a grid function space on Ωh . For any u, v ∈ Vh , define the discrete inner products and corresponding norms
as
                                     M
                                     X −1                          M
                                                                   X
                          hu, vi = h      ui vi , hδx u, δx vi = h   (δx ui−1/2 )(δx vi−1/2 ),
                                              i=1                                i=1

                                                                          1
                                       hδx2 u, vi = −hδx u, δx vi, δx ui = (ui − ui−1 ),
                                                                          h
                                        p                    p
                                  kuk = hu, ui, |u|1 = hδx u, δx ui, kuk∞ = max |u|,
                                                                                                     0≤i≤M

and denote                                    v                                  v
                                              u M                                u M −1
                                              u X                                u X
                                                                         2
                                       |v|1 = th  (δx vi )2 ,          kδx vk∞ = th    (δx2 vi )2 .
                                                     i=1                                    i=1

According to [25], the following inequalities are fulfilled
                                                  √
                                                    L                              L
                                       kuk∞ ≤         |u|1 ,                 kuk ≤ √ |u|1 .                                    (11)
                                                   2                                6
For the analysis of the difference scheme, we will use the following lemmas:
  Lemma 3. Let v ∈ Vh , we have kδx2 vk ≤ h2 |v|1 .
                                                                                                 √
    Lemma 4. [18] Let v ∈ Vh and v0 = vM = 0. Then, we have kvk∞ ≤ 2L |v|1 .
    Lemma 5. [20] For any G = {G1 , G2 , G3 , . . .} and ψ, we obtain
            m
                "         k−1
                                                            #                    m
                                                                      1−αl                  2−αl
           X
                   α
                          X α            αl            αl            tm          X         tm
                  b0 Gk −     (bk−j−1 − bk−j )Gj − bk−1 ψ Gk ≥
                                 l
                                                                               τ   G2k −             ψ2 .
                          j=1
                                                                   2(2 − α l )           2(2 − α l )
                   k=1                                                                               k=1




                                                                      329
   Lemma 6. Gronwall inequality [25]. Suppose that {F k | k ≥ 0} is a nonnegative sequence and satisfies
    k+1
               Pk
F     ≤ A + Bτ l=1 F l , k ≥ 0, for some nonnegative constants A and B. Then F k+1 ≤ A exp(Bkτ ).
   Theorem 2. (Convergence). Let v(x, t) ∈ [0, L]×[−s, T ], be the solution of (2) such that v(xi , tk ) = Vik and
vi (0 ≤ i ≤ M, −n ≤ k ≤ N ) be the solution of the difference scheme (9), denote eki = Vik − vik , for 0 ≤ i ≤ M ,
 k

−n ≤ k ≤ N ,                  r
                                   3L       3T                     1
                          C=           exp       (5c21 + c22 ) , θ = T α−1 Γ(3 − α),
                                  8θK        8θK                     2
then if
                                                   1
                                               3−α                  41
                                               0                       0
                                   τ ≤ τ0 =            , h ≤ h0 =            ,                                (12)
                                              4C                      4C
one has that
                                    k ek k∞ ≤ C τ 3−α + h4 , 0 ≤ k ≤ N,
                                                              

                                                                                                              (13)
where c1 , c2 and 0 are those from (3).
  The proof uses the previous formulated lemmas in the sense of our results in [12].

   To discuss the stability of the difference scheme (9), we also use the discrete energy method in the same way
like the discussion of the convergence. Let {zik | 0 ≤ i ≤ M, 0 ≤ k ≤ N } be the solution of
                                                                    
                          k−1
       1          k−1/2
                          X                       j−1/2
   A  bα  0 δt zi     −     (bα        α
                                k−j−1 − bk−j )δt zi     − bα
                                                           k−1 ψ(xi )
                                                                      
       µ̄                 j=1
                                                                                                       
                                                  k−1/2                   1     1       1         1
                                        = Kδx2 zi       + Af xi , tk−1/2 , zik + zik−1 , zik−n−1 + zik−n , (14a)
                                                                          2     2       2         2
such that 1 ≤ i ≤ M − 1,     1 ≤ k ≤ N, and supplied by appropriate initial and boundary conditions
                                  z0k = φ0 (tk ),     k
                                                     zM = φL (tk ),         1 ≤ k ≤ N,                            (14b)
                                  zik = r(xi , tk ) + ρki ,       0 ≤ i ≤ M,     −n ≤ k ≤ 0,                      (14c)

where ρki is the perturbation of ψ(xi , tk ).
Following the same steps as in the proof of the convergence theorem, the stability of the scheme is obtained.
Theorem 3. (Stability). Let θik = zik − vik , for 0 ≤ i ≤ M, −n ≤ k ≤ N. Then there exist some arbitrary
positive constants c4 , c5 , h0 , τ0 which fulfill
                                                                           v
                                                0                          u M −1
                              k
                                         √    X
                                                     k                k
                                                                           u X
                            kθ k∞ ≤ c4 τ           kρ k, 0 ≤ k ≤ N, kρ k = th     (ρki )2 ,
                                           k=−n                                                 i=1

conditioned by
                                         h ≤ h0 ,     τ ≤ τ0 ,          max |ρki | ≤ c5 .
                                                                      −n≤k≤0
                                                                      0≤i≤M



4     Numerical Verification
Let νik be the solution of the constructed difference scheme (9) with the step sizes τ and h. Define the maximum
norm error by E(τ, h) = max kVik − νik k∞ . Also, define the following error rates
                         0≤i≤M
                         0≤k≤N
                                                                                                 
                                                  E(2τ, h)                               E(τ, 2h)
                              rate1 = log2                        ,   rate2 = log2                      .
                                                  E(τ, h)                                E(τ, h)
Consider the following numerical test example
                   ∂ α u(x, t)   ∂ 2 u(x, t)
                               =             + f (x, t, u(x, t), u(x, t − s)),           t ∈ (0, 1), 0 < x < π,   (15a)
                       ∂tα           ∂x2




                                                                  330
                                                                  
                                                         Γ(4)
f (x, t, u(x, t), u(x, t − s)) = sin(x) (t3 + 2t + 4) + Γ(4−α) t3−α − u(x, t − s) + sin(x)((t − s)3 + 2(t − s) + 4),
with the following initial and boundary conditions

                         u(x, t) = (t3 + 2t + 4) sin(x),   0 ≤ x ≤ π,    t ∈ [−s, 0),   s > 0,                   (15b)

                                          u(0, t) = u(π, t) = 0,   t ∈ [0, 1].                                   (15c)
The exact solution of this problem is
                                            u(x, t) = (t3 + 2t + 4) sin(x).                                        (16)
Results are presented in Tables (1) - (2) (time) and Table 3 (space). From these numerical results, we can see a
good agreement between theoretical and numerical results.

Table 1: Errors and convergence order of the difference scheme (9) for (15) in time variable with α = 1.25,
h = π/3000 and with time delay s = 1.

                                              τ       E(τ, h)       rate1
                                              1
                                             10       0.00015
                                              1
                                             20       0.00004       1.732
                                              1
                                             40       0.00001       1.738
                                              1
                                             80    4.049 × 10−6     1.741
                                              1
                                             160   1.209 × 10−6     1.744
                                              1
                                             320   3.597 × 10−7     1.749



Table 2: Errors and convergence order of the difference scheme (9) for (15) in time variable with α = 1.75,
h = π/3000 and with time delay s = 41 .

                                              τ       E(τ, h)       rate1
                                              1
                                             40       0.00003
                                              1
                                             80       0.00001       1.232
                                              1
                                             160   5.415 × 10−6     1.238
                                              1
                                             320   2.288 × 10−6     1.243
                                              1
                                             640   9.625 × 10−7     1.249



Table 3: Errors and convergence order of the difference scheme (9) for (15) in space variable with α = 1.5,
τ = 1/10000 and with time delay s = 12 .
                                          h       E(τ, h)     rate2
                                          π
                                          4      0.00027
                                          π
                                          8      0.000018     3.872
                                          π                −6
                                         16   1.2525 ×  10    3.880
                                          π
                                         32  8.24587 × 10−8 3.925
                                          π
                                         64  5.28023 × 10−9 3.965
                                          π
                                         128 3.33926 × 10−10 3.983


Acknowledgements
This work was supported by Government of the Russian Federation Resolution N 211 of March 16, 2013.

References
[1] E. A. Abdel-Rehim, A. M. A. El-Sayed and A. S. Hashem. Simulation of the approximate solutions of the
    time-fractional multi-term wave equations. Computers and Mathematics with Applications, 73(6):1134–1154,
    2017.




                                                           331
[2] O. P. Agrawal. A general solution for a fourth-order fractional diffusion wave equation defined in a bounded
    domain. Computers and Structures, 79(16):1497–1501, 2001.

[3] K. Burrage, A. Cardone, R. D’Ambrosio and B. Paternoster. Numerical solution of time fractional diffusion
    systems. Applied Numerical Mathematics, 116:82–94, 2017.

[4] J. Batzel and F. Kappel. Time delay in physiological systems: Analyzing and modeling its impact.
    Mathematical Biosciences, 234(2):61–74, Dec 2011.

[5] R. V. Culshaw, S. Ruan and G. Webb. A mathematical model of cell-to-cell spread of HIV-1 that includes a
    time delay. J. Math. Biol., 46(5):425–444, May 2003.

[6] L. Chang, G. Sun, Z. Wang and Z. Jin. Rich dynamics in a spatial predator-prey model with delay. Applied
    Mathematics and Computation, 256:540–550, 2015.

[7] P. Chen, X. Zhang and Y. Li. Study on fractional non-autonomous evolution equations with delay.
    Mathematics with Applications, 73(5):794–803, 2017.

[8] R. Du, W. R. Cao and Z. Z. Sun. A compact difference scheme for the fractional diffusion-wave equation.
    Applied Mathematical Modelling, 34(10):2998–3007, 2010.

[9] M. Dehghan and R. Salehi. Solution of a nonlinear time-delay model in biology via semi-analytical approaches.
    Computer Physics Communications, 181:1255–1265, 2010.

[10] R. Gorenflo, Y. Luchko and F. Mainardi. Wright functions as scale-invariant solutions of the diffusion-wave
    equation. Journal of Computational and Applied Mathematics, 118(1):175–191, 2000.

[11] Z. Hao, K. Fan, W. Cao and Z. Sun. A finite difference scheme for semilinear space-fractional diffusion
    equations with time delay. Applied Mathematics and Computation, 275:238–254, 2016.

[12] A.S. Hendy, R.H. De Staelen, V.G. Pimenov. A semi-linear delayed diffusion-wave system with distributed
    order in time. Numer Algor, doi:10.1007/s11075-017-0344-7, 2017.

[13] Z. Jackiewicz, H. Liu, B. Li and Y. Kuang. Numerical simulations of traveling wave solutions in a drift
    paradox inspired diffusive delay population model. Mathematics and Computers in Simulation, 96:95–103,
    2014.

[14] Y. Kian and M. Yamamoto. On existence and uniqueness of solutions for semilinear fractional wave
    equations. Fractional Calculus and Applied Analysis, 20(1):117–138, Jan 2017.

[15] P. Liu. Periodic solutions in an epidemic model with diffusion and delay. Applied Mathematics and
    Computation, 265:275–291, 2015.

[16] R. R. Nigmatullin. The realization of the generalized transfer equation in a medium with fractal geometry.
    Physica status solidi (b), 133(1):425–430, Jan 1986.

[17] V. G. Pimenov, A. S. Hendy and R. H. De Staelen. On a class of non-linear delay distributed order fractional
    diffusion equations. Journal of Computational and Applied Mathematics, 318:433–443, 2017.

[18] A. A. Samarskii and V. B. Andreev. Finite difference methods for elliptic equations. Nauka, Moscow, 1976.
    (in Russian) = А.А. Самарский, В.Б. Андреев. Разностные методы для эллиптических уравнений.
    Наука, Москва, 1976.

[19] J. Shi and R. Shivaji. Persistence in reaction diffusion models with weak allee effect. Journal of Mathematical
    Biology, 52(6):807–829, Jun 2006.

[20] Z. Sun and X. Wu. A fully discrete difference scheme for a diffusion-wave system. Applied Numerical
    Mathematics, 56(2):193–209, Feb 2006.

[21] Z. Wang and S. Vong. Compact difference schemes for the modified anomalous fractional sub-diffusion
    equation and the fractional diffusion-wave equation. Journal of Computational Physics, 277:1–15, 2014.




                                                        332
[22] J. Y. Yang, J. F. Huang, D. M. Liang and Y. F. Tang. Numerical solution of fractional diffusion-wave
    equation based on fractional multistep method. Applied Mathematical Modelling, 38(14):3652–3661, 2014.

[23] Q. Zhang, M. Ran and D. Xu. Analysis of the compact difference scheme for the semilinear fractional partial
    differential equation with time delay. Applicable Analysis, 1–18, Jun 2016.
[24] W. Zhang, J. Li and Y. Yang. A fractional diffusion-wave equation with non-local regularization for image
    denoising. Signal Processing, 103:6–15, 2014.
[25] Z. B. Zhang and Z. Z. Sun. A linearized compact difference scheme for a class of nonlinear delay partial
    differential equations. Applied Mathematical Modelling, 37:742–752, 2013.
[26] M. Zheng, F. Liu, Q. Liu, K. Burrage and M. Simpson. Numerical solution of the time fractional reaction-
    diffusion equation with a moving boundary. Journal of Computational Physics, 338:493–510, 2017.




                                                      333