Numerical treatment of non-linear fourth-order distributed fractional sub-diffusion equation with time-delay

Numerical treatment of non-linear fourth-order distributed fractional sub-diffusion equation with time-delay

Journal Pre-proof

Numerical treatment of non-linear fourth-order distributed fractional sub-diffusion equation with time-delay Sarita Nandal, Dwijendra Narain Pandey PII: DOI: Reference:

S1007-5704(19)30465-4 https://doi.org/10.1016/j.cnsns.2019.105146 CNSNS 105146

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Received date: Revised date: Accepted date:

14 February 2019 29 November 2019 3 December 2019

Please cite this article as: Sarita Nandal, Dwijendra Narain Pandey, Numerical treatment of non-linear fourth-order distributed fractional sub-diffusion equation with time-delay, Communications in Nonlinear Science and Numerical Simulation (2019), doi: https://doi.org/10.1016/j.cnsns.2019.105146

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Highlights • Please find herewith a paper, entitled by “Numerical treatment of non-linear fourth order distributed fractional sub-diffusion equation with time-delay”. This paper is concerned with producing high convergence order for fourth-order distributed fractional sub-diffusion time-delay differential equation. For increasing order of convergence in time, space and fractional order dimension, we opted for L2-1σ formula, compact difference operator and Simpson’s 1/3rd rule respectively. This gives higher convergence order in both time and space and better than the numerical techniques present in literature. We did numerical experimentation to show the efficiency and accuracy of the proposed scheme. Our constructed scheme is unconditionally stable and convergent. In remark section, we also gave an extension of the problem where proposed scheme can be implemented very efficiently. • Fourth order spatial derivatives are required in fractional diffusion equation to model groove formation on metal surface with a beam boundary and in beam transverse vibration. We included the effect of time-delay as in most of the real-life models there is always a delay which makes the model more real-life based and accurate. Additionally, for the solution of the considered equation, we opted for compact difference technique which gives the highest convergence order than any other techniques developed so far for fractional differential equations. Therefore, this article gives an efficient numerical technique for a real-world based problem.

NUMERICAL TREATMENT OF NON-LINEAR FOURTH-ORDER DISTRIBUTED FRACTIONAL SUB-DIFFUSION EQUATION WITH TIME-DELAY

Sarita Nandal1 , Dwijendra Narain Pandey2 Department of Mathematics Indian Institute of Technology Roorkee Roorkee-247667, India 1 [email protected], 2 [email protected] 1,2

Abstract In this paper, we proposed a linearized second-order numerical technique for non-linear fourth-order distributed fractional sub-diffusion equation with time delay. Time fractional derivative is represented using Caputo derivative and further approximated using L2 − 1σ formula which gives secondorder temporal convergence and compact difference operator is employed for spatial dimensions. The proposed method is unconditionally stable and convergent to the analytical solution with the order of convergence O(τ 2 + h4 + (∆α)4 ) improvement in earlier work. Non-linear terms are linearized with the help of Taylor’s series. At the last, we provided a few examples to show the efficiency of the compact difference scheme to support the theoretical results and also presented comparison with L1-approximation of Caputo fractional derivative. Key words: Distributed fractional derivative, L2 − 1σ formula, Compact difference scheme, Stability, Convergence. 1. Introduction Nowadays there is rapid growth in the study of fractional differential equations in the formation of various kinds of mathematical models such as in diffusion, finance, material science, etc. Many recently published articles have considered a substantial number of such models for the input-output relationship of a time variable rooted in frequency-domain [1]- [3]. Fractional distributed order differential equations are used widely in the model of elastic mediums and stress-strain behaviors to discuss eigenfunctions for torsional model/dielectric spherical shells which further leads to fractional distributed-order differential equations [4]- [6]. Hartley and Lorenzo [8] used distributed-order differential equations to model thermeorheological behavior. Atanakovic et al [11] studied and analyzed the existence and uniqueness of solutions for fractional distributed-order differential equations having vital applications in viscoelasticity and system identification theory. Diethelm and Ford [12] numerically studied distributed order differential equations by converting it to multi-term fractional differential equations and proved the convergence of the proposed scheme. Najafi et al [13] discussed the stability and convergence of first-order distributed fractional differential equations. Katsikadelis [14] developed the numerical scheme for fractional distributed order (FDO) differential equations and approached the FDO with two different approaches, in one approach by converting to a single FDO differential equation and in the second approach by considering the system of FDO differential equations. Sun and Zhang [15] developed a linearized numerical technique for nonlinear delay differential equations, proved the stability, unique solvability, and convergence of order O(τ 2 +h4 ) in L∞ norm. Katsikadelis [16] constructed the difference scheme for both linear and non-linear FDO differential equations by approximating distributed order derivative term with multi-term fractional derivatives. Hao et al [18] developed the fourth-order numerical approximation scheme for fractional differential equations for 1D and 2D both and discussed stability and convergence as well. Li et al [19] investigate the long-term stability of compact difference scheme for the delayed reaction-diffusion equation. Ye et al [20] derived a compact difference scheme for a distributed time-fractional differential equation. Gao et al [21] 2

3

developed two difference schemes for both 1D and 2D and proved unconditional stability and convergence using L1 and L∞ norm with following convergence orders O(τ 2 + h2 + (∆α)2 ) and O(τ 2 + h4 + (∆α)4 ). Deng [22] derived a high order compact alternating direct implicit difference scheme, which is a combination of fourth-order compact difference approximation to spatial derivatives and second-order backward difference formula for temporal derivatives. Authors in [31] and [32] studied variable order fractional differential equation and neutral delay differential equation of fourth-order and constructed a scheme to enhance the order of convergence for time-dimension. Gao and Sun [23] developed two difference schemes with unconditional stability and convergence of orders O(τ 2 + h21 + h22 + (∆α)2 ) and O(τ 2 + h41 + h42 + (∆α)4 ) for two-dimensional distributed order fractional differential equations. Hao et al [24] proposed a linearized compact difference numerical scheme for semi-linear fractional diffusion-wave equation with time delay. Li and Wu [25] proposed a numerical scheme for solving distributed order fractional diffusion equations with the use of quadrature formulas by converting distributed fractional term into the multi-term fractional derivative. Authors in [26] - [10] studied both the fourth-order and second-order fractional diffusion equations and discussed the analysis of the proposed numerical scheme also. Authors in [33] studied and analyzed the mixed finite element Galerkin method in combination with a discrete scheme for the fourth-order diffusion equation. Fourthorder fractional diffusion equations are very important in model formation for beam vibration and groove formation on a metal surface. There is some innovating work done in [38]- [40] for the fourth-order fractional diffusion equation given in terms of Mittag-Leffler function. There are numerous realistic models based on distributed order fractional diffusion-wave equation. Authors in [44] considered the constitutive equation of a type of distributed fractional-order for the study of waves in a viscoelastic rod of finite length. They prescribed boundary conditions on displacement and obtained stress and displacement in a stress-relaxation test. An author in [45] introduced new equations based on distributed order fractional Fokker-Planck equation of flow in swelling porous media in a material coordinate for infiltration and absorption. For the study of radial groundwater flow to or from well, authors in [46] introduced the theory of distributed order fractional-order diffusion-wave equation and also discussed its application to pumping and slug tests. One of the most exciting recent advances of nonlinear science and theoretical physics has been a development of methods to look for exact solutions for nonlinear differential equations. This is important because many mathematical models are described by nonlinear differential equations. The modified simple equation method is an efficient method for obtaining exact solutions of nonlinear evolution equations. In the paper [43], the authors proposed the modified simple equation method is applied to construct exact solutions of the modified equal width (MEW) equation and the Fisher equation and the nonlinear Telegraph equation and the Cahn-Allen equation. Authors in [49] presented an operator-based framework for the construction of analytical soliton solutions to fractional Riccati equations. To preserve the additivity of base function powers under multiplication, fractional differential equations are mapped from Caputo algebra to Riemann-Liouville algebra. An author in [50] derived a new semi-analytical approach based on Fourier-series expansion for solving a wide class of fractional partial differential equations (FPDEs) such as wave-diffusion equations, modified anomalous fractional sub-diffusion equations, time-fractional telegraph equations. Authors in [53] derived the solutions for the KdV equation by the algebraic operator method based upon the generalized operator of differentiation (has been exploited for the derivation of analytical solutions for KdV equation). The inverse scattering transforms [41] and the direct method by Hirota [42] are known as impressive methods to search solutions of exactly solvable differential equations. The singular manifold method [43], the transformation method, the tanh-function method, the sh-function method and the Weierstrass function method are useful in many applications to search exact solutions of nonlinear partially solvable differential equations. Overall, there are various methods used to find analytical/semi-analytical solution of the various partial differential equations by numerous method listed in [42], [50] [51], [52] and the references therein. As there are only very few partial differential equations with or without delay whose analytical solutions can be expressed. The existence of time-delay term with non-linear source term makes it difficult to carry the theoretical analysis of time-delay partial

4

differential equations, while the numerical method can greatly make up for the lack of theoretical research. Therefore, in this manuscript, we opt to construct a numerical scheme for the below-considered equation. In the present article, we are interested in constructing a linearized numerical technique of improved time-dimension convergence order to solve the following fourth-order distributed-fractional sub-diffusion differential equation with time delay of the following form Z 1 ∂ 4 u(x, t) ∂ α u(x, t) dα + = u(x, t) + f (x, t, u(x, t), u(x, t − s)), (x, t) ∈ (0, L) × (0, T ), (1) w(α) ∂tα ∂x4 0 Subject to following initial conditions and boundary conditions u(x, t) = φ(x, t), u(0, t) = α1 (t), 2

(x, t) ∈ [0, L] × [−s, 0],

u(L, t) = α2 (t),

∂ u(0, t) = β1 (t), ∂x2

2

t ∈ [0, T ],

∂ u(L, t) = β2 (t), ∂x2

t ∈ [0, T ].

(2) (3) (4)

α

u(x,t) where ∂ ∂t is the fractional derivative of order α ∈ (0, 1) defined in the sense of Caputo, w(α) > 0 α R1 [30] is the function acting as weight for the order of differentiation and 0 w(α) dα = W (a constant) > 0, s > 0 is the time delay, φ(x, t), α1 (t), α2 (t), β1 (t), β2 (t) are all given smooth functions and f (x, t, u(x, t), u(x, t − s)) is the non-linear source term with time-delay and sufficiently smooth. Lipschitz continuity of source function is also considered. Consider that the partial derivatives fµ (x, t, µ, ν) and fν (x, t, µ, ν) are continuous in the 0 neighborhood of the solution, where 0 is a positive constant. Define c0 = max |u(x, t)|, 0
c1 = c2 =

sup 0
sup

0
|fµ (x, t, u(x, t) + 1 , u(x, t − s) + 2 )|, |fν (x, t, u(x, t) + 1 , u(x, t − s) + 2 )|.

8,3 We also assume that the function u(x, t) ∈ Cx,t ([0, L] × [0, T ]) i.e. eight times continuously differentiable for spatial variable x and three times continuously differentiable for temporal variable t. Assume m be the integer satisfying ms ≤ T ≤ (m + 1)s. Define Ir = (rs, (r + 1)s), r = −1, 0, ..., m − 1, Im = (ms, T ), Sm I = q=−1 Iq . In equation (1), we have three types of complexities as follows: first is time-delay term, second is non-linear source term and third is distributed fractional derivative (i.e. integral term on the left-hand side (L.H.S.) of equation (1). Here we propose to construct a compact difference scheme using L2 − 1σ [9] giving a best second-order temporal convergence for (0 < α < 1). We opted for Simpson’s rule for converting distributed-fractional derivative term to the multi-fractional derivative. In the numerical experimentation section, we provided a comparison with L1-approximation [29] (a discrete approximation of Caputo fractional derivative). In this article, we aim to develop an efficient numerical technique that has improved the time-dimension convergence order to O(τ 2 ) for (0 < α < 1).

2. Preliminaries and Notations In this section, we express some basic definitions and lemmas to be used in the sequel. Here, we express the most standard definitions of fractional integral and Caputo fractional derivative, on a finite interval of the real line. Definition 1. Fractional Integral: The fractional integral [35] of order α for a function u is defined by Z x 1 J α u(x) = (x − s)α−1 u(s)ds, x > 0, α > 0, (5) Γ(α) a provided the integral in the right hand side exists. Here Γ(·) is the Gamma function.

5

Definition 2. For a function u(x), the Riemann-Liouville (RL) fractional derivative [34]- [35] of order α is given by Rx u(s) ( 1 dm Γ(m−α) dxm a (x−s)α−[α] ds, m − 1 ≤ α < m, m ∈ N, RL α a Dx u(x) = dm u α = m ∈ N. d xm , Definition 3. The Caputo fractional derivative [37] of order α for a function u ∈ C m ([0, ∞)) is defined by C α 0 Dx u(x)

:= J m−α Dm u(x) ( Rx 1 m−α−1 (m) u (s)ds, Γ(m−α) 0 (x − s) := m d u(x), m dt

m − 1 ≤ α < m; α = m.

(6)

For more details on fractional calculus reader may follow [35]- [37]. The below-given Lemma will be used for approximation of integral term in equation (1) and hence converting the distributed order fractional derivative to multi-term fractional derivatives. Lemma 1. ( [20], [30]) Consider equidistant division of the interval [0, 1] into 2J subintervals, assume 1 ∆α = 2J and denote αl = l∆α, 0 ≤ l ≤ 2J. Then the composite Simpson’s rule reads Z 1 2J X (∆α)4 (4) G (ξ), ξ ∈ [0, 1]. (7) G(α)dα = ∆α γl G(αl ) − 180 0 l=0

γl =

(

1 3 2 3 4 3

l = 0, 2J l = 2, 4, ...2J − 4, 2J − 2 l = 1, 3, ...2J − 3, 2J − 1.

For approximation, we consider partition of the mesh Ωh t = Ωh ×Ωt where Ωh = {xj = j h|0 ≤ j ≤ M } L and Ωt = {tk = kτ, k = 0, 1, ..., N ; T = τ N }. Introducing two positive integers M and N , let h = M , T τ = N are the space and temporal step lengths respectively. Denote Vh = {u|u = (u0 , u1 , u2 , ..., uM ), u0 = uM = 0} as the grid function space on Ωh . For any u ∈ Vh we define δx uj = h1 (uj − uj−1 ), δx2 uj = 1 2 2 4 h2 (uj+1 − 2uj + uj−1 ), δx uj = δx (δx uj ). We now introduce definition given by Alikhanov [9] (discrete approximation of Caputo fractional derivative) and a few Lemmas which are important in the construction of the numerical scheme for equations (1)-(4). Definition 4. [9] Let σ = 1 − α2 , then for the Caputo fractional derivative of the order α, 0 < α < 1, of the function u(t) ∈ C 3 [0, T ] at the fixed point tk−1+σ , k ∈ {1, 2, ..., N }, a discrete approximation for Caputo-fractional derivative called L2 − 1σ formula of second order temporal convergence for α ∈ (0, 1) is given as follows: Let a0 = σ 1−α , al = (l + σ)1−α − (l − 1 + σ)1−α , l ≥ 1, 1 bl = 2−α [(l + σ)2−α − (l − 1 + σ)2−α ] − 12 [(l + σ)1−α + (l − 1 + σ)1−α ], l ≥ 1, (k)

when k = 0, c0 = a0 , and further when k ≥ 1,

(k) cj

Defining C α 0 Dt

uk−1+σ j

=

(

a0 + b1 , j = 0, aj + bj+1 − bj , 1 ≤ j ≤ k − 1, aj − bj , j = k.

" # k−1 X (k) τ 1−α (k) k (k) (k) 0 j = c u − (ck−j−1 − ck−j )u − ck−1 u . Γ(2 − α) 0 j=1

(8)

(9)

Below given Lemma estimates the error of the L2 − 1σ formula (a discrete approximation of Caputo fractional derivative).

6

Lemma 2. [9] For any α ∈ (0, 1) and u ∈ C 3 [0, tk+1 ], it holds that C α C α k−1+σ = O(τ 3−α ). 0 Dt u(t)|t=tk−1+σ − 0 Dt uj (k)

Below given Lemma 3 provide properties of constants cj occurred in L2 − 1σ formula and, Lemma 4 provides a property for L2 − 1σ formula which will play a vital role in the analysis of the proposed numerical scheme. (k)

Lemma 3. [9] Assume α ∈ (0, 1). cj (0 ≤ j ≤ k − 1) is defined by equation (8), then, the following holds 1−α (k) cj > (k + σ)−α , (10) 2 (k) (k) (k) (k) (k) ck−1 < ck−2 < ... < c2 < c1 < c0 . (11) Lemma 4. [9] Assume a grid function u = {uk |0 ≤ k ≤ N } defined on Ωt , then the following holds
1C α α D (u2 ) ≤ σuk + (1 − σ)uk−1 C (12) 0 Dtk−1+σ u. 2 0 tk−1+σ Next Lemma gives the truncation error of fourth order for the linear operator (to be implemented for spatial dimensions i.e for variable x) as described in definition 5. Lemma 5. [61], Consider that function g(x) ∈ C 6 [0, L]. It holds that h4 (6) g 00 (xj+1 ) + 10g 00 (xj ) + g 00 (xj−1 ) − δx2 gj = g (ξ), 12 240 g(xj+1 ) − 2g(xj ) + g(xj−1 ) where δx2 gj = , h2 and ξ ∈ [xj−1 , xj+1 ].

1 ≤ j ≤ M − 1, (13)

Below given definition introduces the linear operator to be implemented for obtaining fourth-order convergence in spatial dimension for equation (1). Definition 5. For g = (g0 , g1 , ..., gM ) the linear operator [61] is given as follows ( 1 1≤j ≤M −1 12 (gj+1 + 10gj + gj−1 )
v(xj , tk ) = Vjk ,

0 ≤ j ≤ M, −n ≤ k ≤ N.

(14)

Now, we consider the derivation of finite difference scheme. For the sake of simplifying the equation 2 u(x,t) , then our problem (1) − (4) becomes (1), we use the reduced order method and suppose v(x, t) = ∂ ∂x 2 Z 1 ∂ α u(x, t) ∂ 2 v(x, t) w(α) dα + = u(x, t) + f (x, t, u(x, t), u(x, t − s)), 0 < x < L, 0 < t < T, (15) ∂tα ∂x2 0 ∂ 2 u(x, t) v(x, t) = , 0 < x < L, 0 < t < T, (16) ∂x2 u(x, t) = φ(x, t), 0 ≤ x ≤ L, t ∈ [−s, 0], (17) u(0, t) = α1 (t),

u(L, t) = α2 (t),

v(0, t) = β1 (t),

v(L, t) = β2 (t),

0 ≤ t ≤ T,

0 ≤ t ≤ T.

(18) (19)

7

Now, we approximate the integral term of (15) using Lemma 1: Z 1 2J X ∂ α u(x, t) (∆α)4 ∂ α u(x, t) C αl w(α) dα = ∆α γl w(αl )0 Dt u(x, t) − w(α) ∂tα 180 ∂tα 0 l=0

(20)

, α=ξ

where ξ ∈ [0, 1].

Next, we consider the equations (15) − (16) at the point (xj , tk−1+σ ), and using L2 − 1σ formula (definition 4) and Lemma 2, then we have the following multi-term fractional differential equation

∆α

2J X

αl k−1+σ γl w(αl )C 0 Dt uj

l=0

(∆α)4 ∂2v ∂αu (x , t ) − (x , t ) + w(α) j k−1+σ j k−1+σ ∂x2 180 ∂tα

+ f (xj , tk−1+σ , u(xj , tk−1+σ ), u(xj , tk−1+σ−n )) +

2J X

Rkl ,

l=0

Using Lemma 2, we can write where

|Rkl |

3 ∂ u τ 3−αl ≤ (xj , ξ1 ) , 3σ αl Γ(1 − αl )∆α ∂t3

1 ≤ j ≤ M − 1,

ξ1 ∈ [0, tk ],

= u(xj , tk−1+σ ) α=ξ

1 ≤ k ≤ N − 1.

(21)

l = 0, ..., 2J.

Employing Taylor’s expansion for right hand side gives f (xj , tk−1+σ , u(xj , tk−1+σ ), u(xj , tk−1+σ−n ))

=f (xj , tk−1+σ , 2σUjk−1 + (2 − 3σ)Ujk−2 − (1 − σ)Ujk−3 , σUjk−n + (1 − σ)Ujk−n−1 )

+ (u(xj , tk−1+σ ) − 2σUjk−1 − (2 − 3σ)Ujk−2 + (1 − σ)Ujk−3 )fµ (xj , tk−1+σ , ζjk , χkj )

+ (u(xj , tk−1+σ−n ) − σUjk−n − (1 − σ)Ujk−n−1 )fν (xj , tk−1+σ , ζjk , χkj )

= f (xj , tk−1+σ , 2σUjk−1 + (2 − 3σ)Ujk−2 − (1 − σ)Ujk−3 , σUjk−n + (1 − σ)Ujk−n−1 ) + O(τ 2 ).

where ζjk is between u(xj , tk−1+σ ) and 2σUjk−1 u(xj , tk−1+σ−n ) and σUjk−n + (1 − σ)Ujk−n−1 .

+ (2 −

3σ)Ujk−2

− (1 −

σ)Ujk−3 ,

(22)

and χkj is between

Using definition of v(x, t), we can write

∂2u (xj , tk−1+σ ), 1 ≤ j ≤ M − 1, 1 ≤ k ≤ N − 1. (23) ∂x2 Applying the linear operator < (definition 5) on both sides of the equations (21) and (23) 2J X (∆α)4 ∂αu ∂2v C αl k−1+σ ∆α γl w(αl )<(0 Dt uj w(α) α (xj , tk−1+σ ) =
l=0

α=ξ

+
∂2u
2J X

Rkl , (24)

l=0

(25)

∂2u ∂2u ∂2u (xj , tk−1+σ ) = σ< 2 (xj , tk ) + (1 − σ)< 2 (xj , tk−1 ) + O(τ 2 ), 2 ∂x ∂x ∂x h4 ∂ 6 u k (θ , tk−1+σ ) + O(τ 2 ), θjk ∈ (xj−1 , xj+1 ) = σδx2 Ujk + (1 − σ)δx2 Ujk−1 + 240 ∂x6 j h4 ∂ 6 u k = δx2 Ujk−1+σ + (θ , tk ) + O(τ 2 ), θjk ∈ (xj−1 , xj+1 ). (26) 240 ∂x6 j

8

Similarly, we have <

h4 ∂ 6 u k ∂2v (xj , tk−1+σ ) = δx2 Vjk−1+σ + (θ , tk ) + O(τ 2 ), 2 ∂x 240 ∂x6 j

θjk ∈ (xj−1 , xj+1 ).

(27)

Substitute above equations into (24) and (25), we obtain

∆α

2J X l=0

αl k−1+σ γl w(αl )<(C ) + δx2 Vjk−1+σ =
4 α (∆α) ∂ u − (1 − σ)Ujk−3 , σUjk−n + (1 − σ)Ujk−n−1 ) + w(α) α (xj , tk−1+σ ) 180 ∂t

α=ξ

+<

2J X

+

h4 ∂ 6 u k (θ , tk ) 240 ∂x6 j

Rkl ,

(28)

l=0


h4 ∂ 6 u k (θ , tk ) + O(τ 2 ), 240 ∂x6 j

(29)

or above equation (28) and (29) can be written as

∆α

2J X l=0

αl k−1+σ γl w(αl )<(C ) + δx2 Vjk−1+σ =
− (1 − σ)Ujk−3 , σUjk−n + (1 − σ)Ujk−n−1 ) + (R1 )kj ,


=

δx2 Ujk−1+σ

where

+

(R1 )kj

and (R2 )kj =

(R2 )kj ,

1 ≤ j ≤ M − 1,

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

1 ≤ k ≤ N − 1,

(∆α)4 ∂αu = w(α) α (xj , tk−1+σ ) 180 ∂t 4

(30) (31)

2J

+ α=ξ

6

h ∂ u k (θ , tk ) + O(τ 2 ). 240 ∂x6 j

X h4 ∂ 6 u k Rkl , (θj , tk ) + < 6 240 ∂x l=0

(R)kj is the truncation error given below |(R)kj | = |(R1 )kj + (R2 )kj |

≤ |(R1 )kj | + |(R2 )kj |

ˆ 2 + h4 + (∆α)4 ), ≤ C(τ

1 ≤ j ≤ M − 1,

1 ≤ k ≤ N − 1.

(32)

since the solution u(x, t) is smooth in given bounded domain therefore, there exist a positive constant Cˆ independent of τ and h. Taking notice of initial and boundary conditions Ujk = φ(xj , tk ), U0k V0k

= α1 (tk ), = β1 (tk ),

0 ≤ j ≤ M,

k UM k VM

= α2 (tk ),

= β2 (tk ),

−n ≤ k ≤ 0,

(33)

1 ≤ k ≤ N.

(35)

1 ≤ k ≤ N,

(34)

9

Finally, we have our compact difference scheme for equations (1)-(4), by replacing ’Ujk ’ by ukj and ’Vjk ’ by vjk in above equations and omitting the truncation errors (R1 )kj and (R2 )kj , we get as follows ∆α

2J X l=0

αl k−1+σ γl w(αl )<(C ) + δx2 vjk−1+σ =
− (1 − σ)uk−3 , σuk−n + (1 − σ)uk−n−1 ), j j j

1 ≤ j ≤ M − 1,


≤ k ≤ N − 1,

(36)

1 ≤ k ≤ N − 1,

(37) (38) (39) (40)

4. Solvability, stability and convergence In this section, solvability, stability and convergence will be proved for the equations (1)-(4). Given below two Lemmas will be used to achieve the target of the analysis of the scheme. The expression (., .) is used for defining inner product. For any u, v ∈ Vh , we introduce the discrete inner product

(u, v) =

M −1 X

and

uj vj ,

norms

j=1

||u||2 = (u, u),

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

||
max

|uj |.

1≤j≤M −1

Lemma 6. [7] For any grid function u ∈ Vh , it holds that L2 ||δx u||2 , 6 L2 2 2 ||δx u||2 ≤ ||δ u|| . 6 x ||u||2 ≤

(41) (42)

Proof: For the proof of equation (41), reader can follow [7] and equation (42) can be easily derived with the help of (41), discrete Green formula and the Cauchy-Schwartz inequality. Lemma 7. [17] For any grid function u ∈ Vh , it holds that 1 ||u||2 ≤ ||
(43)

Proof: Expression in below given equation for operator < can be followed in [48] and references therein.     h2 2 h2 2 ||


(44)

Noticing that u0 =uM =0 and applying discrete Green formula, we can write (δx2 u, u) = −(δx u, δx u)

(45)

10

Using inverse estimates ||δx2 u||2 ≤

4 2 h2 ||δx u||

and ||δx u||2 ≤

4 2 h2 ||u|| ,

we can write

2 1 ||
(46)

Now, we begin with an analysis of the scheme in terms of stability, convergence, and solvability. Firstly, we will analyze the numerical stability of the compact difference scheme (36)-(40). The numerical stability means that if a small perturbation is introduced in the system then it implies the numerical solution has a small perturbation. For this purpose, we suppose that {zjk , zˆjk |0 ≤ j ≤ M, −n ≤ k ≤ N } be the solution of

∆α

2J X l=0

αl k−1+σ γl w(αl )<(C ) + δx2 zˆjk−1+σ =
− (1 − σ)zjk−3 , σzjk−n + (1 − σ)zjk−n−1 ) + fˆjk ,

<ˆ zjk−1+σ = δx2 zjk−1+σ + gˆjk , zjk = φ(xj , tk ), 0 ≤ j ≤ M, k z0k = α1 (tk ), zM = α2 (tk ),

1 ≤ j ≤ M − 1,

zˆ0k

1 ≤ k ≤ N.

= β1 (tk ),

k zˆM

= β2 (tk ),

−n ≤ k ≤ 0,

1 ≤ j ≤ M − 1,

1 ≤ k ≤ N − 1,

1 ≤ k ≤ N − 1,

(47) (48) (49) (50)

1 ≤ k ≤ N,

(51)

where fˆjk , gˆjk are small perturbations. Let ρkj = zjk − ukj , ρˆkj = zˆjk − vjk , 0 ≤ j ≤ M , −n ≤ k ≤ N . Definition 6. A numerical method for (1)-(4) is stable if the discrete numerical solutions uk , v k satisfying (36)-(40) and z k , zˆk satisfying (47)-(51) are such that solution is bounded in terms of perturbations. Theorem 1. (Stability) The compact difference scheme (36)-(40) is stable with respect to the perturbations fˆjk , gˆjk , i.e. ! α 18 w(1)T ||<ρk ||2 ≤ µ2 ||<ρ0 ||2 + 2 max ||fˆk−1+σ ||2 + max ||ˆ g k−1+σ ||2 . (52) −n≤k≤N 3(1 − α)Γ(2 − α) L4 −n≤k≤N

where µ2 is a constant given in equation (76).

Subtracting equations (36)-(40) from equations (47)-(51) respectively. Then, we obtain the following

∆α

2J X l=0

αl k−1+σ γl w(αl )<(C ) + δx2 ρˆk−1+σ = <ρk−1+σ + fˆjk−1+σ , 0 Dt ρj j j

<ˆ ρk−1+σ j k ρj = 0, ρk0 = 0,

1 ≤ j ≤ M − 1,

1 ≤ k ≤ N − 1, (53)

=

δx2 ρk−1+σ j

0 ≤ j ≤ M, ρkM

= 0,

+

ρˆk0

gˆjk−1+σ ,

1 ≤ j ≤ M − 1,

−n ≤ k ≤ 0,

= 0,

ρˆkM

= 0,

1 ≤ k ≤ N.

1 ≤ k ≤ N − 1,

(54) (55) (56)

Taking inner-product of (53), (54) with <ρk−1+σ and <ˆ ρk−1+σ , respectively, and summing up for j from j j 1 to M − 1, we have

11

(∆α

2J X l=0

αl k−1+σ γl w(αl )<(C ), <ρk−1+σ ) + (δx2 ρˆk−1+σ , <ρk−1+σ ) = (<ρk−1+σ , <ρk−1+σ ) 0 Dt ρ

+ (fˆk−1+σ , <ρk−1+σ ), k−1+σ

(<ˆ ρ

k−1+σ

, <ˆ ρ

)=

(57) (δx2 ρk−1+σ , <ˆ ρk−1+σ )

k−1+σ

+ (ˆ g

Noticing ρ , ρˆ ∈ Vh , then we can write k

k

k−1+σ

, <ˆ ρ

).

h2 2 k−1+σ δ ρ ) 12 x 2 h = −(δx ρˆk−1+σ , δx ρk−1+σ ) + (δx2 ρˆk−1+σ , δx2 ρk−1+σ ) 12 = (δx2 ρk−1+σ , <ˆ ρk−1+σ ).

(δx2 ρˆk−1+σ , <ρk−1+σ ) = (δx2 ρˆk−1+σ , uk−1+σ +

On adding (57) and (58), we get (∆α

2J X l=0

+ (ˆ g k−1+σ , <ˆ ρk−1+σ ).

2J X

= ||<ρ

|| + (fˆk−1+σ , <ρk−1+σ ) + (ˆ g k−1+σ , <ˆ ρk−1+σ ) − (δx2 ρk−1+σ , gˆk−1+σ ).

Using -Cauchy inequality P Q ≤ 2J X l=0

(60) (61)

(62)

(63)

1 1 1 k−1+σ 2 αl k−1+σ γl w(αl )<(C ), <ρk−1+σ ) + ||<ˆ ρk−1+σ ||2 + ||δx2 ρk−1+σ ||2 + ||ˆ g || 0 Dt ρ 2 2 2

l=0 k−1+σ 2

(∆α

(59)

αl k−1+σ γl w(αl )<(C ), <ρk−1+σ ) + ||<ˆ ρk−1+σ ||2 = ||<ρk−1+σ ||2 + (fˆk−1+σ , <ρk−1+σ ) 0 Dt ρ

From (58), we have 1 1 k−1+σ 2 1 ||<ˆ ρk−1+σ ||2 = ||δx2 ρk−1+σ ||2 + ||ˆ g || + (δx2 ρk−1+σ , gˆk−1+σ ). 2 2 2 Using (63) into (62) gives (∆α

(58)

 2 2P

+

1 2 2 Q

(64)

(0 <  ≤ 1), above equation can be written as

1 1 1 k−1+σ 2 αl k−1+σ γl w(αl )<(C ), <ρk−1+σ ) + ||<ˆ ρk−1+σ ||2 + ||δx2 ρk−1+σ ||2 + ||ˆ g || 0 Dt ρ 2 2 2

1 k−1+σ 2 ≤ ||<ρk−1+σ ||2 + (fˆk−1+σ , <ρk−1+σ ) + ||ˆ g || + 2 1 k−1+σ 2 ≤ ||<ρk−1+σ ||2 + (fˆk−1+σ , <ρk−1+σ ) + ||ˆ g || + 2

1 ||<ˆ ρk−1+σ ||2 − (δx2 ρk−1+σ , gˆk−1+σ ) 2 1 1 ||<ˆ ρk−1+σ ||2 + ||ˆ g k−1+σ ||2 + ||δx2 ρk−1+σ ||2 . 2 4 (65)

Rewriting above inequality (∆α

2J X l=0

1 αl k−1+σ γl w(αl )<(C ), <ρk−1+σ ) + ||δx2 ρk−1+σ ||2 ≤||<ρk−1+σ ||2 + 2(fˆk−1+σ , <ρk−1+σ ) 0 Dt ρ 2 + ||ˆ g k−1+σ ||2 .

According to Lemma 6

L4 2 k−1+σ 2 ||δ ρ || . 36 x 18 L4 2(fˆk−1+σ , <ρk−1+σ ) ≤ 4 ||fˆk−1+σ ||2 + ||<ρk−1+σ ||2 , L 18 18 1 ≤ 4 ||fˆk−1+σ ||2 + ||δx2 ρk−1+σ ||2 . L 2 ||<ρk−1+σ ||2 ≤

(66)

(67)

(68)

12

Using above inequalities in (66) (∆α

2J X l=0

αl k−1+σ γl w(αl )<(C ), <ρk−1+σ ) ≤ 0 Dt ρ

18 ˆk−1+σ 2 ||f || + ||ˆ ρk−1+σ ||2 + ||ˆ g k−1+σ ||2 . L4

(69)

The L.H.S. of above equation can be bounded from below by ∆α

2J X l=0

αl k−1+σ γl w(αl )
(k)

#

−
!

2J X γl w(αl )τ −αl

"

(k) c0 ||<ρk ||2

!

k−1 X

=

∆α

2J X γl w(αl )τ −αl l=0

(k) <(ck−j−1

Γ(2 − αl )

(k) ck−j )||<ρk ||2

"

(k)


k−1 X j=1


− − − Γ(2 − αl ) j=1 " # k−1 X γ2J w(α2J )τ −α2J (k) (k) (k) (k) k 2 k 2 k 2 c0 ||<ρ || − <(ck−j−1 − ck−j )||<ρ || −
(k)

(k)

<(ck−j−1 − ck−j )ρj

#

(70)

or " # k−1 X w(1)τ −α (k) 18 (k) (k) (k) c ||<ρk ||2 − <(ck−j−1 − ck−j )||<ρk ||2 −
ck0 ||<ρk ||2 ≤

k−1 X j=1

(k)

(k)

(k)

(ck−j−1 − ck−j )||<ρk ||2 + ck−1 ||<ρk ||2 +µ1

18 ˆk−1+σ 2 ||f || + ||ˆ ρk−1+σ ||2 L4 !

+ ||ˆ g k−1+σ ||2 ,

(71)

(72)

where µ1 =

w(1)τ −α . 3Γ(2 − α)

From Lemma 3, we have (k) ck−1

 −α  −α 1−α α 1−α α > k−1− > k− , 1 ≤ k ≤ N, 2 2 2 2

Therefore,

 −α w(1)τ −α w(1)T α N −α w(1)T α α µ1 = = < k− 3Γ(2 − α) 3Γ(2 − α) 3Γ(2 − α) 2 α w(1)T (k) <2 c . 3(1 − α)Γ(2 − α) k−1

(73)

(74)

13

Using (74) into (72) gives ck0 ||<ρk ||2 <

k−1 X j=1

 (k) (k) (k) (ck−j−1 − ck−j )||<ρk ||2 + ck−1 ||<ρk ||2 + 2 ! || ,

w(1)T α 3(1 − α)Γ(2 − α)

18 ˆk−1+σ 2 ||f || + ||ˆ ρk−1+σ ||2 L4

k−1+σ 2

+ ||ˆ g ≤

k−1 X j=1

 (k) (k) (k) (ck−j−1 − ck−j )||<ρk ||2 + ck−1 µ2 ||<ρ0 ||2 + 2 ! || .

w(1)T α 3(1 − α)Γ(2 − α)

18 ˆk−1+σ 2 ||f || L4

k−1+σ 2

+ ||ˆ g

where

(75)

µ2 = 1 + 2

w(1)T α . 3(1 − α)Γ(2 − α)

D = µ2 ||<ρ0 ||2 + 2

Let

(76)

w(1)T α 3(1 − α)Γ(2 − α)

!

18 ˆk−1+σ 2 ||f || + ||ˆ g k−1+σ ||2 . L4

Below given inequality can be proved easily by induction: ||<ρk ||2 ≤ D, ck0 ||<ρk ||2

≤

1 ≤ k ≤ N,

ck0 D.

It holds obviously when k = 0, and assume that the conclusion is valid for k = 1, 2, ..., N − 1, i.e., Then for 1 ≤ k ≤ N , we have

||<ρk ||2 ≤ D,

ck0 ||<ρk ||2 ≤ ≤

1 ≤ k ≤ N − 1.

k−1 X j=1

(ckk−j−1 − ckk−j )||<ρk ||2 + ckk−1 D,

k−1 X j=1

(ckk−j−1 − ckk−j )D + ckk−1 D,

= ck0 D.

Which completes the proof.

We now analyze the convergence of scheme. Theorem 2. (Convergence) Consider that Ujk and Vjk with 0 ≤ j ≤ M, −n ≤ k ≤ N are solutions of the constructed difference scheme given in (36)-(40). Subtracting equations (36)-(40) from equations (1)-(4) respectively. Then, we obtain the following error difference scheme given below ∆α

2J X l=0

 αl k−1+σ γl w(αl )<(C ) + δx2 eˆk−1+σ =
+ (2 − 3σ)Ujk−2 − (1 − σ)Ujk−3 , σUjk−n + (1 − σ)Ujk−n−1 ) − f (xj , tk−1+σ , 2σuk−1 j  k−2 k−3 k−n k−n−1 + (2 − 3σ)uj − (1 − σ)uj , σuj + (1 − σ)uj )

<ˆ ek−1+σ j ekj = 0, ek0 = 0, where,

=

δx2 ejk−1+σ ,

0 ≤ j ≤ M,

ekM = 0, ekj = Ujk −

eˆk0 ukj

1 ≤ j ≤ M − 1,

−n ≤ k ≤ 0,

= 0,

and

eˆkM eˆkj

= 0,

(77) (78) (79)

1 ≤ k ≤ N,

= Vjk − vjk .

(80)

14

Here, (U, V ) and (u, v) are numerical and exact solutions respectively. Using stability theorem we arrive at the following ˆ 3−α + (∆α)4 + h4 ). ||ek ||∞ ≤ C(τ

(81)

where Cˆ is a positive constant. Theorem 3. (Solvability) The established compact difference scheme (36)-(40) is uniquely solvable. ∆α

2J X l=0

αl k−1+σ γl w(αl )<(C ) + δx2 vjk−1+σ =
− (1 − σ)uk−1 , σuk−n + (1 − σ)uk−n−1 ), j j j


=

δx2 uk−1+σ , j

0 ≤ j ≤ M, ukM

= 0,

v0k

1 ≤ j ≤ M − 1,

1 ≤ j ≤ M − 1,

−n ≤ k ≤ 0,

= 0,

k vM

= 0,

1 ≤ k ≤ N − 1,

(82) (83) (84)

1 ≤ k ≤ N.

(85)

Proof: We will use mathematical induction to prove this theorem. It is clear that for −n ≤ k ≤ 0, uk is determined according to the initial condition (38). Now suppose that ul has to be determined. Let k = l in (36) and (37), then we can get the system of linear algebraic equations. The coefficient matrix of this system is strictly diagonally dominant, thus there exists a unique solution for above given system. By inductive principle, the proof ends. 5. Remark Remark 1. The proposed compact difference scheme in our paper can be extended to the multi-delay fractional partial differential equations as follows Z 1 ∂ 4 u(x, t) ∂ α u(x, t) dα + β = u(x, t) + f (x, t, u(x, t), u(x, t − s1 ), u(x, t − s2 ), ..., u(x, t − sq )), w(α) ∂tα ∂x4 0 (x, t) ∈ (0, L) × (0, T ), (86) u(x, t) = φ(x, t), u(0, t) = α1 (t),

(x, t) ∈ [0, L] × [−s, 0],

(87)

∂ 2 u(L, t) = β2 (t), ∂x2

(89)

u(L, t) = α2 (t),

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

t ∈ [0, T ], t ∈ [0, T ].

(88)

Here, si > 0 and i = 1, 2, ..., q and s = max si . When we consider the case of constant coefficient β to be 1≤i≤q

zero, then, we obtain a simplified distributed order fractional differential equation of multiple delays. For the case of constant β6= 0, we can use the estimates of stability and convergence used for single time-delay. Specifically, for β = 1 we get the case discussed in this paper. For obtaining the high convergence order numerical technique of (86) − (89), we consider the constrained time step size as follows τ = s1 /M1 = s2 /M2 = ... == sq /Mq . In order to ensure that Mi , i = 1, 2, ..., q are all integers we can take the largest common time step of si . ˆ c1 and c2 used in numerical scheme analysis part are generic constants. Remark 2. Constants such as C, 6. Numerical experiments In this section, we give a few examples to testify our theoretical results and to check the convergence behavior of the numerical solution in fractional order direction, spatial direction, and temporal direction. We define the error between exact solution and numerical solution {ukj |0 ≤ j ≤ M, −n ≤, k ≤ N } in L2

15

and L∞ norm respectively as follows v u M u X EL2 (h, τ, ∆α) = max th |u(xj , tk ) − ukj |2 , 0≤j≤M

j=1

E∞ (h, τ, ∆α) =

max

0≤j≤M,−n≤k≤N

|u(xj , tk ) − ukj |.

The convergence orders in L∞ and L2 norm are computed by ! ! EL2 (h, 2τ, ∆α) EL2 (2h, τ, ∆α) Order(τ ) = log2 , Order(h) = log2 , EL2 (h, τ, ∆α) EL2 (h, τ, ∆α) ! EL2 (h, τ, 2∆α) Order(∆α) = log2 , EL2 (h, τ, ∆α) ! ! EL∞ (h, 2τ, ∆α) EL∞ (2h, τ, ∆α) Order(τ ) = log2 , Order(h) = log2 , EL∞ (h, τ, ∆α) EL∞ (h, τ, ∆α) ! EL∞ (h, τ, 2∆α) . Order(∆α) = log2 EL∞ (h, τ, ∆α) All numerical experiments are done with the help of MATLAB software. Example 1. First we consider the following problem (1)-(4) with non-zero initial conditions, [30]:  α 9 ∂ u(x, t) ∂ 4 u(x, t) −α dα = + (x2 − x)u(x, t) + f (x, t, u(x, t), u(x, t − s)), (x, t) ∈ (0, 1) × (0, t), α 2 ∂t ∂x4 0 √ 105 π(t − 1)t3/2 ex x(x − 1) − t7/2 ex (x2 + 3x) where f (x, t, u(x, t), u(x, t − s)) = t7/2 x(x + 3) exp(x) + 16 log t Z

1

Γ



+ u2 (x, t − 0.1) − (x2 − x)2 e2x (t − 0.1)7 − (x2 − x)2 t7/2 ex − t7/2 ex (x2 + 7x + 8)

u(x, t) = t7/2 ex (x2 − x),

u(0, t) = 0, 2

u(1, t) = 0,

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

2

x ∈ [0, 1], t ∈ [0, 1],

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

t ∈ [−0.1, 0],

t ∈ [0, 1].

It is easy to verify that exact solution of above problem is u(x, t) = t7/2 ex (x2 − x). 1 To testify the temporal convergence order, firstly, we fix h = 40 and keep varying τ for various values of α = 0.2, 0.4, 0.6, 0.8. In table 1, error is presented in L2 norm EL2 (τ, h, ∆α) and L∞ norm EL∞ (τ, h, ∆α). It is observed that temporal accuracy of order O(τ 2 ) is achieved conforming the theoretical results given in Theorem 2. 1 and keep varying h for distinct values of α. Next to verify the spatial convergence order, we fix τ = 40 In table 2, error is presented in L2 and L∞ norms and it can be observed that spatial accuracy of order O(h4 ) is achieved. Table 3 shows the fractional distributed-order accuracy of fourth order O(∆α)4 obtained on fixing h = τ = 1/40.

16

Table 1. The computational error and convergence orders in time for example 1. τ

EL2 (h, τ, ∆α)

0.2 1/20 1/50 1/80 0.4 1/20 1/50 1/80 0.6 1/20 1/50 1/80 0.8 1/20 1/50 1/80

7.0315e-004 4.5103e-004 5.3647e-005 7.3138e-004 2.3844e-004 3.1944e-005 4.6659e-004 2.4806e-005 5.5179e-006 1.0027e-004 2.5144e-005 6.8951e-006

using L1 formula [29] Order(τ ) 1.7662 1.7698 1.7810 1.5748 1.5776 1.5849 1.3761 1.3809 1.3856 1.1842 1.1864 1.1889

using L2 − 1σ formula [9] Order(τ ) 1.9768 1.9828 1.9964 1.9836 1.9907 1.9987 1.9844 1.9953 1.9996 1.9867 1.9923 2.0023

E∞ (h, τ, ∆α) 5.1258e-004 9.0047e-004 6.4362e-005 7.5833e-004 4.7967e-004 5.8349e-005 2.9711e-004 8.4763e-005 3.1544e-006 6.0016e-004 2.2358e-005 1.9915e-006

using L1 formula [29] Order(τ ) 1.7671 1.7701 1.7823 1.5754 1.5801 1.5851 1.3775 1.3812 1.3867 1.1871 1.1873 1.1891

using L2 − 1σ formula [9] Order(τ ) 1.9844 1.9901 1.9965 1.9817 1.9932 1.9984 1.9930 1.9967 2.0025 1.9971 1.9991 1.9999

Table 2. The computational error and convergence orders in space for example 1. α 0.2 0.4 0.6 0.8

h EL2 (h, τ, ∆α) Order(h) EL∞ (h, τ, ∆α) Order(h) 1/10 5.5112e-006 3.9812 1.7411e-006 4.0041 1/20 6.6244e-007 3.9873 7.6101e-007 3.9917 1/40 3.1088e-008 3.9937 4.2564e-008 3.9905 1/10 8.6158e-006 3.9830 3.4816e-006 4.0136 1/20 6.7366e-007 3.9888 1.9481e-007 3.9914 1/40 4.9421e-008 3.9947 8.3952e-008 3.9901 1/10 9.6871e-006 3.9832 6.1877e-006 4.0033 1/20 6.3962e-007 3.9901 9.6255e-007 3.9911 1/40 7.2764e-008 3.9964 2.3607e-008 3.9878 1/10 3.6833e-006 3.9851 7.9901e-006 4.0007 1/20 5.9105e-007 3.9944 5.2401e-007 3.9910 1/40 1.4687e-008 3.9985 9.1966e-008 3.9825

Table 3. The computational error and convergence orders for fractional distributedorder for example 1. ∆α EL2 (h, τ, ∆α) Order(∆α) 1/2 2.1341e-002 3.9852 1/4 5.0184e-003 3.9876 1/8 1.6251e-003 3.9881 1/16 8.3822e-004 3.9893 Figure 1a

Figure 1b h=

Exact solution

0

= 1/40 0

Numerical solution

α

-0.2 -0.4 -0.6 1

-0.2 -0.4 -0.6 1

1 0.5

t-axis

0.5 0

0

x-axis

1 0.5

t-axis

0.5 0

0

x-axis

17

Table 4. The computational error and convergence orders in time for example 2. α

τ

EL2 (h, τ, ∆α)

0.2 1/20 1/50 1/80 0.4 1/20 1/50 1/80 0.6 1/20 1/50 1/80 0.8 1/20 1/50 1/80

6.1250e-004 2.4033e-005 3.6842e-006 7.6303e-004 3.7181e-005 1.5700e-006 9.6158e-004 4.4617e-005 2.6324e-006 6.4199e-004 7.2006e-005 9.3587e-006

using L1 formula [29] Order(τ ) 1.7628 1.7711 1.7845 1.5686 1.5743 1.5818 1.3779 1.3784 1.3805 1.1889 1.1935 1.1967

using L2 − 1σ formula [9] Order(τ ) 1.9865 1.9982 1.9957 1.9947 1.9969 1.9972 1.9983 1.9996 2.0001 1.9987 1.9991 2.0003

E∞ (h, τ, ∆α) 2.4711e-003 5.3647e-004 1.6438e-005 3.2891e-003 6.0361e-004 2.4782e-005 7.5186e-004 4.4309e-004 9.8217e-005 1.6372e-004 5.6918e-004 2.3001e-005

using L1 formula [29] Order(τ ) 1.7740 1.7785 1.7802 1.5700 1.5768 1.5819 1.3779 1.3818 1.3811 1.1901 1.1948 1.1971

using L2 − 1σ formula [9] Order(τ ) 1.9904 1.9931 1.9965 1.9817 1.9932 1.9984 1.9930 1.9967 2.0025 1.9971 1.9994 2.0019

1 Figure 1: Exact and numerical solution surface with α = 0.3, h=τ = 40

For visualization, in Figures 1a and 1b we plot the exact and numerical solution considering τ = h = indicating an excellent closeness between numerical and exact solution. Therefore, we conclude that our scheme is efficient for the fourth-order distributed-fractional differential equations with time-delay. 1 40

Example 2. Consider another example ( [20], [30]) given below with following initial and boundary conditions. Z

1

∂ 4 u(x, t) ∂ α u(x, t) dα = + exp(x) u(x, t) − u(x, t)2 − u(x, t − 0.1) + G(x, t), α ∂t ∂x4 0 (x, t) ∈ (0, π) × (0, 1), Γ(4 − α)

(6t2 − 6t) + sin x((t − 0.1)3 + 2(t − 0.1) + 4) log t t ∈ [−0.1, 0]

G(x, t) = −(t3 + 2t + 4) sin x(exp +1) + sin x u(x, t) = (t3 + 2t + 4) sin x, u(0, t) = 0, 2

u(π, t) = 0,

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

2

x ∈ [0, π],

t ∈ [0, 1],

∂ u(π, t) = 0, t ∈ [0, 1]. ∂x2

The exact solution of the above problem is u(x, t) = (t3 + 2t + 4) sin x. We implemented the scheme (36)-(40) for finding the numerical solution of above problem. The efficiency of the scheme is numerically examined by taking sufficiently small spatial and temporal steps. 1 Here, we fixed h = 40 and keep varying τ for various values of α = 0.2, 0.4, 0.6, 0.8 like we did in Example 1. In table 4, error is presented in L2 and L∞ norms and confirms the temporal accuracy of order O(τ 2 ) is achieved validating the theoretical results in Theorem 2. 1 Again, in table 5 we fixed τ = 40 and keep varying h for different values of α to testify the spatial convergence and Table 6 validates the fractional distributed order convergence of order O(∆α)4 . 1 Figures 2a and 2b depicts the plot for exact and numerical solution considering τ = h = 40 verifying an excellent closeness between numerical and exact solution.

18

Table 5. The computational error and convergence orders in space for example 2. α 0.4 0.6 0.8

h EL2 (h, τ, ∆α) Order(h) EL∞ (h, τ, ∆α) Order(h) 1/30 4.8460e-005 3.9845 3.8001e-005 3.9879 1/60 6.7658e-007 3.9915 5.5299e-007 3.9923 1/80 3.3255e-008 3.9979 9.4630e-008 3.9979 1/30 1.8651e-005 3.9858 1.4688e-005 3.9941 1/60 9.7650e-007 3.9975 4.0273e-007 3.9980 1/80 8.8850e-008 3.9994 7.9611e-008 4.0010 1/30 2.9273e-005 3.9875 9.6370e-005 3.9980 1/60 6.7965e-007 3.9941 4.2988e-007 3.9991 1/80 4.4620e-008 4.0001 6.3685e-008 4.0011

Table 6. The computational error and convergence orders for fractional distributedorder for example 2. ∆α EL2 (h, τ, ∆α) Order(∆α) 1/2 6.7012e-002 3.9860 1/4 3.5304e-003 3.9873 1/8 9.3830e-003 3.9885 1/16 2.8296e-004 3.9892 Figure 2a

Figure 2b h=

= 1/40

Numerical solution

Exact solution

6 4 2 0 1

6 4 2 0 1

1 0.5

t-axis

0.5 0

0

x-axis

1 0.5

t-axis

0.5 0

0

x-axis

1 Figure 2: Exact and numerical solution surface with α = 0.3, h=τ = 40

7. Conclusions Here we proposed a new efficient linearized numerical scheme for the solution of fourth-order distributed fractional sub-diffusion equation with time-delay on a bounded domain. We approximated the integral using Simpson’s 1/3rd rule which converts distributed fractional derivative term to a multi fractional derivative. After that, the multi-term fractional differential equation is solved numerically with the proposed difference scheme which is better than the schemes available in the literature in terms of the order of convergence for both time and spatial dimensions. Our scheme gives second-order convergence for temporal dimension and fourth-order convergence for spatial dimensions. Analysis of the scheme is done using a discrete energy method. We also demonstrated the same with the help of a few examples given in the numerical experimentation section. 8. Acknowledgement We thank the Ministry of Human Resource Development, Government of India for its financial support.

19

References [1] R.L. Bagley, P.J. Torvik, On the existence of the order domain and the solution of distributed order equations Part I, International Journal of Applied Mathematics, vol. 2, (2000), pp. 865-882. [2] R.L. Bagley, P.J. Torvik, On the existence of the order domain and the solution of distributed order equations Part II, International Journal of Applied Mathematics, vol. 2, (2000), pp. 965-987. [3] M. Caputo, Linear models of dissipation whose Q is almost frequency independent II, Geophysical Journal International, vol. 13, (1967), pp. 529-539. [4] M. Caputo, Elasticitá e Dissipazione (Elasticity and anelastic dissipation), Zanichelli, Bologna, (1969). [5] M. Caputo, Distributed order differential equations modelling dielectric induction and diffusion, Fractional Calculus and Applied Analysis, vol. 4, (2001), pp. 421-442. [6] M. Caputo, Diffusion with space memory modelled with distributed order space fractional differential equations, Annals of Geophysics, vol. 46, (2003), pp. 223-234. [7] A. A. Samarskii, V. B. Andreev, Difference methods for elliptic equations(Russian book). Moscow, Nauka (1976). [8] T.T. Hartley, C.F. Lorenzo, Fractional system identification: An approach using continuous order-distributions, (1999). [9] A. A. Alikhanov, A new difference scheme for the time fractional diffusion equation, Journal of Computational Physics, vol. 280 (2015), pp. 424-438. [10] V.G. Pimenov, A.S. Hendy, R.H. De Staelen, On a class of non-linear delay distributed order fractional diffusion equations, Journal of Computational and Applied Mathematics, vol. 318, pp. 433-443, (2017). [11] T.M. Atanackovic, L. Oparnica, S. Pilipovic, On a nonlinear distributed order fractional differential equation, Journal of Mathematical Analysis and Applications, vol. 328, (2007), pp. 590-608. [12] K. Diethelm, N.J. Ford, Numerical analysis for distributed-order differential equations, Journal of Computational and Applied Mathematics 225 (2009) 96-104. [13] H.S. Najafi, A. R. Sheikhani, A. Ansari, Stability Analysis of Distributed Order Fractional Differential Equations, Abstract and Applied Analysis, vol. 2011, (2011), 12 pages. [14] J.T. Katsikadelis, The fractional distributed order oscillator: a numerical solution, Journal of the Serbian Society for Computational Mechanics, vol. 6, (2012) pp. 148-159. [15] Z. Sun, Z. Zhang, A linearized compact difference scheme for a class of nonlinear delay partial differential equations, Applied Mathematical Modelling, vol. 37, (2013), pp. 742-752. [16] J. T.Katsikadelis, Numerical solution of distributed order fractional differential equations, Journal of Computational Physics vol. 259, (2014), pp. 11-22. [17] P. Zhang, H. Pu, A second-order compact differnce scheme for the fourth order fractional sub-diffusion equation, Numerical Algorithms, vol. 76, (2017), pp. 573-598. [18] Z. Hao, Z. Sun, W. Cao, A fourth-order approximation of fractional derivatives with its applications, Journal of Computational Physics, vol. 281, (2015), pp. 787-805. [19] D. Li, C. Zhang, J. Wen, A note on compact finite difference method for reaction-diffusion equations with delay, Applied Mathematical Modelling, vol. 39, (2015), pp. 1749-1754. [20] H. Ye, F. Liu, V. Anh, Compact difference scheme for distributed-order time-fractional diffusion-wave equation on bounded domains, Journal of Computational Physics, vol. 298, (2015), pp. 652-660. [21] G.H. Gao, H. Sun, Z. Sun, Some high-order difference schemes for the distributed-order differential equations, Journal of Computational Physics, vol. 298, (2015), pp. 337-359. [22] D. Deng, The study of a fourth-order multistep ADI method applied to nonlinear delay reaction-diffusion equations, Applied Numerical Mathematics, vol. 96, (2015), pp. 118-133. [23] G.H. Gao, Z. Sun, Two alternating direction implicit difference schemes with the extrapolation method for the twodimensional distributed-order differential equations, Computers and Mathematics with Applications, vol. 69, (2015), pp. 926-948. [24] Z. Hao, K. Fan, W. Cao, Z. Sun, A finite difference scheme for semilinear space-fractional diffusion equations with time delay, Applied Mathematics and Computation, vol. 275, (2016), pp. 238-254. [25] X.Y. Li, B.Y. Wu, A numerical method for solving distributed order diffusion equations, Applied Mathematics Letters, vol. 53, (2016), pp. 92-99. [26] H. Chen, S. Lu, W. Chen, Finite difference/spectral approximations for the distributed order time fractional reactiondiffusion equation on an unbounded domain, Journal of Computational Physics, vol. 315, (2016), pp. 84-97. [27] Z. Hao, W. Cao, G. Lin, A second-order difference scheme for the time fractional substantial diffusion equation, Journal of Computational and Applied Mathematics, vol. 313, (2017), pp. 54-69. [28] M.L. Morgadoa, M. Rebelo, L.L. Ferras, N. J. Ford, Numerical solution for diffusion equations with distributed order in time using a Chebyshev collocation method, Applied Numerical Mathematics, vol. 114, (2017), pp. 108-123. [29] X. Hua, L. Zhang, A new implicit compact difference scheme for the fourth-order fractional diffusion-wave system, International Journal of Computer Mathematics, vol. 91, (2014), pp. 2215-2231. [30] A. S. Hendy, R. H. De Staelen, V. G. Pimenov, A semi-linear delayed diffusion-wave system with distributed order in time, Numerical Algorithms, vol. 77, (2018), pp. 885-903. [31] S. Nandal, D. N. Pandey, Numerical Solution of Time Fractional Non-linear Neutral Delay Differential Equations of Fourth-Order, Malaya Journal of Mathematik, vol. 7(3), (2019), pp. 579-589.

20

[32] S. Nandal, D. N. Pandey, Numerical solution of non-linear fourth order fractional sub-diffusion wave equation with time delay, Applied Mathematics and Computation, vol. 369, (2019), https://doi.org/10.1016/j.amc.2019.124900. [33] N. Liu, Y. Liu, H. Li, J. Wang, Time second-order finite difference/finite element algorithm for nonlinear time-fractional diffusion problem with fourth-order derivative term, Computers and Mathematics with Applications, vol. 75, (2018), pp. 3521-3536. [34] V. Kiryakova, Generalised Fractional Calculus and Applications, (1993). [35] I. Podlubny, Fractional Differential Equations, (1999), San Diego, Academic Press. [36] K.S. Miller and B. Ross, An introduction to the fractional calculus and fractional differential equations, (1993). [37] K. Diethelm, The analysis of fractional differential equations: An application oriented exposition using differential operators of Caputo type, Springer Science & Buisness Media, (2010). [38] O.P. Agrawal, A general solution for the fourth-order fractional diffusion-wave equation, Fractional Calculus and Applied Anaysis, vol. 3, (2000), pp. 1-12. vol. 3, (2000), pp. 1-12. [39] O. P. Agrawal, A general solution for a fourth-order fractional diffusion-wave equation defined in a bounded domain, Computers & Structures, vol. 79, (2001), pp. 1497-1501. [40] K.B. Oldham, J. Spainer, The fractional calculus, vol. 111 of Mathematics in science and engineering, Academic Press, New York, (1974). [41] R.M. Ganji, H. Jafari, A Numerical Approach for Multi-variable Orders Differential Equations Using Jacobi Polynomials, International Journal of Applied and Computational Mathematics, vol. 5, (2019), p.34. [42] N.A. Kudryashov, Simplest equation method to look for exact solutions of nonlinear differential equations, Chaos, Solitons and Fractals, vol. 24, (2005), pp. 1217-1231. [43] N. Taghizadeh, M. Mirzazadeh, A. S. Paghaleh, J. Vahidi, Exact solutions of nonlinear evolution equations by using the modified simple equation method, Ain Shams Engineering Journal, vol. 3, (2012), pp. 321-325. [44] T.M., Atanackovic, S. Pilipovic, D. Zorica, Distributed-order fractional wave equation on a finite domain. Stress relaxation in a rod, International Journal of Engineering Science, vol. 49, (2011a), pp. 175-190. [45] N. Su, Mass-time and space-time fractional partial differential equations of water movement in soils: theoretical framework and application to infiltration, Journal of Hydrology, vol. 519, (2014), 1792-1803. [46] N. Su, P.N. Nelson, S. Connor, The distributed-order fractional diffusion-wave equation of groundwater flow: theory and application to pumping and slug tests, Journal of Hydrology, vol. 529, (2015), pp. 1262-1273. [47] Z.Z. Sun, The Method of Order Reduction and Its Application to the Numerical Solutions of Partial Differential Equations, Science Press, Beijing, (2009). [48] H.L. Liao, Z.Z. Sun, Maximum norm error bounds of ADI and compact ADI methods for solving parabolic equations, Numerical Methods for Partial Differential Equations: An International Journal, vol. 26, (2010), pp. 37-60. [49] Z. Navickas, T. Telksnys, R. Marcinkevicius, M. Ragulskis, Operator-based approach for the construction of analytical soliton solutions to nonlinear fractional-order differential equations, Chaos, Solitons and Fractals, vol. 104, (2017), pp. 625-634. [50] S.Y. Reutskiy, A new semi-analytical collocation method for solving multi-term fractional partial differential equations with time variable coefficients, Applied Mathematical Modelling, vol. 45, (2017), pp. 238-254. [51] A. Irshad, S. T. Mohyud-Din, N. Ahmed, U. Khan, A New Modification in Simple Equation Method and its applications on nonlinear equations of physical nature, Results in Physics, vol. 7, (2017), pp. 4232-4240. [52] T. A. Nofal, Simple equation method for nonlinear partial differential equations and its applications, Journal of the Egyptian Mathematical Society, vol. 24, (2016), pp. 204-209. [53] Z. Navickas, L. Bikulciene, M. Rahula, M. Ragulskis, Algebraic operator method for the construction of solitary solutions to nonlinear differential equations, Communications in Nonlinear Science and Numerical Simulation, vol. 18, (2013), pp. 1374-1389. [54] H. Su, C. Guo, The global solution of anisotropic fourth-order Schrödinger equation, Advances in Difference Equations, vol. 2019, (2019), p. 2019. [55] S. R. Grace, J. Dˇ z urina, I. Jadlovsk´ a, T. Li, On the oscillation of fourth-order delay differential equations, Advances in Difference Equations, vol. 2019, (2019), p. 118. [56] L. He, J. Lv, Efficient finite element numerical solution of the variable coefficient fractional subdiffusion equation, Advances in Difference Equations, vol. 2019, (2019), pp. 116. [57] A. Singh, S. Das, S. H. Ong, H. Jafari, Numerical Solution of Nonlinear Reaction-Advection-Diffusion Equation, Journal of Computational and Nonlinear Dynamics, vol. 14(4), (2019), p. 041003. [58] S. Huang, B. Wang, Stabilization Conditions for a Class of Fractional-Order Nonlinear Systems, Journal of Computational and Nonlinear Dynamics, vol. 14(5), (2019), p. 054501. [59] H. Hassani, Z. Avazzadeh, A.T. Machado, Solving Two-Dimensional Variable-Order Fractional Optimal Control Problems With Transcendental Bernstein Series, Journal of Computational and Nonlinear Dynamics, vol. 14(6), (2019), p. 061001. [60] T. Abdeljawad, J. Alzabut, On Riemann-Liouville fractional q-difference equations and their application to retarded logistic type model, Mathematical Methods in the Applied Sciences, vol. 41, (2018), pp. 8953-8962. [61] Z.Z. Sun, An unconditionally stable and O(τ 2 +h4 ) order L∞ convergent difference scheme for linear parabolic equations with variable coefficients, Numerical Methods for Partial Differential Equations: An International Journal, vol. 17(6), (2001), pp. 619-631.

Author contributions Use this form to specify the contribution of each author of your manuscript. A distinction is made between five types of contributions: Conceived and designed the analysis; Collected the data; Contributed data or analysis tools; Performed the analysis; Wrote the paper. For each author of your manuscript, please indicate the types of contributions the author has made. An author may have made more than one type of contribution. Optionally, for each contribution type, you may specify the contribution of an author in more detail by providing a one-sentence statement in which the contribution is summarized. In the case of an author who contributed to performing the analysis, the author’s contribution for instance could be specified in more detail as ‘Performed the computer simulations’, ‘Performed the statistical analysis’, or ‘Performed the text mining analysis’. If an author has made a contribution that is not covered by the five pre-defined contribution types, then please choose ‘Other contribution’ and provide a one-sentence statement summarizing the author’s contribution.

Manuscript title: Numerical treatment of non-linear fourth-order distributed fractional sub-diffusion equation with time-delay

Author 1: SARITA NANDAL ☒

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☒

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☒

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☒

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☒

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

Author 2: Dwijendra N Pandey ☐

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☒

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☒

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

Author 3: Enter author name ☐

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☐

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☐

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

Author 4: Enter author name ☐

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☐

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☐

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

Author 5: Enter author name ☐

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☐

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☐

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

Author 6: Enter author name ☐

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☐

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☐

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

Author 7: Enter author name ☐

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☐

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☐

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

Author 8: Enter author name ☐

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☐

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☐

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

Author 9: Enter author name ☐

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☐

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☐

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

Author 10: Enter author name ☐

Conceived and designed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Collected the data Specify contribution in more detail (optional; no more than one sentence)

☐

Contributed data or analysis tools Specify contribution in more detail (optional; no more than one sentence)

☐

Performed the analysis Specify contribution in more detail (optional; no more than one sentence)

☐

Wrote the paper Specify contribution in more detail (optional; no more than one sentence)

☐

Other contribution Specify contribution in more detail (required; no more than one sentence)

1 Author declaration 1. Conflict of Interest No potential conflict of interest exists. We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome. 2. Funding First author is grateful to MHRD for providing fellowship to carry her research work under grant code OH-31-23-200-428. 3. Intellectual Property We confirm that we have given due consideration to the protection of intellectual property associated with this work and that there are no impediments to publication, including the timing of publication, with respect to intellectual property. In so doing we confirm that we have followed the regulations of our institutions concerning intellectual property.

4. Research Ethics We further confirm that any aspect of the work covered in this manuscript that has involved human patients has been conducted with the ethical approval of all relevant bodies and that such approvals are acknowledged within the manuscript. IRB approval was obtained (required for studies and series of 3 or more cases) Written consent to publish potentially identifying information, such as details or the case and photographs, was obtained from the patient(s) or their legal guardian(s).

5. Authorship We attest that all authors contributed significantly to the creation of this manuscript and we confirm that the manuscript has been read and approved by all named authors. Also, we confirm that the order of authors listed in the manuscript has been approved by all named authors.

6. Contact with the Editorial Office

2 Sarita Nandal is the Corresponding Author declared on the title page of the manuscript. We understand that this Corresponding Author is the sole contact for the Editorial process. She is responsible for communicating with the other authors about progress, submissions of revisions and final approval of proofs. Furthermore, we confirm that the email address shown below is accessible by the corresponding author, is the address to which corresponding author’s account is linked, and has been configured to accept email from the editorial office of Elsevier Journal. [email address: [email protected], [email protected]] We the undersigned agree with all of the above.

Author’s name (First, Last)

Signature

Date

1. SARITA NANDAL

29/09/2019

2. DWIJENDRA PANDEY

29/09/2019