Purposive, qualitative, active vision

Journal of Computational Physics 298 (2015) 652–660

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Compact difference scheme for distributed-order time-fractional diffusion-wave equation on bounded domains H. Ye a , F. Liu b,∗ , V. Anh b a b

Department of Applied Mathematics, Donghua University, Shanghai 201620, PR China Mathematical Sciences, Queensland University of Technology, GPO Box 2434, Brisbane, Qld. 4001, Australia

a r t i c l e

i n f o

Article history: Received 21 March 2014 Received in revised form 22 May 2015 Accepted 7 June 2015 Available online 2 July 2015 Keywords: Distributed-order fractional derivative Diffusion-wave equation Compact difference scheme Stability Convergence

a b s t r a c t In this paper, we derive and analyse a compact difference scheme for a distributed-order time-fractional diffusion-wave equation. This equation is approximated by a multi-term fractional diffusion-wave equation, which is then solved by a compact difference scheme. The unique solvability of the difference solution is discussed. Using the discrete energy method, we prove the compact difference scheme is unconditionally stable and convergent. Finally, numerical results are presented to support our theoretical analysis. © 2015 Elsevier Inc. All rights reserved.

1. Introduction An important application of distributed-order equations is to model ultraslow diffusion where a plume of particles spreads at a logarithmic rate [1–3]. When the order of the fractional derivative is distributed over the unit interval, it is useful for modeling a mixture of delay sources [4]. Also, distributed-order equations may be viewed as consisting of viscoelastic and visco-inertial elements when the order of the fractional derivative varies from zero to two [5,6]. Motivated by these applications, some attention has been paid to the fractional partial differential equations (FPDEs) with distributed order [7–10]. Chechkin et al. [11] proposed diffusionlike equations with time fractional derivatives of the distributed order for the kinetic description of anomalous diffusion and relaxation phenomena and proved the positivity of the solutions of their proposed equations. They demonstrated that retarding subdiffusion and accelerating superdiffusion were governed by distributed-order fractional diffusion equations [12]. The fundamental solutions for the one-dimensional time fractional diffusion equation and multi-dimensional diffusion-wave equation of distributed order were obtained by Mainardi et al. [13, 14] and Atanackovic et al. [15], respectively. Atanackovic et al. [16] also proved the existence of the solution to the Cauchy problem for the time distributed order diffusion equation and calculated it by the use of Fourier and Laplace transformations. Furthermore, they studied waves in a viscoelastic rod of finite length, where viscoelastic material was described by a constitutive equation of fractional distributed-order type (see [17]). Luchko [18] proved the uniqueness and continuous dependence on initial conditions for the generalized time-fractional diffusion equation of distributed order on bounded domains. Meerschaert et al. [4] provided explicit strong solutions and stochastic analogues for distributed-order time-fractional diffusion equations on bounded domains, with Dirichlet boundary conditions.

*

Corresponding author. E-mail address: [email protected] (F. Liu).

http://dx.doi.org/10.1016/j.jcp.2015.06.025 0021-9991/© 2015 Elsevier Inc. All rights reserved.

H. Ye et al. / Journal of Computational Physics 298 (2015) 652–660

653

On the other hand, different numerical methods for solving FPDEs have been proposed [19–22]. Recently, Liu et al. [23] proposed some computationally effective numerical methods for simulating the multi-term time-fractional diffusion-wave equations. There are also some papers discussing numerical methods of the distributed-order equations. For example, Diethelm and Ford [24] developed a numerical scheme for the solution of a distributed-order Fractional ordinary DE and gave a convergence theory for their method. Based on the matrix form representation of discretized fractional operators (see [25]), Podlubny et al. [26] extended the range of applicability of the matrix approach to discretization of distributed-order derivatives and integrals, and to numerical solution of distributed-order differential equations (both ordinary and partial). Katsikadelis [27] presented an efficient numerical method to solve linear and nonlinear distributed-order FODEs. However, published papers on numerical methods of the distributed-order FPDEs are sparse. This motivates us to consider effective numerical methods for distributed-order time-fractional diffusion-wave equations. In this paper, we first approximate the integral term in the distributed-order diffusion-wave equation using numerical approximation. Then the given distributed-order equation can be written as a multi-term time fractional diffusion-wave equation. We derive a compact difference scheme which is uniquely solvable for the multi-term fractional diffusion-wave equation. Using the discrete energy method, we prove the compact difference scheme is unconditionally stable and convergent. Finally, two numerical examples are provided to show the effectiveness of our method. The rest of the paper is organized as follows. In Section 2, a compact difference scheme is derived. Section 3 presents the solvability, stability and convergence for the compact difference scheme. Two examples are given in Section 4 and some conclusions are drawn in Section 5. 2. Compact difference scheme Consider the following distributed-order time-fractional diffusion-wave equations  (α )

Dt

u (x, t ) = K

∂ 2 u (x, t ) + f (x, t ) ∂ x2

(2.1)

in an open bounded domain 0 < x < L , 0 < t < T . Here K > 0, x and t are the space and time variables. The time-fractional  (α ) derivative Dt of distributed order is defined by  (α )

Dt

2 c α 0 D t u (x, t )

u (x, t ) =

 (α )dα

(2.2)

1

with the left-side Caputo fractional derivative c0 D tα defined as (see [28])



c α 0 D t u (x, t )

1

t

(n−α ) 0 (t ∂n u ∂ t n (x, t )

=

n

− τ )n−α −1 ∂∂ τun (x, τ )dτ , n − 1 < α < n, , α = n.

and with a continuous non-negative weight function [1, 2], such that the conditions

(2.3)

 : [1, 2] → R that is not identically equal to zero on the interval

2 0 ≤  (α ),  = 0, α ∈ [1, 2],

 (α )dα = W > 0

(2.4)

1

hold true, where W is a positive constant. In this paper, the initial-boundary conditions

u (x, 0) = φ1 (x),

ut (x, 0) = φ2 (x),

u (0, t ) = ψ1 (t ),

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

0 ≤ x ≤ L, 0≤t ≤ T

(2.5) (2.6)

for Eq. (2.1) are considered. Now, we state our numerical method as follows. Step 1: Discretize the integral term in the distributed-order equation. Let us discretize the interval [1, 2], in which the order α is changing, using the grid 1 = ξ0 < ξ1 < · · · < ξq = 2(q ∈ N ), with the steps ξs not necessarily equidistant. We obtain  (α )

Dt

u (x, t ) ≈

q  s =1

where

  c αs  α d s c0 D t s u (x, t ), 0 D t u (x, t ) ξs = q

 (αs )

s =1

αs ∈ (ξs−1 , ξs ], ds =  (αs ) ξs , ξs = ξs − ξs−1 , s = 1, 2, · · · , q.

(2.7)

654

H. Ye et al. / Journal of Computational Physics 298 (2015) 652–660

For simplicity of the presentation, but without loss of generality, we take ξs =

αs =

use the mid-point quadrature rule for approximating the integral (2.2). Let  (α )

Dt

u (x, t ) =

q 

1 q

ξs−1 +ξs

= σ (q ∈ N ) and ds =  (qαs ) . We can

2

=1+

2s−1 ,s 2q

= 1, 2, · · · , q. Then,

α

d s c0 D t s u (x, t ) + R 1 ,

(2.8)

s =1

where R 1 = O (σ 2 ) (see [29]). Consider the following multi-term fractional diffusion-wave equation q 

ds

s =1

 c αs  ∂ 2 u (x, t ) D u ( x , t ) + R = K + f (x, t ), 1 0 t ∂ x2

(2.9)

with the initial-boundary conditions (2.5)–(2.6). Step 2: Solve the multi-term equation. We assume that we are working on a uniform grid xi = ih, i = 0, 1, · · · , M; Mh = L ; tk = kτ , k = 0, 1, · · · , N ; N τ = T . The domain [0, L ] × [0, T ] is covered by h × τ , where h = {xi |xi = ih, 0 ≤ i ≤ M } and τ = {tk |tk = kτ , 0 ≤ k ≤ N }. Suppose u = {uki |0 ≤ i ≤ M , 0 ≤ k ≤ N } is a grid function on h × τ . Introduce the following notations: k− 12

ui

1

k− 12

= (uki + uki −1 ),

δt u i

2

1

τ

(uki − uki −1 ),

1

δx uki− 1 = (uki − uki−1 ),

δx2 uki = (δx uki+ 1 − δx uki− 1 ).

h

2

1

=

h

2

2

Now we show a compact difference scheme for solving the multi-term equation (2.9) with the initial-boundary conditions (2.5)–(2.6). Define, on h × τ , the following grid functions U ik = u (xi , tk ), f ik = f (xi , tk ), 0 ≤ i ≤ M , 0 ≤ k ≤ N. Suppose u (x, t ) ∈

Cx6,,t3 ([0, L ] × [0, T ]). For 1 < αs < 2, using a fully discrete difference scheme derived by Sun and Wu [30] (see also Theorem 8.2.5 in [31]) and noting the initial condition (2.5), we have

⎡

k− 12

c αs 0 Dt U i

⎤

 τ ⎣ k− 12  αs j− 1 αs = − b k − j −1 − b k − δt U i 2 − bkα−s 1 φ2 (xi )⎦ + R 2s , δt U i j μ¯ s j =1 k −1

(2.10)

where

bk s = (k + 1)2−αs − k2−αs , α



1

| R 2s | ≤

(3 − αs ) s = 1, 2, · · · , q .

2 − αs 12

μ¯ s = τ αs (3 − αs ), 2

+

3 −α s

3 − αs

− (1 + 21−αs ) +

1 12

(2.11)

 max |

0≤t ≤tk

3

∂ u 3 −α s |τ , ∂t3 k− 12

Considering equation (2.9) at the point (xi , tk− 1 ) and using the notation U i q  τ ds s =1

=

μ¯ s K 2

2

⎡

1 ⎣δt U k− 2 i

k −1

 α − b s

αs

− bk − j k − j −1

j =1

[U xx (xi , tk ) + U xx (xi , tk−1 )] +

1



⎤

j− 1 δt U i 2

− bk−1 φ2 (xi )⎦ + αs

q 

(2.12)

= 12 (U ik + U ik−1 ), it is natural to have d s R 2s + R 1

s =1

[ f (xi , tk ) + f (xi , tk−1 )] .

(2.13)

1 ≤ i ≤ M − 1, 0 ≤ k ≤ N , i = 0, M , 0 ≤ k ≤ N .

(2.14)

2

Denote an average operator P as follows:



P uki =

1 12 uki





uki−1 + 10uki + uki+1 ,

Acting the operator P on both sides of (2.13) and noticing that

1 2

k− 12

[P U xx (xi , tk ) + P U xx (xi , tk−1 )] = δx2 U i

we obtain q  τ ds s =1

=

μ¯ s

⎡ P

1 ⎣δt U k− 2 i

k− 1 K δx2 U i 2

k −1

 α − b s j =1

+P

k− 1 fi 2

k − j −1 k− 12

+ Ri

,

− bkα−s j



+ O (h4 ),

(2.15)

⎤ j− 1 δt U i 2

− bkα−s 1 φ2 (xi )⎦

1 ≤ i ≤ M − 1, 1 ≤ k ≤ N ,

(2.16)

H. Ye et al. / Journal of Computational Physics 298 (2015) 652–660

655

where k− 12

1

fi

= ( f ik + f ik−1 ),

k− 1 Ri 2

=−

2

q 

d s P R 2s + O (h4 ) + O (σ 2 ).

s =1

From (2.12), we know that there exists a positive constant C 11 such that

|−

q 

1

d s P R 2s | ≤ C 11 τ 1+ 2 σ

q 

s =1

ds .

s =1

Since q 

ds =

s =1

2

q   (αs )

→

q

s =1

 (α )dα = W > 0, 1

we have q 

d s ≤ C 12 ,

s =1

where C 12 is a positive constant. It follows that k− 12

|Ri

1

| ≤ C 1 (τ 1+ 2 σ + h4 + σ 2 ),

(2.17)

where C 1 is a positive constant. In addition, it follows from the initial and boundary value conditions that

Let

U i0 = φ1 (xi ),

0 ≤ i ≤ M,

U 0k = ψ1 (tk ),

U kM = ψ2 (tk ),

uki

0 ≤ k ≤ N.

(2.19)

be the numerical approximation to u (xi , tk ). We can derive the following compact difference numerical scheme

q  τ ds s =1

μ¯ s

⎡

P

1 ⎣δt uk− 2 i

k −1

 α − b s

k − j −1

j =1 k− 12

= K δx2 u i u 0i uk0

(2.18)

k− 12

+ P fi

,

αs

− bk − j



⎤

j− 1 δt u i 2

− bk−1 φ2 (xi )⎦ αs

1 ≤ i ≤ M − 1, 1 ≤ k ≤ N ,

= φ1 (xi ),

0 ≤ i ≤ M,

= ψ1 (tk ),

ukM = ψ2 (tk ),

(2.20) (2.21)

0 ≤ k ≤ N.

(2.22)

3. Analysis of the compact difference scheme 3.1. Solvability It is clear to see that at each time level, (2.20)–(2.22) is a linear tridiagonal system that need to be solved. Since the coefficient matrix is strictly diagonally dominant, the difference scheme (2.20)–(2.22) has a unique solution. This can be written as the following result. Theorem 3.1. The compact difference scheme (2.20)–(2.22) is uniquely solvable. 3.2. Stability Denote the grid function space on h by Vh = { v | v = ( v 0 , v 1 , · · · , v M −1 , v M ), v 0 = v M = 0}. For any u , v ∈ Vh , we define the discrete inner product

(u , v ) = h

M −1  i =1

ui v i

656

H. Ye et al. / Journal of Computational Physics 298 (2015) 652–660

and denote L 2 norm u  =

δx u , δx v = h

√

(u , u ). The H 1 seminorms | · |1 ,  ·  A and the maximum norm  · ∞ are as follows:

M −1





δx u i + 1

2

i =0

δx u , δx v A = δx u , δx v −

  δx v i + 1 , |u |1 = δx u , δx u , u ∞ = max |u i | 0≤ i ≤ M

2

h2 12

(δx2 u , δx2 v ),

δx u  A =

 δx u , δx u A .

Lemma 3.1. If the grid function { v ki |0 ≤ i ≤ M , 0 ≤ k ≤ N } satisfies v k0 = 0, v kM = 0, 0 ≤ k ≤ N, then

  1 1 1

− δx2 v k− 2 , P δt v k− 2 = δx v k 2A − δx v k−1 2A , 2τ

1 ≤ k ≤ N.

Proof. See Lemma 4.2 in [32]. Theorem 3.2. Let { v ki |0 ≤ i ≤ M , 0 ≤ k ≤ N } be the solution of the following difference system q  τ ds s =1

μ¯ s

⎡

P

1 ⎣δt v k− 2 i

k − j −1

j =1 k− 12

k− 12

= K δx2 v i v 0i v k0

k −1

 α − b s

+ gi

= φ1 (xi ), v kM

= 0,

αs

− bk − j



⎤

j− 1 δt v i 2

− bk−1 φ2 (xi )⎦ αs

1 ≤ i ≤ M − 1, 1 ≤ k ≤ N ,

,

(3.1)

0 ≤ i ≤ M,

= 0,

(3.2)

0 ≤ k ≤ N.

(3.3)

Then it holds that

δx v k 2A ≤ δx φ1 2A + +

1 K

1 1− 2q

τ

max{tk

K

1 1− 2q

max{tk

1

, tk2q }

1

k 

1

, tk2q } q

 (α s ) s=1 q(2−αs ) l=1

1

 g l− 2 2

q 

 (αs ) P φ2 2 , 1 ≤ k ≤ N . q(3 − αs ) s =1

(3.4)

1

Proof. Taking the inner product of (3.1) with P δt v k− 2 , we obtain

 q  τ ds

s =1

μ¯ s

k −1

  



 1 1 αs αs αs j − 12 k− 12 k− 12 P δt v k− 2 , P δt v k− 2 − bk − − b P δ v , P δ v b P φ ( x ), P δ v − t t 2 t j −1 k− j k −1 j =1

= K (δx2 v

k− 12

, P δt v

k− 12

) + (g

k− 12

1

, P δt v k− 2 ). α

α

α

s s Applying Lemma 3.1 and Cauchy–Schwarz inequality, noticing that both bk− and bk− − bk−s j are positive, we have 1 j −1

q  τ ds s =1

μ¯ s

1

P δt v k− 2 2 +

K

2τ

δx v k 2A − δx v k−1 2A



⎡ ⎤ q k −1 

  1 1 τ ds ⎣ αs αs = bk− j −1 − bk− j P δt v j − 2 , P δt v k− 2 ⎦ μ¯ s s =1

j =1

q   1 1 1 τ d s αs

bk−1 P φ2 (x), P δt v k− 2 + ( g k− 2 , P δt v k− 2 ) + ¯s μ s =1 ⎡ ⎤ q k −1



  τ ds  1 1 α α j − 2 k − 2 s ⎣ bk − − bk−s j P δt v 2  + P δt v 2  ⎦ ≤ j −1 ¯s 2μ s =1

j =1

  τ ds α

1 1 1 s bk − P φ2 2 + P δt v k− 2 2 + |( g k− 2 , P δt v k− 2 )|. 1 ¯s 2μ q

+

s =1

It follows that

H. Ye et al. / Journal of Computational Physics 298 (2015) 652–660

q  τ ds s =1

≤

μ¯ s

1

P δt v k− 2 2 + ⎡

δx v k 2A − δx v k−1 2A

τ

k −1

 α ⎣ b s

q  τ ds

μ¯ s

s =1

K

αs

− bk− j P δt v k − j −1

j =1

Denote F 0 = Kτ δx v 0 2A and

Fk =

K

τ

δx v k 2A +

q  τ ds

μ¯ s

s =1



⎤



657

j − 12 2 ⎦



+

q  1 1 τ d s αs b P φ2 2 + 2|( g k− 2 , P δt v k− 2 )|. μ¯ s k−1 s =1

⎛ ⎞ k  1 αs ⎝ bk − P δt v j − 2 2 ⎠ . j j =1

Then,

F k ≤ F k −1 +

q  1 1 τ d s αs b P φ2 2 + 2|( g k− 2 , P δt v k− 2 )|. μ¯ s k−1 s =1

Replacing k by l and summing up for l from 1 to k, we obtain

Fk ≤ F0 +

q k  τ d s  αs bl−1 P φ2 2 μ¯ s s =1

+

k 

l =1

1

q

τ d s αs s =1 μ ¯ s bk−l

l =1

g

 q  k   τ ds α 1 s  + b P δt v l− 2 2 . μ¯ s k−l

l− 12 2

l =1

s =1

Therefore,

δx v k 2A ≤ δx φ1 2A + +

k 1 

K

τ

g τ d s αs s =1 μ ¯ s bk−l

q

l =1

q k 1  τ 2 d s  αs bl−1 P φ2 2 , K μ¯ s s =1

l− 12 2



1 ≤ k ≤ N.

l =1

s Observing that bk− ≥ (2 − αs )(k − l + 1)1−αs ≥ (2 − αs )k1−αs , we have l

α

q q q 1  −1 − 1  τ ds αs  1−αs  (αs )  (αs ) bk−l ≥ tk ≥ min{tk2q , tk 2q } . μ¯ s q(2 − αs ) q(2 − αs ) s =1

Also, since

k

s =1

αs

l=1 bl−1

=k

s =1

2−αs

, it follows that

q q k 1  τ 2 d s  αs 1  2−αs  (αs ) bl−1 P φ2 2 = tk P φ2 2 K μ¯ s K q(3 − αs ) s =1

s =1

l =1

≤

1 K

1 1− 2q

max{tk

1

, tk2q }

q  s =1

 (αs ) q(3 − αs )

P φ2 2 .

Finally, we obtain the inequality (3.4) and the theorem is proved. Using Theorem 3.2, we obtain the following stability statement. Theorem 3.3. The compact difference numerical scheme (2.20)–(2.22) is unconditionally stable to the initial values φ1 (x) and φ2 (x) and the forcing term f . 3.3. Convergence We now consider the convergence of the difference approximation. Noticing that U ik is the exact solution of the system

(2.1), (2.5)–(2.6) and uki is the numerical solution of the compact difference approximation (2.20)–(2.22), we denote the error

eki = U ik − uki ,

0 ≤ i ≤ M,

0 ≤ k ≤ N.

658

H. Ye et al. / Journal of Computational Physics 298 (2015) 652–660

Subtracting (2.20)–(2.22) from (2.16), (2.18)–(2.19), we obtain the error equations q  τ ds s =1

= e 0i ek0

μ¯ s

⎡

P

1 ⎣δt ek− 2 i

k − j −1

j =1

k− 1 K δx2 e i 2

= 0, = 0,

k −1

 α − b s

+

k− 1 Ri 2 ,

− bkα−s j



⎤

j− 1 δt e i 2 ⎦

1 ≤ i ≤ M − 1, 1 ≤ k ≤ N ,

(3.5)

0 ≤ i ≤ M, ekM

= 0,

(3.6)

0 ≤ k ≤ N.

(3.7)

Theorem 3.2 implies the error satisfies

δx ek 2A ≤

1 1− 2q

τ K

max{tk

1

1

, tk2q } q

k 

 (α s ) s=1 q(2−αs ) l=1

1

 R l− 2 2 ,

1 ≤ k ≤ N.

(3.8)

Since k− 12

|Ri

1

| ≤ C 1 (τ 1+ 2 σ + h4 + σ 2 ),

and q 

 (αs ) → q(2 − αs ) s =1

2 1

 (α ) d α = C 2 > 0, (2 − α )

there exists a positive C , such that 1

δx ek 2A ≤ C max{ T 2 , T }(τ 1+ 2 σ + h4 + σ 2 )2 . To proceed further, we need the following lemmas. Lemma 3.2. (See [33].) For any mesh function u ∈ Vh , it holds that

2 3

|u |21 ≤ δx u 2A ≤ |u |21 .

Lemma 3.3. (See [34].) For any mesh function u ∈ Vh , it holds that

√

u ∞ ≤

L

2

| u |1 .

Now, we can state the following result. 6, 3

Theorem 3.4. Suppose that the continuous problem (2.1), (2.5)–(2.6) has a smooth solution u (x, t ) ∈ C x,t ( ), and let uki be the uki

solution of the difference scheme (2.20)–(2.22). Then the solution unconditionally converges to u (xi , tk ) as h, τ and σ tend to zero. Furthermore, there is a positive constant C such that the error satisfies k

e ∞ ≤

max{ T ,

√

T}√

4

1

6LC (τ 1+ 2 σ + h4 + σ 2 ),

0 ≤ k ≤ N.

4. Numerical results In order to illustrate the behavior of our numerical method and demonstrate the effectiveness of our theoretical analysis, two examples are now presented. Example 1. Consider the following distributed-order diffusion-wave equation:

2

ν α c0 D tα u (x, t )dα = K

∂ 2 u (x, t ) , ∂ x2

0 < x < 1, 0 < t ≤ T ,

(4.1)

1

where ν is a positive constant that can be physically interpreted as the relaxation time, K is also a positive constant representing the diffusion coefficient. Here, the initial-boundary conditions

H. Ye et al. / Journal of Computational Physics 298 (2015) 652–660

Fig. 1. The numerical approximation of u (x, t ) for the system (4.1)–(4.3) when

659

ν = 3, 4.5, 6, respectively.

Fig. 2. Exact solutions (lines) and numerical solutions (symbols) at t = 0.3 (triangles), t = 0.6 (stars) and t = 0.9 (squares).

u (x, 0) = x2 (1 − x2 ), u (0, t ) = 0,

ut (x, 0) = 0,

u (1, t ) = 0,

0 ≤ x ≤ 1,

0≤t ≤ T

(4.2) (4.3)

for Eq. (4.1) are considered. Using the numerical method described in Section 2, we obtain the numerical solutions (Fig. 1) of the fractional diffusion equation for K = 1, T = 0.5, ν = 3, 4.5, 6 respectively, with h = 0.02, τ = 0.01, σ = 0.1. Example 2. Consider the following distributed-order time-fractional diffusion-wave equation:

⎧2 ⎪ (4 − α )c0 D tα u (x, t )dα = ⎪ ⎨ 1

∂ 2 u (x,t ) ∂ x2

+ f (x, t ), 0 < x < π, 0 < t ≤ T , ⎪ ⎪ ⎩ u (x, 0) = 4 sin x, ut (x, 0) = 2 sin x, 0 ≤ x ≤ π , u (0, t ) = u (π , t ) = 0, 0 ≤ t ≤ T , where

 3

f (x, t ) = sin x t + 2t + 4 +

6t 2 − 6t log t

(4.4)

 .

The exact solution of the above problem is u (x, t ) = (t 3 + 2t + 4) sin x. A comparison of the exact solution and the numerical solution with h = π /100, τ = 0.0176, σ = 0.1 at t = 0.3 (triangles), t = 0.6 (stars) and t = 0.9 (squares) is shown in Fig. 2. From Fig. 2, it can be seen that the numerical solution is in good agreement with the exact solution.

660

H. Ye et al. / Journal of Computational Physics 298 (2015) 652–660

5. Conclusion In this paper, a compact difference scheme for the distributed-order time-fractional diffusion-wave equations on bounded domains has been described. By approximating the integral term in the distributed-order equation using the mid-point quadrature rule, we obtained a multi-term time fractional diffusion-wave equation. It should be noted that Simpson formula could replace the mid-point quadrature rule to obtain the high integration convergence order. For solving the multi-term equation, we constructed a compact difference scheme which has the advantage of high spatial accuracy and unconditionally stability. We proved that the compact difference scheme was unconditionally stable and convergent using the discrete energy method. Two numerical examples demonstrated the effectiveness of theoretical results. The method presented in this paper was designed for one-dimensional problem of distributed-order fractional differential equation. If the method is applied to a high-dimensional problem, the computational complexity and CPU time are considerable. In our future work, we would like to construct effective and easy-to-use numerical schemes for high-dimensional problems. Acknowledgements We would like to thank the two anonymous referees for providing us with constructive comments and suggestions. This work is supported in part by the National Nature Science Foundation of China (No. 11371087). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]

Y.G. Sinai, The limiting behavior of a one-dimensional random walk in a random medium, Theory Probab. Appl. 27 (1982) 256–268. A.V. Chechkin, J. Klafter, I.M. Sokolov, Fractional Fokker–Planck equation for ultraslow kinetics, Europhys. Lett. 63 (2003) 326. A.N. Kochubei, Distributed order calculus and equations of ultraslow diffusion, J. Math. Anal. Appl. 340 (2008) 252–281. M.M. Meerschaert, E. Nane, P. Vellaisamy, Distributed-order fractional diffusions on bounded domains, J. Math. Anal. Appl. 379 (2011) 216–228. C.F. Lorenzo, T.T. Hartley, Variable order and distributed order fractional operators, Nonlinear Dyn. 29 (2002) 57–98. T.T. Hartley, C.F. Lorenzo, Fractional-order system identification based on continuous order-distributions, Signal Process. 83 (2003) 2287–2300. M. Naber, Distributed order fractional sub-diffusion, Fractals 12 (2004) 23–32. I.M. Sokolov, A.V. Chechkin, J. Klafter, Distributed-order fractional kinetics, Acta Phys. Pol. B 35 (2004) 1323–1341. C.H. Eab, S.C. Lim, Fractional Langevin equations of distributed order, Phys. Rev. 83 (2011) 031136. Z. Jiao, Y. Chen, I. Podlubny, Distributed-Order Dynamic Systems: Stability, Simulation, Applications and Perspectives, Springer, 2012. A.V. Chechkin, R. Gorenflo, I.M. Sokolov, V.Y. Gonchar, Distributed order time fractional diffusion equation, Fract. Calc. Appl. Anal. 6 (2003) 259–280. A.V. Chechkin, R. Gorenflo, I.M. Sokolov, Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations, Phys. Rev. E 66 (2002) 046129. F. Mainardi, G. Pagnini, R. Gorenflo, Some aspects of fractional diffusion equations of single and distributed order, Appl. Math. Comput. 187 (2007) 295–305. F. Mainardi, A. Mura, G. Pagnini, R. Gorenflo, Time-fractional diffusion of distributed order, J. Vib. Control 14 (2008) 1267–1290. T.M. Atanackovic, S. Pilipovic, D. Zorica, Time distributed-order diffusion-wave equation. II. Applications of Laplace and Fourier transformations, Proc. R. Soc. A, Math. Phys. Eng. Sci. 465 (2009) 1893–1917. T.M. Atanackovic, S. Pilipovic, D. Zorica, Existence and calculation of the solution to the time distributed order diffusion equation, Phys. Scr. 2009 (2009) 014012. T.M. Atanackovic, S. Pilipovic, D. Zorica, Distributed-order fractional wave equation on a finite domain. Stress relaxation in a rod, Int. J. Eng. Sci. 49 (2011) 175–190. Y. Luchko, Boundary value problems for the generalized time-fractional diffusion equation of distributed order, Fract. Calc. Appl. Anal. 12 (2009) 409–422. F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker–Planck equation, J. Comput. Appl. Math. 166 (2004) 209–219. F. Liu, P. Zhuang, V. Anh, I. Turner, K. Burrage, Stability and convergence of the difference methods for the space–time fractional advection–diffusion equation, Appl. Math. Comput. 191 (2007) 12–20. F. Liu, P. Zhuang, K. Burrage, Numerical methods and analysis for a class of fractional advection–dispersion models, Comput. Math. Appl. 64 (2012) 2990–3007. 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. F. Liu, M.M. Meerschaert, R.J. McGough, P. Zhuang, Q. Liu, Numerical methods for solving the multi-term time-fractional wave-diffusion equation, Fract. Calc. Appl. Anal. 16 (2013) 9–25. K. Diethelm, N.J. Ford, Numerical analysis for distributed-order differential equations, J. Comput. Appl. Math. 225 (2009) 96–104. I. Podlubny, Matrix approach to discrete fractional calculus, Fract. Calc. Appl. Anal. 3 (2000) 359–386. I. Podlubny, T. Skovranek, B.M.V. Jara, I. Petras, V. Verbitsky, Y. Chen, Matrix approach to discrete fractional calculus III: non-equidistant grids, variable step length and distributed orders, Philos. Trans. R. Soc., Math. Phys. Eng. Sci. 371 (2013). J.T. Katsikadelis, Numerical solution of distributed order fractional differential equations, J. Comput. Phys. 259 (2014) 11–22. I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. J.D. Faires, R.L. Burden, Numerical Methods, Brooks/Cole, Cengage Learning, 2013. Z.-z. Sun, X. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193–209. Z. Sun, The Method of Order Reduction and Its Application to the Numerical Solutions of Partial Differential Equations, Science Press, Beijing, 2009. J. Ren, Z.-z. Sun, Numerical algorithm with high spatial accuracy for the fractional diffusion-wave equation with Neumann boundary conditions, J. Sci. Comput. (2013) 1–28. G. Gao, Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. A. Samarskii, V. Andreev, Difference Methods for Elliptic Equations, Nauka, Moscow, 1976.