Numerical solution of non-linear fourth order fractional sub-diffusion wave equation with time delay

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 19, 2019;3:44]

Applied Mathematics and Computation xxx (xxxx) xxx

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Numerical solution of non-linear fourth order fractional sub-diffusion wave equation with time delay Sarita Nandal∗, Dwijendra Narain Pandey Department of Mathematics, Indian Institute of Technology Roorkee, Roorkee 247667, India

a r t i c l e

i n f o

Article history: Received 1 November 2018 Revised 18 October 2019 Accepted 28 October 2019 Available online xxx Keywords: Fourth order fractional sub-diffusion equation L2 − 1σ formula Compact difference scheme Stability Convergence

a b s t r a c t In this paper, we constructed a linearized compact difference scheme for fourth order nonlinear fractional sub-diffusion equation with time delay and variable coefficients. The primary purpose of our work is to use the idea of the L2 − 1σ formula for temporal dimension and compact linear operator for spatial dimension. The proposed method is unconditionally stable and convergent to the analytical solution with the order of convergence O(τ 2 + h4 ), where τ and h are temporal and spatial lengths, respectively. Numerical experimentation is carried out to show the efficiency and accuracy of the proposed scheme. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Nowadays, researchers have placed more attention on the development of integral equations and fractional differential equations as these equations can describe numerous chemical, physical phenomenons more precisely than in comparison to integer-order ones and have been used widely in processes like anomalous diffusion process [1], finance, and biological systems [2], hydrology [3] and in many more. Though, various ways are available to solve them theoretically using methods like the Laplace method, Fourier method, Green function method [4,5]. Sometimes, it’s challenging to obtain an analytical solution of such equations theoretically, therefore, development of numerical methods for finding a solution of integral and fractional differential equations seems to be vital and essential. In recent years, a class of fractional differential equations has gained more attention in different domains is timefractional sub-diffusion wave equations (FSDWE). This type of equation is like integrodifferential equations which might be obtained by replacing integer order time derivative with the fractional order time derivative of order α , where α ∈ (0, 1) [6–9] in classical sub-diffusion equations. Time fractional sub-diffusion wave equations have critical applications in signal processing, anomalous diffusion, turbulence, wave propagation [10–12]. Many authors have studied FSDWE theoretically and computationally [13,14]. As we know that finding an analytical solution is not always convenient for practical applicationoriented problems, most of such type of problems can not be solved analytically. Time delay occurs in many real-life applications such as population ecology, cell biology, control theory which are modeled mathematically [15–17]. Many authors had implemented L1 formula (a discrete approximation of Caputo fractional derivative) which have the truncation error of order O(τ 2−α ), [10–14]. Authors in [29] studied neutral delay differential ∗

Corresponding author. E-mail addresses: [email protected], [email protected] (S. Nandal), [email protected] (D. Narain Pandey).

https://doi.org/10.1016/j.amc.2019.124900 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC 2

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

equation of fourth-order and constructed a scheme to enhance the order of convergence for spatial dimension. Cui [19], Liu et al. [18] and Zhuang et al. [32] constructed an implicit difference scheme for an anomalous sub-diffusion wave equation with a source term and discussed the analysis of scheme using discrete energy method. Authors et al. in [13] constructed an implicit difference scheme for first and second-order temporal accuracy and fourth-order improved spatial accuracy for the fractional sub-diffusion wave equation and proved the stability and convergence analysis using Fourier analysis. Gao and Sun [27] constructed a compact difference scheme for fractional sub-diffusion wave equation using L1 approximation for the time-fractional term and proved the unique solvability of the proposed scheme and its stability and convergence using energy method. Authors in [38] presented a novel variable-order fractional nonlinear Klein Gordon model and in [39] proposed a novel numerical technique for approximating the solutions of the two-dimensional linear and nonlinear fractional cable equation. They defined fractional derivative in terms of the Atangana–Baleanu–Caputo sense and implemented a nonstandard implicit compact finite difference method to solve numerically the proposed equation. The structure of the paper is as follows: In Section 3 definitions, lemmas and notations to be used in the derivation of the compact difference scheme are given. In next Section 4, derivation of proposed numerical scheme is done. In the next section, the proposed compact difference scheme is analyzed in terms of stability and convergence. The proposed numerical scheme is unconditionally stable and uniquely solvable. In the next section, numerical experimentation is done through examples by comparing analytical and numerical results. Results are presented in terms of graphs and tables, clearly indicating the order of scheme in spatial and temporal terms having a validation with theoretical results. An excellent agreement is found between numerical and analytical results. The fourth order time fractional sub-diffusion equation is given as follows

∂ α u(x, t ) ∂ 4 u(x, t ) =K + f (x, t ), 0 < x < L, 0 ≤ t < T . α ∂t ∂ x4

(1)

∂ α u(x, t ) ∂ 4 u(x, t ) = a (x ) + b(x ) u(x, t ) + f (x, t, u(x, t ), u(x, t − s )), 0 < x < L, 0 < t < T , α ∂t ∂ x4

(2)

Here, α ∈ (0, 1) is the fractional-order, for α = 1 classical fourth-order diffusion-wave equation is obtained. Various numerical solution approaches such as finite difference method [23], finite element method [26], collocation methods [24,25] etc. has been constructed in fast few years for equations such as (1). To overcome the complexity of the analytical solutions; the need for numerical solutions is increased numerously nowadays. Fractional sub-diffusion delay differential equations are of great importance these days for their real life-based applications such as in the field of medicine [20], economics [21], Physics [22] etc. Present work is the extension to the work of Zhang and Pu [31]. In the current work, we seek to generalize Eq. (1) by including a time-delayed non-linear source term and variable coefficients. Here, we derive a linearized efficient numerical scheme for the class of fourth-order non-linear time-delay sub-diffusion fractional differential equation with variable coefficients. Consider the fourth-order non-linear fractional sub-diffusion equation with time delay and variable coefficients

Initial and Dirichlet conditions

u(x, t ) = ψ (x, t ), u(0, t ) = g1 (t ),

0 ≤ x ≤ L,

t ∈ [−s, 0],

u(L, t ) = g2 (t ),

∂ 2 u ( 0, t ) = h1 (t ), ∂ x2

0 ≤ t ≤ T,

∂ 2 u(L, t ) = h2 (t ), 0 ≤ t ≤ T . ∂ x2

(3) (4)

(5)

Here s > 0 is delay, a(x) and b(x) are variable coefficients and f (x, t, u(x, t ), u(x, t − s )) stands for non-linear time delayed source term. a(x), b(x), ψ (x, t), g1 (t), g2 (t), h1 (t), h2 (t) are all given and sufficiently smooth functions. ψ (0, t ) = g1 (t ), ψ (L, t ) = g2 (t ), ψ  (0, t ) = h (t ), and ψ  (L, t ) = h2 (t ). α 1 Fractional derivative ∂ ∂ut(αx,t ) is considered in Caputo sense as follows C α 0 Dt u

(x, t ) ≡

∂ α u(x, t ) = ∂tα



0

t

(t − ξ )−α

∂ u(x, ξ ) d ξ , 0 < α < 1. ∂ξ

(6)

In (2), we have two type of complexities one is delay term and second is non-linear source term. Here, we seek to obtain a linearized numerical scheme for (2)–(5) which is uniquely solvable, stable and convergent. Throughout the paper, we consider that the obtained solution u(x, t) along with source function f(x, t, μ, ν ) are sufficiently smooth likewise considered in [37] in the following sense:  • Let m be the integer satisfying ms ≤ T ≤ (m + 1 )s. Define Ir = (rs, (r + 1 )s ), r = -1, 0, ..., m-1, Im = (ms, T ), I = m q=−1 Iq 8,3 and assume that u(x, t ) ∈ Cx,t ([0, L] × [−s, T ] ), • the partial derivatives fμ (x, t, μ, ν ) and fν (x, t, μ, ν ) are continuous in the 0 neighborhood of the solution. Define c1 = sup | fμ (x, t, μ, ν ) + 1 , u(x, t − s ) + 2 |, c2 = sup | fν (x, t, μ, ν ) + 1 , u(x, t − s ) + 2 |.

Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

3

2. Notations For approximation, we consider partition of the mesh ht = h × t . Introducing two positive integers M and N, let h = τ = ns (n > 0 is a positive constant) are the space and temporal step lengths, respectively. Define x j = j h, 0 ≤ j ≤ M, tk = k τ , −n ≤ k ≤ N, h = {x j |0 ≤ j ≤ M}, t = {tk | − n ≤ k ≤ N}. In addition, σ = 1 − α2 , denote tk−1+σ = (k − 1 + σ )τ . Denote Vh = {u|u = (u0 , u1 , u2 , . . . , uM ), u0 = uM = 0} as the grid function space on h . For any grid function u ∈ Vh , we introduce the notations below L M,

1 (u j+1 − u j ), 0 ≤ j ≤ M − 1, h   1 δx u j+ 12 − δx u j− 12 , 1 ≤ j ≤ M − 1. δx2 u j = h

δx u j+ 12 =

Introducing the discrete inner products and and the corresponding norms for any u, v ∈ Vh as follows:

(u, v ) = h

||u|| =

M−1 



(u, v ) = h

(δx u, δx v ) = h

u jv j,

j=1

(u, u ),

M−1 

||δx u|| =

( u j )v j ,



M−1 

1 2

δx u j + v j+ 12 ,

j=0

(δx u, δx u ),  ||u|| = (u, u ).

j=1

3. Preliminaries Definition 1. Alikhanov [28] constructed a new discrete formula called L2 − 1σ formula as discrete approximation of Caputofractional derivative of second order temporal convergence for α ∈ (0, 1). Defining C α 0 Dt

ukj −1+σ

  k−1 

j (k ) 0 τ −α (k ) k (k ) (k ) = c u − ck− j−1 − ck− j u − ck−1 u , (2 − α ) 0 j=1

(7)

Here, ( · ) is the Gamma function. Where, c0(k ) = a0 , when k = 0, and for k ≥ 1,

c(jk ) =

⎧ ⎨a0 + b1 ,

j = 1,

ak + bk+1 − bk ,

1 < j ≤ k − 2,

ak − bk ,

j = k − 1.

⎩

(8)

where

a l = ( l + σ ) 1 −α − ( l − 1 + σ ) 1 −α , l ≤ 1,  1   1 2 −α 2 −α bl = (l + σ ) − (l − 1 + σ ) − ( l + σ ) 1 −α + ( l − 1 + σ ) 1 −α , 2−α 2

a0 =

σ 1 −α ,

l ≤ 1,

l ≤ 1,

We now introduce the following Lemmas which will be used in the analysis of the difference scheme. Lemma 1. Alikhanov [28] gave the error estimation of the L2 − 1σ formula to approximate the Caputo fractional derivative. Suppose u(t) ∈ C3 [0, tk ], it holds that

  −α C α  0 Dt u(t )|t=tk−1+σ − C0 Dtα uk−1+σ  ≤ (4σ − 1 )σ |u(3) (t )|τ (3−α ) .   12 (2 − α ) tmax 0 ≤t≤tk

Below given Lemmas 2 and 3 describes a few properties of L2 − 1σ formula. Lemma 2 [28]. Consider that α ∈ (0, 1) and σ = 1 − α2 . c (jk ) (0 ≤ j ≤ k − 1 ) is defined by (8), then following inequality holds

c(jk ) >

1−α ( k + σ ) −α , 2

) ) ck(k−1 < ck(k−2 < . . . < c2(k ) < c1(k ) < c0(k ) .

(9) (10)

Lemma 3 [28]. Consider u(x, t ) = {ukj |0 ≤ j ≤ M, 0 ≤ k ≤ N} is a grid function defined on τ , then below given inequality holds

1C α D ( uk )2 ≤ 2 0 tk−1+σ j

  k k−1 C α σ u j + ( 1 − σ )u j 0 Dtk−1+σ u j .

(11)

Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC 4

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

Below given definition and Lemma are to be implemented for spatial dimension. Definition 2 describes the compact linear operator. Definition 2 [30]. The compact linear operator is given as follows

( u ) j =

⎧1 ⎨ (u j+1 + 10u j + u j−1 ), ⎩

12

u j,

1 ≤ j ≤ M − 1, j = 0,

Lemma 4 [28]. Denote ξ (t ) =

g (xi ) = δx2 g(xi ) +

(1 − t )3 [5 − 3(1 − t )2 ]

h4 360

 1

or

M.

and consider function g(x) ∈ C6 [0, L]. It holds that

g(6) (xi − th ) + g(6) (xi + th )

0

 ξ (t )d t , 1 ≤ i ≤ M − 1.

4. Derivation of the numerical scheme ( 8,3 ) Define the grid functions U kj = u(x j , tk ), V jk = v(x j , tk ), 0 ≤ j ≤ M, −n ≤ k ≤ N. Suppose u(x, t ) ∈ Cx,t ([0, L] × [−s, T ] ). Now, we consider the derivation of finite difference scheme. For the sake of simplifying the considered Eq. (2), we use the 2 reduced order method and suppose v(x, t) = ∂ ∂ux(2x,t ) , then the problems (2)–(5) is equivalent to

∂ α u(x, t ) ∂ 2 v(x, t ) = a (x ) + b(x ) u(x, t ) + f (x, t, u(x, t ), u(x, t − s )), 0 < x < L, 0 < t < T , α ∂t ∂ x2 v(x, t ) =

∂ 2 u(x, t ) , 0 < x < L, 0 < t < T , ∂ x2

u(x, t ) = ψ (x, t ), u(0, t ) = g1 (t ),

0 ≤ x ≤ L,

t ∈ [−s, 0],

u(L, t ) = g2 (t ),

(12) (13) (14)

0 ≤ t ≤ T,

(15)

v(0, t ) = h1 (t ), v(L, t ) = h2 (t ), 0 ≤ t ≤ T .

(16)

Now, we consider the (12) and (13) at the mesh point (x j , tk−1+σ ), then, we have

∂ α u(x j , tk−1+σ ) ∂ 2 v(x j , tk−1+σ ) = a (x j ) + b(x j ) u(x j , tk−1+σ ) α ∂t ∂ x2 + f (x j , tk−1+σ , u(x j , tk−1+σ ), u(x j , tk−1+σ −n )), 1 ≤ j ≤ M − 1, 1 ≤ k ≤ N − 1, v(x j , tk−1+σ ) = U jk = ψ (x j , tk ),

∂ 2 u(x j , tk−1+σ ) , 1 ≤ j ≤ M − 1, 1 ≤ k ≤ N − 1, ∂ x2 0 ≤ j ≤ M,

−n ≤ k ≤ N,

(17)

(18) (19)

U0k = g1 (tk ),

k UM = g2 (tk ),

1 ≤ k ≤ N,

(20)

V0k = h1 (tk ),

VMk = h2 (tk ),

1 ≤ k ≤ N.

(21)

Implementing the compact difference operator  (2) on both sides of Eqs. (17) and (18) gives



 ∂ α u(x j , tk−1+σ ) ∂ 2 v(x j , tk−1+σ )  =  a ( x ) + (b(x j ) u(x j , tk−1+σ )) j ∂tα ∂ x2 +  f (x j , tk−1+σ , u(x j , tk−1+σ ), u(x j , tk−1+σ −n )), 1 ≤ j ≤ M − 1, 1 ≤ k ≤ N − 1,

v(x j , tk−1+σ ) = 

∂ 2 u(x j , tk−1+σ ) , 1 ≤ j ≤ M − 1, 1 ≤ k ≤ N − 1. ∂ x2

(22)

(23)

Using Taylor’s series following expressions can be easily obtained

u(x j , tk−1+σ ) = σ u(x j , tk ) + (1 − σ )u(x j , tk−1 ).

(24)

Denoting

U jk−1+σ = σ U jk + (1 − σ )U jk−1 ,

(25)

Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

5

and

U jk−1+σ −n = σ U jk−n + (1 − σ )U jk−n−1 .

(26)

Linearizing the non-linear time-delay source term f(x, t, μ, ν ) using Taylor’s series gives

f jk−1+σ = f (x j , tk−1+σ , 2σ U jk+1 + (2 − 3σ )U jk − (1 − σ )U jk−1 , σ U jk−n + (1 − σ )U jk−n−1 ),

(27)

denoting

Uˆ jk−1+σ = 2σ U jk+1 + (2 − 3σ )U jk − (1 − σ )U jk−1 .

(28)

Using Lemma 4, we have



∂ 2 u(x j , tk−1+σ ) ∂ 2 u ( x j , tk ) ∂ 2 v(x j , tk−1 ) = σ + ( 1 − σ ) + O(τ 2 ) = δx2U jk−1+σ + O(τ 2 + h4 ). ∂ x2 ∂ x2 ∂ x2

(29)

Similarly, we have

v(x j , tk−1+σ ) = V jk−1+σ + O(τ 2 ), 

(30)

∂ 2 v(x j , tk−1+σ ) = δx2V jk−1+σ + O(τ 2 + h4 ). ∂ x2

(31)

Using Definition 1 and substituting Eqs. (29)–(31) into (22) and (23); we obtain

Dtαk−1+σ U j = (a(x j )δx2V jk−1+σ ) + (b(x j ) U jk−1+σ ) +  f (x j , tk−1+σ , 2σ U jk+1 + (2 − 3σ )U jk − (1 − σ )U jk−1 , σ U jk−n + (1 − σ )U jk−n−1 ) + (R1 )kj , V jk−1+σ = δx2U jk−1+σ + (R2 )kj , where

1 ≤ j ≤ M − 1,

1 ≤ j ≤ M − 1,

1 ≤ k ≤ N − 1,

1 ≤ k ≤ N − 1,

|(R1 )kj | + |(R2 )kj | ≤ c1 (τ 2 + h4 ), 1 ≤ j ≤ M − 1, 1 ≤ k ≤ N − 1.

(32) (33) (34)

where c1 is a positive constant independent of τ and h. Omitting the small terms (R1 )kj and (R2 )kj in (32) and (33), we obtain the compact difference scheme for the problems (2)–(5). By replacing U kj by ukj and V jk by vkj we get the proposed scheme as follows

  k−1  τ −α (k ) (k ) (k ) k j k k c u − (ck− j−1 − ck− j )u − ck−1 ψ j (2 − α ) 0 j=1

= a(x j )δx2 vkj −1+σ + b(x j )ukj −1+σ +  f (x j , tk−1+σ , 2σ ukj +1 + (2 − 3σ )ukj − (1 − σ )ukj −1 , σ ukj −n + (1 − σ )ukj −n−1 ), vkj −1+σ = δx2 ukj −1+σ , ukj = ψ (x j , tk ), uk0 = g1 (tk ),

1 ≤ j ≤ M − 1, 1 ≤ k ≤ N − 1, 1 ≤ j ≤ M − 1,

0 ≤ j ≤ M,

ukM = g2 (tk ),

1 ≤ k ≤ N − 1,

−n ≤ k ≤ N,

(35) (36) (37)

1 ≤ k ≤ N,

(38)

vk0 = h1 (tk ), vkM = h2 (tk ), 1 ≤ k ≤ N.

(39)

5. Solvability, stability and convergence In this section, solvability, stability and convergence will be proved for the Eqs. (2)–(5). Following Lemmas and Theorems will be used to achieve the target of the analysis of the scheme. Lemma 5 [34]. For any grid function u ∈ Vh , it holds that

||u||2 ≤

L2 ||δx u||2 , 6

||δx u||2 ≤

L2 2 2 ||δx u|| . 6

(40)

(41)

Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC 6

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

Algorithm 1 Algorithm of the numerical scheme in (35)–(39). L Input: Write L, M, s, n and N. Compute spatial and temporal step sizes respectively; h = M , and τ = ns . Output: Solution of Eqs. (2)–(5) i.e. function u(x, t ) 1. Set the initial and boundary conditions, i.e. ψ (x, t ), g1 (t ), g2 (t ), h1 (t ), and h2 (t ) for j = 1, M, and k = 1, N. 2. for j=2:M-1 and k=2:N-1, compute each term of Eq. (35) separately 3. In next step, keep every term associated with function u(x, t ) on Left Hand Side and rest terms associated with x and t only on Right Hand Side 4. We will obtain a matrix for Left Hand Side and a column vector for Right Hand Side End for 5. Now compute A−1 b, where matrix A is concerned for Left Hand Side and column vector b represents Right Hand Side.

Proof. For the proof of Eq. (40), reader can follow [34] and Eq. (41) can be easily derived with the help of (40), discrete Green formula and the Cauchy–Schwartz inequality.  Lemma 6 [31]. For any grid function u ∈ Vh , it holds that

1 ||u||2 ≤ ||u||2 ≤ ||u||2 . 3

(42)

Proof. For the proof of Eq. (42), reader can follow [31]. Inequality (42) can be obtained using discrete Green formula and inverse estimates given below

4 ||δx u||2 , h2 4 ≤ 2 ||u||2 . h

||δx2 u||2 ≤ ||δx u||2

(43) 

ukj

Lemma 7. Consider that (35)–(39)

and

vkj

with 0 ≤ j ≤ M, −n ≤ k ≤ N are solutions of the constructed difference scheme given in

C0 Dtαk−1+σ u j = (a(x j )δx2 vkj −1+σ ) + (b(x j )ukj −1+σ ) +  f (x j , tk−1+σ , 2σ ukj +1 + (2 − 3σ )ukj − (1 − σ )ukj −1 , σ ukj −n + (1 − σ )ukj −n−1 ), vkj −1+σ = δx2 ukj −1+σ + H1k−1+σ , ukj = ψ jk ,

0 ≤ j ≤ M,

uk0 = g1 (tk ),

1 ≤ j ≤ M − 1,

1 ≤ k ≤ N − 1,

1 ≤ j ≤ M − 1,

(45)

−n ≤ k ≤ 0,

ukM = g2 (tk ),

(46)

1 ≤ k ≤ N,

(47)

vk0 = h1 (tk ), vkM = h2 (tk ) 1 ≤ k ≤ N. Then, following inequality holds

||u || ≤ ||u || k

2

0

2



+ 2 T α (1 − α )

+ 8(1 − σ )2 max

−n≤k≤N

(44)

(48)

L4 max ||gk−1+σ ||2 + 4 max ||H1k−1+σ ||2 18 −n≤k≤N −n≤k≤N

 ||ukj ||2 + max ||ukj −n+σ −1 ||2 + K1 max ||ukj −1+σ ||2 . −n≤k≤N

−n≤k≤N

(49)

Where H1k−1+σ is the initial perturbation and stability analysis will be done concerning the perturbation. Proof. Taking inner product of (44), (45) by ukj −1+σ and vkj −1+σ respectively and adding up for j from 1 to M − 1, then we obtain

(C0 Dtαk−1+σ u j , ukj −1+σ ) = (a(x j )δx2 vkj −1+σ , ukj −1+σ ) + (b(x j )u(x j , tk−1+σ ), ukj −1+σ ) + ( f (x j , tk−1+σ , 2σ ukj +1 + (2 − 3σ )ukj − (1 − σ )ukj −1 , σ ukj −n + (1 − σ )ukj −n−1 ), ukj −1+σ ), (50)

(vkj −1+σ , vkj −1+σ ) = (δx2 ukj −1+σ , vkj −1+σ ) + (H1k−1+σ , vkj −1+σ ).

(51)

Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

7

Below given Eq. (52) can be obtained easily using the fact that uk , vk ∈ Vh

(δx2 vkj −1+σ , ukj −1+σ ) =

  h2 δx2 vkj −1+σ , ukj −1+σ + δx2 ukj −1+σ 12

= − (δx vkj −1+σ , δx ukj −1+σ ) +

h2 2 k−1+σ 2 k−1+σ (δ v , δx u j ) 12 x j

= (δx2 ukj −1+σ , vkj −1+σ )

(δx2 vkj −1+σ , ukj −1+σ ) = (δx2 ukj −1+σ , vkj −1+σ ).

(52)

So adding, (50) and (51), we get

(C0 Dtαk−1+σ u, uk−1+σ ) + ||vk−1+σ ||2 = (a(x j )δx2 vkj −1+σ , ukj −1+σ ) + (b(x j )u(x j , tk−1+σ ), ukj −1+σ ) + ( f (x j , tk−1+σ , 2σ ukj +1 + (2 − 3σ )ukj − (1 − σ )ukj −1 , σ ukj −n + (1 − σ )ukj −n−1 ), ukj −1+σ ) + (δx2 ukj −1+σ , vkj −1+σ ) + (H1k−1+σ , vk−1+σ ).

(53)

Using Eq. (45), we can write the following inequality

1 1 1 ||vk−1+σ ||2 = ||δx2 ukj −1+σ ||2 + ||H1k−1+σ ||2 + (δx2 ukj −1+σ , H1k−1+σ ). 2 2 2

(54)

Substituting Eq. (54) into (53), we get

(C0 Dtαk−1+σ u, uk−1+σ ) + ||δx2 ukj −1+σ ||2 + ||H1k−1+σ ||2 + 2(δx2 ukj −1+σ , H1k−1+σ ) = (a(x j )δx2 vkj −1+σ , ukj −1+σ ) + (b(x j )u(x j , tk−1+σ ), ukj −1+σ ) + ( f (x j , tk−1+σ , 2σ ukj +1 + (2 − 3σ )ukj − (1 − σ )ukj −1 , σ ukj −n + (1 − σ )ukj −n−1 ), ukj −1+σ ) + (δx2 ukj −1+σ , vkj −1+σ ) + (H1k−1+σ , vk−1+σ ).

(55)

Using Lemma 3, Eq. (55) becomes

1C α D ||u||2 + ||δx2 ukj −1+σ ||2 + ||H1k−1+σ ||2 + 2(δx2 ukj −1+σ , H1k−1+σ ) 2 0 tk−1+σ ≤ (a(x j )δx2 vkj −1+σ , ukj −1+σ ) + (b(x j )u(x j , tk−1+σ ), ukj −1+σ ) + ( f (x j , tk−1+σ , 2σ ukj +1 + (2 − 3σ )ukj − (1 − σ )ukj −1 , σ ukj −n + (1 − σ )ukj −n−1 ), ukj −1+σ ) + (δx2 ukj −1+σ , vkj −1+σ ) + (H1k−1+σ , vk−1+σ ).

(56)

Using -Cauchy inequality P Q ≤ 2 P 2 + 21 Q 2 , where > 0 (more precisely we have considered = 1), and |a(xj )| ≤ M1 and |b(xj )| ≤ M2 in h , then Eq. (56) becomes

1 1C α 1 D ||u||2 + 2||δx2 ukj −1+σ ||2 + 2||H1k−1+σ ||2 ≤ ||a(x j )δx2 vkj −1+σ ||2 + ||ukj −1+σ ||2 2 0 tk−1+σ 2 2 1 1 1 2 k−1+σ 2 1 2 k−1+σ 2 k−1+σ 2 + ||b(x j )u(x j , tk−1+σ )|| + ||u j || + ||δx u j || + ||v j || 2 2 2 2 + ( f (x j , tk−1+σ , 2σ ukj +1 + (2 − 3σ )ukj − (1 − σ )ukj −1 , σ ukj −n + (1 − σ )ukj −n−1 ), ukj −1+σ ) +

1 1 ||H1k−1+σ ||2 + ||vk−1+σ ||2 . 2 2

Above equation can be re-written as Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC 8

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

M 1C α 1 D ||u||2 + 2||δx2 ukj −1+σ ||2 + 2||H1k−1+σ ||2 ≤ 1 ||δx2 vkj −1+σ ||2 + ||ukj −1+σ ||2 2 0 tk−1+σ 2 2 M2 1 1 2 k−1+σ 2 1 2 k−1+σ 2 k−1+σ 2 + ||u(x j , tk−1+σ )|| + ||u j || + ||δx u j || + ||v j || 2 2 2 2 + ( f (x j , tk−1+σ , 2σ ukj +1 + (2 − 3σ )ukj − (1 − σ )ukj −1 , σ ukj −n + (1 − σ )ukj −n−1 ), ukj −1+σ ) +

1 1 ||H1k−1+σ ||2 + ||vk−1+σ ||2 . 2 2

(57)

Using Eqs. (24)–(28) and -Cauchy inequality, where > 0, we can easily obtain the following

( f (x j , tk−1+σ , 2σ ukj +1 + (2 − 3σ )ukj − (1 − σ )ukj −1 , σ ukj −n + (1 − σ )ukj −n−1 ), ukj −1+σ ) ≤ ||gkj −1+σ ||2 + 8(1 − σ )2 ||ukj ||2 + 5||ukj −1+σ ||2 + ||ukj −n+σ −1 ||2 .

(58)

In Eq. (58), g(x, t) is function of time and space dimension only a part of non-linear time delay source term f (x, t, u(x, t ), u(x, t − s )). Rewriting Eq. (57) using (52) and (58) we get C α 0 Dtk−1+σ

||u||2 ≤ K1 ||ukj −1+σ ||2 + 4||H1k−1+σ ||2 + ||gkj −1+σ ||2 + 8(1 − σ )2 ||ukj ||2 + ||ukj −n+σ −1 ||2 .

(59)

here K1 is a constant independent of τ and h. Using Lemmas 2 and 3, we have

||uk−1+σ ||2 ≤

L4 ||δ 2 uk−1+σ ||2 . 36 x j

(gk−1+σ , uk−1+σ ) ≤

(60)

18 1 L4 L4 ||uk−1+σ ||2 + ||gk−1+σ ||2 ≤ ||uk−1+σ ||2 + ||gk−1+σ ||2 . 4 18 2 18 L

(61)

Using Eqs. (60) and (61), we get C α 0 Dtk−1+σ

||u||2 ≤ K1 ||ukj −1+σ ||2 +

L4 k−1+σ 2 ||g || + 4||H1k−1+σ ||2 + 8(1 − σ )2 ||ukj ||2 + ||ukj −n+σ −1 ||2 . 18

(62)

Eq. (62) can be rewritten as follows

c0k ||uk ||2 ≤

k−1 

 L4 (ckk− j−1 − ckk− j )||u j ||2 + ckk−1 ||uk ||2 + μ K1 ||ukj −1+σ ||2 + ||gk−1+σ ||2 18

j=1



+ 4||H1k−1+σ ||2 + 8(1 − σ )2 ||ukj ||2 + ||ukj −n+σ −1 ||2 ,

(63)

where,

μ = τ α (2 − α ) = T α (1 − α ) (1 − α ) N−α   α −α < T α (1 − α ) (1 − α ) k − 2

2ckk−1 T α

<

(1 − α ).

(64)

From Lemma 2, we know that



ckk−1

1−α α > k−1− 2 2





−α

1−α α > k− 2 2

 −α

,

1 ≤ k ≤ N.

(65)

Substituting (64) into (63), we get

c0k ||uk ||2 ≤

k−1 

  L4 (ckk− j−1 − ckk− j )||u j ||2 + ckk−1 ||u0 ||2 + 2T α (1 − α ) K1 ||ukj −1+σ ||2 + ||gk−1+σ ||2 18

j=1



+ 4||H1k−1+σ ||2 + 8(1 − σ )2 ||ukj ||2 + ||ukj −n+σ −1 ||2 Consider



D = ||u0 ||2 + 2T α (1 − α ) K1 ||ukj −1+σ ||2 +

.

(66)



L4 k−1+σ 2 ||g || + 4||H1k−1+σ ||2 + 8(1 − σ )2 ||ukj ||2 + ||ukj −n+σ −1 ||2 . 18

Below given inequality can be proved easily by induction: Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

9

||uk ||2 ≤ D, 1 ≤ k ≤ N, c0k ||uk ||2 ≤ c0k D. It holds obviously when k = 0, and assume that the conclusion is valid for k = 1, 2, . . . , M − 1, i.e.,

||uk ||2 ≤ D, 1 ≤ k ≤ M − 1. Then, for 2 ≤ k ≤ N from (57), we have

c0k ||uk ||2 ≤

k−1 

(ckk− j−1 − ckk− j )||uk ||2 + ckk−1 D

j=1

≤

k−1 

(ckk− j−1 − ckk− j )D + ckk−1 D = c0k D.

j=1



Which completes the proof. 5.1. Stability

Theorem 1 (Stability). With the help of Lemma 7, following theorem can be written as follows: The proposed compact difference scheme (35)–(39) is unconditionally stable for given function ψ (x, t) and the source function f(x, t, μ, ν ). 5.2. Solvability Theorem 2 (Solvability). The established compact difference scheme (35)–(39) is uniquely solvable. From [33], it needs only to prove uniqueness of linear homogeneous system as given below has only trivial solution.

C0 Dtαk−1+σ u j − (a j δx2 vkj −1+σ ) − (b j ukj −1+σ ) = 0,

1 ≤ j ≤ M − 1, 1 ≤ k ≤ N − 1,

(67)

vkj −1+σ = δx2 ukj −1+σ ,

1 ≤ j ≤ M − 1, 1 ≤ k ≤ N − 1,

(68)

ukj = 0,

0 ≤ j ≤ M,

−n ≤ k ≤ 0,

(69)

uk0 = 0,

ukM = 0,

1 ≤ k ≤ N,

(70)

vk0 = 0, vkM = 0, 1 ≤ k ≤ N.

(71)

Proof. Using Theorem 1 from [31] and Lemma 7; Eqs. (35)–(39) possesses a unique solution.



5.3. Convergence Theorem 3 (Convergence) ([31], Lemma 6). Suppose {U kj |0 ≤ j ≤ M, −n ≤ k ≤ N} and {ukj |0 ≤ j ≤ M, −n ≤ k ≤ N} are solutions of the (2)–(5) and the constructed difference scheme given in (35)–(39), then it holds that

||ek || ≤ c(τ 3−α + h4 ). where

ekj

=

U kj

− ukj

and

ekj

=

(72) V jk

− vkj

Proof. Subtracting Eqs. (3)–(8) from (35)–(39), we get the following error equations as given below

etαk−1+σ u j = (a(x j )δx2 ekj −1+σ ) + (b(x j )ekj −1+σ ) +  f (x j , tk−1+σ , 2σ ekj +1 + (2 − 3σ )ekj − (1 − σ )ekj −1 , σ ekj −n + (1 − σ )ekj −n−1 ), ekj −1+σ = δx2 ekj −1+σ ,

1 ≤ j ≤ M − 1,

1 ≤ j ≤ M − 1,

1≤k≤N−1

(73)

1 ≤ k ≤ N − 1.

(74)

Using Lemma 7, we get the following inequality



||e || k

2

≤ 2 T α (1 − α ) 4 max

−n≤k≤N

||

H1k−1+σ

 || + 8(1 − σ ) max ||e || + max ||e 2

2

k

−n≤k≤N

2

−n≤k≤N

k−n+σ −1

|| + K1 max ||e 2

−n≤k≤N

k−1+σ

||

2

.

(75) Using Lemma 6 and Eq. (75), the claimed inequality (72) is obtained.



Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC 10

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx Table 1 The computational error and convergence order in temporal dimension for Example 1 by varying τ .

α

τ

Using Eq. (35) EL2 (h, τ )

Using [36] EL2 (h, τ )

Using [36] Order(τ )

Using Eq. (35) Order(τ )

Using Eq. (35) EL∞ (h, τ )

Using [36] EL∞ (h, τ )

Using [36] Order(τ )

Using Eq. (35) Order(τ )

0.2

1/20 1/50 1/80 1/20 1/50 1/80 1/20 1/50 1/80 1/20 1/50 1/80

7.0315e-004 1.3840e-005 1.0069e-006 7.3138e-004 5.1235e-005 6.7321e-006 4.6659e-004 2.4806e-005 5.5179e-006 1.0027e-004 2.5144e-005 6.8951e-006

4.1255e-003 3.5693e-003 8.4287e-004 2.2360e-003 6.4638e-003 4.7416e-004 1.8647e-003 1.2005e-004 6.9611e-004 9.4621e-003 5.3685e-004 7.4692e-004

1.7698 1.7881 1.7914 1.5784 1.5883 1.5911 1.3418 1.3647 1.3878 1.1576 1.1769 1.1839

1.9768 1.9828 1.9964 1.9836 1.9907 1.9987 1.9844 1.9953 1.9996 1.9867 1.9923 2.0023

5.1258e-004 4.3832e-005 4.1435e-006 7.5833e-004 7.1021e-005 6.3554e-006 2.9711e-004 8.4763e-005 3.7392e-006 6.0016e-004 2.2358e-005 1.9915e-006

5.2141e-003 1.3625e-003 9.5246e-004 4.3620e-003 6.4921e-004 2.3725e-004 7.2468e-003 8.3693e-004 6.8044e-005 2.9631e-003 3.6002e-004 9.7689e-004

1.7714 1.7886 1.7919 1.5818 1.5892 1.5963 1.3510 1.3713 1.3883 1.1607 1.1778 1.1845

1.9844 1.9901 1.9965 1.9817 1.9932 1.9984 1.9930 1.9967 2.0025 1.9971 1.9994 2.0019

0.4

0.6

0.8

Table 2 The computational error and convergence order in spatial dimension for Example 1 by varying h.

α

h

Using [36] EL2 (h, τ )

Using Eq. (35) EL2 (h, τ )

Using Eq. (35) Order(h)

Using [36] EL∞ (h, τ )

Using Eq. (35) EL∞ (h, τ )

Using Eq. (35) Order(h)

0.2

1/10 1/20 1/40 1/10 1/20 1/40 1/10 1/20 1/40 1/10 1/20 1/40

2.3471e-006 1.6328e-007 6.2845e-007 3.1874e-006 7.2399e-007 1.0468e-007 8.8911e-006 9.3967e-006 1.7468e-007 5.5676e-006 2.4601e-007 6.3942e-007

5.5112e-006 6.6244e-007 3.1088e-008 8.6158e-006 6.7366e-007 4.9421e-008 9.6871e-006 6.3962e-007 7.2764e-008 3.6833e-006 5.9105e-007 1.4687e-008

3.9812 3.9873 3.9937 3.9830 3.9888 3.9947 3.9832 3.9901 3.9964 3.9851 3.9944 3.9985

9.1255e-006 6.6341e-007 7.5208e-007 1.4392e-006 2.7106e-007 6.1857e-007 8.3069e-006 3.4682e-007 7.3988e-007 9.7651e-006 2.9386e-007 8.4620e-007

1.7411e-006 7.6101e-007 4.2564e-008 3.4816e-006 1.9481e-007 8.3952e-008 6.1877e-006 9.6255e-007 2.3607e-008 7.9901e-006 5.2401e-007 9.1966e-008

4.0041 3.9917 3.9905 4.0136 3.9914 3.9901 4.0033 3.9911 3.9878 4.0007 3.9910 3.9825

0.4

0.6

0.8

Table 3 The computational error and convergence order in temporal dimension for Example 2 by varying τ .

α

τ

Using Eq. (35) EL2 (h, τ )

Using [36] EL2 (h, τ )

Using [36] Order(τ )

Using Eq. (35) Order(τ )

Using Eq. (35) EL∞ (h, τ )

Using [36] EL∞ (h, τ )

Using [36] Order(τ )

Using Eq. (35) Order(τ )

0.2

1/20 1/50 1/80 1/20 1/50 1/80 1/20 1/50 1/80 1/20 1/50 1/80

1.1137e-004 6.3101e-005 2.1306e-006 4.4288e-004 5.7024e-005 3.4367e-006 3.1043e-004 7.5013e-005 5.5179e-006 6.4331e-004 8.2390e-005 1.5813e-006

8.3988e-003 1.6973e-004 2.1451e-004 4.2582e-003 1.4629e-004 6.6892e-004 7.7235e-003 9.1159e-003 4.7546e-004 6.8570e-003 2.9202e-004 9.6683e-004

1.7316 1.7458 1.7784 1.5623 1.5782 1.5818 1.3781 1.3803 1.3910 1.1874 1.1880 1.1892

1.9907 1.9934 1.9961 1.9938 1.9979 2.0005 1.9988 1.9999 2.0012 1.9989 1.9991 2.0016

7.2154e-004 2.6431e-005 3.4863e-006 9.7284e-004 3.0466e-005 4.7160e-006 1.9122e-004 6.4133e-005 2.3158e-006 4.4189e-004 6.4657e-005 8.1697e-006

3.5243e-003 9.1461e-003 1.2058e-004 5.4639e-003 6.7642e-004 5.5088e-004 7.3689e-003 9.9694e-003 8.6171e-004 7.0083e-003 3.9111e-004 2.8687e-004

1.7412 1.7685 1.7878 1.5681 1.5746 1.5910 1.3776 1.3811 1.3912 1.1871 1.1884 1.1895

1.9907 1.9923 1.9974 1.9890 1.9952 1.9997 1.9949 1.9975 2.0010 1.9991 1.9998 2.0013

0.4

0.6

0.8

6. Numerical experiments In this section, efficiency of proposed difference scheme is demonstrated in terms of computational error and order of convergence for spatial and temporal dimensions through a few examples for the problem, (2)–(5). All numerical experiments are done with help of MATLAB software. Example 1 [35]. Consider problem (2)–(5) with non-zero initial conditions.

∂ α u(x, t ) ∂ 4 u(x, t ) = ( x2 + 1 ) + exp(x ) u(x, t ) + u(x, t )3 + u(x, t − 0.2 ) + g(x, t ), 0 < x < 1, 0 ≤ t < 1, α ∂t ∂ x4   (α + 4 ) here g(x, t ) = (1 − 16π 4 x2 − exp(x )) sin(2π x )t α +3 + t 3 sin(2π x ) + 16π 4 t α 6

− sin(2π x )(t − 0.2 )α +3 − (sin(2π x )t α +3 )3 , Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

11

Table 4 The computational error and convergence order in spatial dimension for Example 2 by varying h.

α

h

Using [36] EL2 (h, τ )

Using Eq. (35) EL2 (h, τ )

Using Eq. (35) Order(h)

Using [36] EL∞ (h, τ )

Using Eq. (35) EL∞ (h, τ )

Using Eq. (35) Order(h)

0.2

1/20 1/50 1/80 1/20 1/50 1/80 1/20 1/50 1/80 1/20 1/50 1/80

1.8163e-006 7.0584e-007 2.3643e-007 5.2408e-006 3.4694e-007 2.9722e-007 8.2390e-006 3.5841e-007 6.8762e-007 3.2333e-006 4.9105e-007 1.7621e-007

4.1952e-006 2.4260e-007 5.7621e-008 6.2360e-006 1.4786e-007 6.8248e-008 7.4628e-006 2.0355e-007 7.8321e-008 8.9624e-006 3.7326e-007 1.4601e-008

3.9866 3.9912 3.9973 3.9889 3.9917 3.9980 4.0000 3.9984 4.0003 3.9985 3.9996 4.0010

7.0172e-006 2.3622e-006 5.8567e-007 2.3644e-006 8.0037e-007 1.7638e-007 9.9213e-006 3.3894e-007 2.4628e-007 4.1899e-006 6.4627e-007 8.3961e-007

6.2481e-006 3.4612e-006 4.7811e-008 1.1293e-006 2.2645e-007 6.3694e-008 5.4687e-006 2.3691e-007 7.7350e-008 4.5017e-006 2.4631e-007 5.7912e-008

3.9885 3.9924 3.9965 3.9888 3.9939 4.9989 3.9948 3.9987 4.0009 3.9983 3.9996 4.0011

0.4

0.6

0.8

(a)

(b)

Numerical solution

Exact solution

1

0

-1 1

1

0

-1 1

1 0.5

t-axis

1 0.5

0.5 0 0

x-axis

t-axis

0.5 0 0

Fig. 1. Exact and numerical solution surface with α = 0.3, h = τ =

u(x, t ) = sin(2π x )t α +3 , u ( 0, t ) = 0,

u ( 1, t ) = 0,

0 ≤ x ≤ 1,

1 40

x-axis for Example 1.

t ∈ [−0.2, 0],

∂ u ( 0, t ) = 0, ∂ x2 2

∂ 2 u ( 1, t ) = 0, 0 ≤ t ≤ 1. ∂ x2

The exact solution of the problem is u(x, t ) = sin(2π x )t α +3 . We have solved the above problem with proposed scheme (35)–(39). The efficiency of the scheme is numerically examined by taking sufficiently small spatial and temporal steps. The results illustrate that our scheme has a spatial accuracy of O(h4 ) and temporal convergence of order O(τ 2 ). A good agreement of theoretical and numerical  results is  obtained. Let EL2 (h, τ ) = max1≤n≤N ||un − U n ||, E∞ (h, τ ) = max1≤n≤N ||un − U n ||∞ , Order (τ ) = log2



log2

EL (2h,τ ) 2

EL (h,τ )



EL (h,2τ ) 2

EL (h,τ ) 2

and Order (h ) =

.

2

1 First we fix h = 40 and keep varying τ . In Table 1, error is presented in L2 norm E2 (τ , h) and L∞ norm E∞ (τ , h). It is observed that temporal accuracy of order O(τ 2 ) is achieved conforming the theoretical results. 1 Next we fix τ = 40 and keep varying h. In Table 2, error is presented in L2 norm E2 (τ , h) and L∞ norm E∞ (τ , h). It is observed that spatial accuracy of order O(h4 ) is achieved conforming the theoretical results. 1 For visualization, in Fig. 1a and b we plot the exact numerical solution considering τ = h = 40 indicating an excellent closeness between numerical and exact solution. Error comparison of scheme using L1-approximation of Caputo fractional derivative is also represented in Tables 1 and 2 for temporal and spatial dimensions.

Example 2 [33]. Consider problem (2)–(5) with non-zero initial and boundary conditions as follows

∂ α u(x, t ) ∂ 4 u(x, t ) = ( x3 + x ) + u(x, t )2 + u(x, t − 0.1 ) + xu(x, t ) + g(x, t ), 0 < x < 1, 0 ≤ t < 1, ∂tα ∂ x4 here

g(x, t ) = 10x3 (2 − 3x )(t 4 + 1 ) +

4!

(5 − α )

t 7/2 x5 (1 − x ) + 10x4 (t 4 + 1 ) − 20(t 4 + 1 )x3 (1 − x )

Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC 12

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

(a)

(b)

Numerical solution

Exact solution

0.15 0.1 0.05 0 1

0.15 0.1 0.05 0 1

1 0.5 0 0

t-axis

1 0.5

0.5

x-axis

0.5 0 0

t-axis

Fig. 2. Exact and numerical solution surface with α = 0.3, h = τ =

x-axis for Example 2.

(b)

4 2 0 1

1

0.5

t-axis

Numerical Solution

Exact solution

(a)

1 40

4

2

0 1

0 0

1

0.5

0.5

x-axis

t-axis

Fig. 3. Exact and numerical solution surface with α = 0.5, h = τ =

0.5 0 0 1 40

x-axis

for Example 3.

− (x3 + x )(120x − 360x2 )(t 4 + 1 ) − x5 (1 − x )((t − 0.1 )4 + 1 ) − (x5 (1 − x )(t 4 + 1 ))2 , u(x, t ) = (x5 − x6 )(t 4 + 1 ), u ( 0, t ) = 0,

u ( 1, t ) = 0,

0 ≤ x ≤ 1, 0 ≤ t ≤ 1,

t ∈ [−0.1, 0],

∂ 2 u ( 0, t ) = 0, ∂ x2

∂ 2 u ( 1, t ) = −10(t 4 + 1 ), 0 ≤ t ≤ 1. ∂ x2

The exact solution of the problem is u(x, t ) = (x5 − x6 )(t 4 + 1 ). We have solved the above problem with proposed scheme (35)–(39). The efficiency of the scheme is numerically examined by taking sufficiently small spatial and temporal steps. The results illustrate that our scheme has spatial convergence of order O(h4 ) and temporal convergence of order O(τ 2 ). A good agreement of theoretical and numerical results is obtained. 1 First, we fix h = 40 and keep varying τ . In Table 3, error is presented in L2 norm E2 (τ , h) and L∞ norm E∞ (τ , h). It is observed that temporal accuracy of order O(τ 2 ) is achieved conforming the theoretical results. 1 Next, we fix τ = 40 and keep varying h. In Table 4, error is presented in L2 norm E2 (τ , h) and L∞ norm E∞ (τ , h). It is observed that spatial accuracy of order O(h4 ) is achieved conforming the theoretical results. 1 For visualization, in Fig. 2a and b we plot the exact and numerical solution considering τ = h = 40 indicating an excellent closeness between numerical and exact solution. Numerical results in Tables 2 and 4 show that for the same spatial step size, the scheme given in Eq. (35) is slightly more accurate than the scheme in [36]. Although both of these schemes give fourthorder convergence in spatial dimensions; which implies the advantage of our method. Error improvement from Tables 1 and 3 can be observed. Therefore, overall our method has an advantage in temporal convergence order and error improvement. We now give an example with non-zero initial and boundary conditions. Example 3 [18]. Consider problems (2)–(5) with non-zero initial and boundary conditions.

∂ α u(x, t ) ∂ 4 u(x, t ) =x + exp(x ) u(x, t ) + u(x, t − 0.5 ) + g(x, t ), 0 < x < 1, 0 ≤ t < 1, α ∂t ∂ x4  (2 + α ) 2α  here g(x, t ) = exp(x ) (1 + α )t α − t (1 + 2α ) Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

ARTICLE IN PRESS

JID: AMC

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

13

Table 5 The computational error and convergence order in temporal dimension for Example 3 by varying τ .

α

τ

Using Eq. (35) EL2 (h, τ )

Using [36] EL2 (h, τ )

Using [36] Order(τ )

Using Eq. (35) Order(τ )

Using Eq. (35) EL∞ (h, τ )

Using [36] EL∞ (h, τ )

Using [36] Order(τ )

Using Eq. (35) Order(τ )

0.3

1/20 1/50 1/80 1/20 1/50 1/80 1/20 1/50 1/80 1/20 1/50 1/80

3.5361e-004 9.2030e-005 2.1984e-006 4.4936e-004 6.7285e-005 1.6917e-006 8.3420e-004 3.0429e-005 6.6437e-006 4.8511e-004 9.6903e-005 1.7003e-006

7.2113e-003 1.6171e-003 3.3892e-004 6.7328e-003 2.4299e-003 8.2540e-004 4.0912e-003 9.4637e-004 2.7300e-004 5.4928e-003 3.3967e-004 9.0666e-004

1.7688 1.7881 1.7923 1.5787 1.5879 1.5923 1.3430 1.3652 1.3890 1.1581 1.1774 1.1835

1.9768 1.9828 1.9964 1.9836 1.9907 1.9987 1.9844 1.9953 1.9996 1.9867 1.9923 2.0023

2.6135e-004 7.8221e-005 1.3946e-006 9.2879e-004 3.5173e-005 5.2985e-006 6.1360e-004 9.7618e-005 2.2493e-006 4.6184e-004 6.8061e-005 8.4792e-006

2.4623e-003 4.0249e-003 7.1762e-004 1.9004e-003 3.8711e-004 9.6173e-004 2.4632e-003 5.5617e-004 4.9243e-005 8.4682e-003 2.1338e-004 5.5527e-004

1.7726 1.7888 1.7929 1.5831 1.5891 1.5973 1.3529 1.3726 1.3894 1.1630 1.1789 1.1864

1.9863 1.9910 1.9973 1.9825 1.9938 1.9981 1.9937 1.9971 2.0017 1.9975 1.9991 2.0011

0.5

0.7

0.9

Table 6 The computational error and convergence order in spatial dimension for Example 3 by varying h.

α

h

Using [36] EL2 (h, τ )

Using Eq. (35) EL2 (h, τ )

Using Eq. (35) Order(h)

Using [36] EL∞ (h, τ )

Using Eq. (35) EL∞ (h, τ )

Using Eq. (35) Order(h)

0.3

1/10 1/20 1/40 1/10 1/20 1/40 1/10 1/20 1/40 1/10 1/20 1/40

4.9621e-006 7.4628e-007 1.1341e-007 6.5894e-006 8.3478e-007 3.6930e-007 5.5281e-006 2.3943e-006 9.8247e-007 4.3903e-006 6.5583e-007 2.6731e-008

8.1124e-006 2.3872e-007 5.6328e-008 6.7231e-006 8.0681e-007 1.3462e-008 4.2843e-006 7.5757e-007 2.3699e-008 8.4800e-006 3.9013e-007 9.8367e-008

3.9816 3.9877 3.9935 3.9834 3.9882 3.9949 3.9831 3.9907 3.9969 3.9853 3.9948 3.9987

3.3924e-006 5.2516e-007 2.3981e-007 6.7634e-006 1.4862e-007 8.6149e-007 3.2982e-006 2.3747e-007 9.0010e-007 5.3421e-006 6.9166e-007 1.5386e-007

5.3005e-006 2.8221e-007 7.3627e-008 8.6992e-006 2.7385e-007 4.4623e-008 9.3762e-006 2.7764e-007 1.9111e-008 3.3905e-006 6.3317e-007 7.6876e-008

4.0038 3.9911 3.9905 4.0011 3.9916 3.9900 4.0024 3.9909 3.9879 4.0003 3.9913 3.9829

0.5

0.7

0.9

u(x, t ) = exp(x )t 1+α , u(0, t ) = t 1+α ,

0 ≤ x ≤ 1,

u(1, t ) = e t 1+α ,

t ∈ [−0.5, 0],

∂ 2 u ( 0, t ) = t 1+α , ∂ x2

∂ 2 u ( 1, t ) = e t 1+α , 0 ≤ t ≤ 1. ∂ x2

The exact solution of the problem is u(x, t ) = exp(x )t 1+α . 1 In Fig. 3a and b we plot the exact and numerical solution considering τ = h = 40 and more numerical results are presented in Tables 5 and 6 providing temporal and spatial convergence orders along with error, respectively.

7. Conclusion In the present paper, we constructed a compact difference scheme for a class of non-linear time-delayed fractional subdiffusion wave equation with variable coefficients. Time fractional term is approximated using L2 − 1σ which improves the temporal accuracy up-to order O(τ 3−α ) and spatial order accuracy to fourth-order using compact difference operator. The non-linear source term is linearized using Taylor’s series. The proposed numerical scheme is more important than the ones proposed earlier by many authors. Unique solvability is attained, stability and convergence are also proved using a discrete energy method for the proposed compact difference scheme. Numerical experimentation is done and an excellent agreement is found between numerical and analytical values of given examples. Acknowledgment We thank the Ministry of Human Resource Development, Government of India for its financial support (grant no. OH-3123-200-428). References [1] R. Klages, G. Radons, I.M. Sokolov, Anomalous Transport: Foundations and Applications, Wiley-VCH, Weinheim, 2008. [2] F. Mainardi, M. Raberto, R. Gorenflo, E. Scalas, Fractional calculus and continuous-time finance II: the waiting-time distribution, Physica A 287 (20 0 0) 468–481. [3] D. Benson, R. Schumer, M. Meerschaert, S. Wheatcraft, Fractional dispersion, Levy motion, and the MADE tracer tests, Transp. Porous Media 42 (2001) 211–240.

Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900

JID: AMC 14

ARTICLE IN PRESS

[m3Gsc;November 19, 2019;3:44]

S. Nandal and D. Narain Pandey / Applied Mathematics and Computation xxx (xxxx) xxx

[4] I. Podlubny, Fractional Differential Equations, Academic Press, New York, 1999. [5] A. Kilbas, H. Srivastava, J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier Science and Technology, 2006. [6] W. Gu, Y. Zhou, X. Ge, A compact difference scheme for solving fractional neutral parabolic differential equation with proportional delay, Hindawi J. Funct. Spaces 2017 (2017) 8. Article ID 3679526 [7] Y.M. Wang, T. Wang, A compact ADI method and its extrapolation for time fractional sub-diffusion equations with nonhomogeneous Neumann boundary conditions, Comput. Math. Appl. 75 (3) (2018) 721–739. [8] G.H. Gao, Z.Z. Sun, A compact difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. [9] T.A.M. Langlands, B.I. Henry, The accuracy and stability of an implicit solution method for the fractional diffusion equation, J. Comput. Phys. 205 (2005) 719–736. [10] J.C. Ren, Z.Z. Sun, X. Zhao, Compact difference scheme for the fractional sub-diffusion equation with Neumann boundary conditions, J. Comput. Phys. 232 (2013) 456–467. [11] S. Vong, Z. Wang, High order difference schemes for a time-fractional differential equation with Neumann boundary conditions, East Asian J. Appl. Math. 4 (3) (2014) 222–241. [12] Y.N. Zhang, Z.Z. Sun, H.L. Liao, Finite difference methods for the time fractional diffusion equation on non-uniform meshs, J. Comput. Phys. 265 (2014) 195–210. [13] C. Chen, F. Liu, V. Anh, I. Turner, Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equations, SIAM J. Sci. Comput. 32 (4) (2010) 1740–1760. [14] Y.N. Zhang, Z.Z. Sun, H.W. Wu, Error estimates of Crank-Nicolson type difference schemes for the subdiffusion equation, SIAM J. Numer. Anal. 49 (6) (2011) 2302–2322. [15] Q. Zhang, C. Zhang, A compact difference scheme combined with extrapolation techniques for solving a class of neutral delay parabolic differential equations, Appl. Math. Lett. 26 (2) (2013) 306–312. [16] W. Gu, A compact difference scheme for a class of variable coefficient quasilinear parabolic equations with delay, Abstract Appl. Anal. 2014 (2014) 8. Article ID 810352 [17] D. Li, C. Zhang, J. Wen, A note on compact finite difference method for reaction-diffusion equations with delay, Appl. Math. Model. 39 (2015) 1749–1754. [18] F. Liu, C. Yang, K. Burrage, Numerical method and analytical technique of the modified anomalous subdiffusion equation with a nonlinear source term, J. Comput. Appl. Math. 231 (2009) 160–176. [19] M. Cui, Compact finite difference method for the fractional diffusion equation, J. Comput. Phys. 228 (2009) 7792–7804. [20] M.M. Khader, T.S. Danaf, A.S. Hendy, A computational matrix method for solving systems of high order fractional differential equations, Appl. Math. Model. 37 (2013) 4035–4050. [21] F. Mainardi, G. Pagnini, A. Mura, R. Gorenflo, Time-fractional diffusion of distributed order, J. Vib. Control 14 (2008) 1267–1290. [22] M.L. Morgado, M. Rebelo, Numerical approximation of distributed order reaction-diffusion equations, J. Comput. Appl. Math. 275 (2015) 216–227. [23] N.J. Ford, M.L. Morgado, Distributed order equations as boundary value problems, Comput. Math. Appl. 64 (2012) 2973–2981. [24] J.T. Katsikadelis, Numerical solution of distributed order fractional differential equations, J. Comput. Phys. 259 (2014) 11–22. [25] G.H. Gao, Z.Z. Sun, Two alternating direction implicit difference schemes with the extrapolation method for the two-dimensional distributed-order differential equations, Comput. Math. Appl. 69 (2015) 926–984. [26] G.H. Gao, H.W. Sun, Z.Z. Sun, Some high order difference schemes for distributed-order differential equations, J. Comput. Phys. 289 (2015) 337–359. [27] G. Gao, Z.Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. [28] A.A. Alikhanov, A new difference scheme for the time fractional diffusion equation, J. Comput. Phys. 280 (2015) 424–438. [29] S. Nandal, D.N. Pandey, Numerical solution of time fractional non-linear neutral delay differential equations of fourth-order, Malaya J. Math. 7 (3) (2019) 579–589. [30] Q. Zhang, M. Ran, D. Xu, Analysis of compact difference scheme for the semilinear fractional partial differential equation with time delay, Appl. Anal. 96 (11) (2017) 1867–1884. [31] P. Zhang, H. Pu, A second-order compact difference scheme for the fourth-order fractional sub-diffusion equation, Numer. Algorithms 76 (2017) 573–598. [32] 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 (2009) 1760–1781. [33] I. Karatay, S.R. Bayramoglu, High-order compact difference scheme for the numerical solution of time fractional heat equations, Sci. World J. 2014 (2014) 8. Article ID 642989 [34] A.A. Samarskii, V.B. Andreev, Difference Methods for Elliptic Equations (Russian Book), Izdatel’stvo Nauka, Moscow, 1976, p. 352. [35] Z. Hao, W. Cao, G. Lin, A second-order difference scheme for the time fractional substantial diffusion equation, J. Comput. Appl. Math. 313 (2017) 54–69. [36] X. Hu, L. Zhang, A new implicit compact difference scheme for the fourth-order fractional diffusion-wave system, Int. J. Comput. Math. 91 (10) (2014) 2215–2231. [37] V.G. Pimenov, A.S. Hendy, R.h. de Staelen, On a class of non-linear delay distributed order fractional diffusion equations, J. Comput. Appl. Math. 318 (2017) 433–443. [38] N.H. Sweilam, S.M. Al-Mekhlafi, A.O. Albalawi, A novel variable-order fractional nonlinear Klein Gordon model: a numerical approach, Numer. Methods Partial Differ. Equ. 35 (5) (2019) 1617–1629. [39] N.H. Sweilam, S.M. Al-Mekhlafi, A novel numerical method for solving the 2-d time fractional cable equation, Eur. Phys. J. Plus 134 (7) (2019) 323.

Please cite this article as: S. Nandal and D. Narain Pandey, Numerical solution of non-linear fourth order fractional subdiffusion wave equation with time delay, Applied Mathematics and Computation, https://doi.org/10.1016/j.amc.2019.124900