New compact difference scheme for solving the fourth-order time fractional sub-diffusion equation of the distributed order

Applied Numerical Mathematics 129 (2018) 58–70

Contents lists available at ScienceDirect

Applied Numerical Mathematics www.elsevier.com/locate/apnum

New compact difference scheme for solving the fourth-order time fractional sub-diffusion equation of the distributed order Maohua Ran a,∗ , Chengjian Zhang b a b

School of Mathematics, Sichuan Normal University, Chengdu 610068, China School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China

a r t i c l e

i n f o

Article history: Received 24 March 2017 Received in revised form 25 December 2017 Accepted 3 March 2018 Available online 8 March 2018 Keywords: Distributed-order derivative Boundary value method Finite difference method Quadrature rule Stability and convergence High accuracy

a b s t r a c t In this paper, a class of new compact difference schemes is presented for solving the fourth-order time fractional sub-diffusion equation of the distributed order. By using an effective numerical quadrature rule based on boundary value method to discretize the integral term in the distributed-order derivative, the original distributed order differential equation is approximated by a multi-term time fractional sub-diffusion equation, which is then solved by a compact difference scheme. It is shown that the suggested compact difference scheme is stable and convergent in L ∞ norm with the convergence order O (τ 2 + h4 + (γ ) p ) when a boundary value method of order p is used, where τ , h and γ are the step sizes in time, space and distributed-order variables, respectively. Numerical results are reported to verify the high order accuracy and efficiency of the suggested scheme. Moreover, in the example, comparisons between some existing methods and the suggested scheme is also provided, showing that our method doesn’t compromise in computational time. © 2018 Published by Elsevier B.V. on behalf of IMACS.

1. Introduction Consider the following initial-boundary value problem

Dtw u (x, t ) +

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

(1.1)

u (x, 0) = 0, 0 ≤ x ≤ L ,

(1.2)

u (0, t ) = u 0 (t ), u ( L , t ) = u L (t ), 0 ≤ t ≤ T ,

(1.3)

∂ 2u ∂ 2u (0, t ) = ηa (t ), ( L , t ) = ηb (t ), 0 ≤ t ≤ T , 2 ∂x ∂ x2

(1.4)

where u 0 , u L , ηa , ηb are given functions which satisfy u 0 (0) = u L (0) = 0. The fractional derivative Dtw of distributed order is defined by

*

Corresponding author. E-mail addresses: [email protected] (M. Ran), [email protected] (C. Zhang).

https://doi.org/10.1016/j.apnum.2018.03.005 0168-9274/© 2018 Published by Elsevier B.V. on behalf of IMACS.

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

1 Dtw u (x, t ) =

γ

w (γ ) 0C D t u (x, t )dγ

59

(1.5)

0

with a non-negative continuous function w satisfying

1 w (γ )dγ = c 0 > 0,

(1.6)

0

and with the Caputo fractional derivative

C γ 0 D t u (x, t )

=

1

(1 − γ )

t 0

C γ 0 Dt

defined as

1 ∂u (x, ν )dν , 0 < γ < 1. (t − ν )γ ∂ ν

(1.7)

It is worth noting that the zero initial value is only considered here. The numerical method proposed in this paper is also applicable for non-zero initial value problem since the latter can be reduced to a zero initial value one by the use of a proper transformation (cf. [4]). Equation (1.1) is employed in physics to simulate ultraslow diffusion with a logarithmic growth. Recently much attention, such as existence and uniqueness of fundamental solution, has been paid on such diffusion-like equations with time distributed-order derivative because their application potential in physics, system control and signal processing, see for example [5,12,13,15,17,19] and references therein. Since the analytical solution of the distributed-order equation is not easy to obtain in general, developing efficient numerical method for the solution of distributed-order problem, especially high-order accuracy numerical algorithm, becomes necessary and important. Nevertheless, relatively few works have been done on this topic. And for all we know, most of the existing numerical work deal with the ordinary differential equation of distributed order, see for example [1,7,8]. Here we will briefly review some of these recent advances focusing mainly on the numerical solution of partial differential equation of distributed order. In [23], Ye, Liu and Anh derived and analyzed a compact difference scheme for the time distributed-order diffusion-wave equation on a bounded domain, where the integral term in the distributed-order derivative is discretized by the mid-point quadrature rule and the Caputo fractional derivatives in the resulting multi-term fractional diffusion equation are approximated by the known L1 formula (cf. [21]). The similar treatment ways can also be founded in [16] and [14], but where the Caputo fractional derivatives are solved by the backward finite difference formula and the reproducing kernel method, respectively. And very recently, Gao et al. [9,10] presented some high order difference schemes for both one-dimensional and two-dimensional fractional diffusion equations of distributed order, first by the composite trapezoid formula or the composite Simpson formula to discretize the distributed integral and then by the weighted and shifted Grünwald formula to treat the time-fractional derivatives in the resultant multi-term fractional diffusion equation; among which the alternating direction and the Richardson extrapolation techniques are also involved. One can observe that the above mentioned numerical methods are all constructed on the basis of the classical numerical quadrature formulae, which have no more than accuracy of order 4 in the distributed-order variable. In contrast, our attention in this paper will be paid on a new class of stable quadrature rules which can achieve accuracy of order arbitrary in theory, and then apply the new quadrature rule to the fourth-order time fractional sub-diffusion equation of the distributed order which is more complex than the time distributed-order sub-diffusion and diffusion-wave equations. The rest of the paper is organized as follows. In Section 2, we shall review the derivation of the new quadrature rule and introduce some lemmas needed in what follows. In Section 3, a class of new compact difference schemes for the problem (1.1)–(1.4) is derived. In Section 4, the stability and convergence of the suggested scheme are proved by utilizing discrete energy method. Finally, a numerical example is given to illustrate the effectiveness and high accuracy of the suggested scheme. 2. New quadrature rule and some lemmas In this paper we shall use a new quadrature rule based on boundary value method to approximate the integral term in the distributed-order derivative (1.5). To describe this idea, we consider the quadrature problem

I  (γ ) = φ(γ ),

γ ∈ [0, γs ], I (0) = φ0 ,

where φ is given smooth function. The solution of the problem (2.1) at

(2.1)

γn is

γn I (γn ) = φ0 +

φ(γ ) dγ , γn ∈ [0, γs ]. 0

(2.2)

60

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

In order to find the numerical approximation I n to I (γn ), we apply an unconventional k-step linear multistep method (LMM) of order p with k1 initial conditions and k2 final conditions to the problem (2.1), where k = k1 + k2 . Then the following discrete problem is obtained, k −k1

αi+k1 I n+i = γ

i =−k1

k −k1

βi +k1 φn+i , n = k1 , . . . , s − k2 ,

(2.3)

i =−k1

where φn = φ(γn ) with γn = n · γ , γ = γs /s. Since only I 0 (= φ0 ) is provided by the continuous problem, the remaining k1 − 1 initial values I 1 , . . . , I k1 −1 and k2 final values I s−k2 +1 , . . . , I s may be obtained by imposing k − 1 additional methods as follows k 

αi( j) I i = γ

k 

i =0 k  i =0

( j)

β i φ i , j = 1, . . . , k 1 − 1,

(2.4)

i =0

αk(−j)i I s−i = γ

k  i =0

( j)

βk−i φs−i , j = s − k2 + 1, . . . , s,

( j)

(2.5)

( j)

where the coefficients {αi } and {βi } are chosen such that the schemes (2.4) and (2.5) have the same order as the main scheme (2.3). The method (2.3)–(2.5) obtained in this way are often called boundary value method (BVM) with (k1 , k2 )-boundary conditions. In this case the given continuous initial value quadrature problem (2.1) is approximated by a discrete boundary value problem. It is evident that the BVM (2.3)–(2.5) contains the traditional LMMs as the particular case of k1 = k and k2 = 0. Let I = ( I 1 , I 2 , . . . , I s ) T , φ = (φ0 , φ1 , . . . , φs ) T , A e = [a A ] ∈ Rs×(s+1) with

⎡

α1(1)

⎢ ⎢ .. ⎢. ⎢ ⎢ (k −1) ⎢α 1 ⎢ 1 ⎢ ⎢ α1 ⎢ ⎢ ⎢ α0 ⎢ A=⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

α2(1)

···

αk(1)

.. .

.. .

.. .

⎤

α2(k1 −1) · · · α2 ··· α1 ···

αk(k1 −1) αk αk−1

αk

..

..

..

.

..

.

α0 (s−k2 +1)

.

α1 (s−k2 +1)

.

..

.

· · · αk−1 (s−k +1)

α0

α1

· · · αk−1 2

.. .

.. .

.. .

α0(s)

α1(s)

· · · αk−1

.. . (s)

(k −1)

(1)

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ s×s ⎥∈ , ⎥ ⎥ ⎥ ⎥ αk ⎥ ⎥ (s−k2 +1) ⎥ αk ⎥ ⎥ ⎥ .. ⎥ . ⎦

αk(s)

and a = [α0 , . . . , α0 1 , α0 , 0, . . . , 0] T ∈ Rs . Similarly, B e can be defined by replacing (2.3)–(2.5) can be rewritten in the matrix form as

α with β in A e . Then the BVM

A I = γ B e φ − I 0 a = γ B e φ − φ 0 a .

(2.6)

If the BVM (2.3)–(2.5) is consistent and A is nonsingular one has

A −1 a = −e ,

(2.7)

where e = (1, 1, . . . , 1) T ∈ Rs . This together with (2.6) leads to

I = γ A − 1 B e φ + φ 0 e .

(2.8)

The above equality implies that

I n = γ

s 

ρnj φ(γ j ) + φ0 , n = 1, 2, . . . , s,

j =0

where ρnj denotes the nth row and ( j + 1)st column of the matrix A −1 B e . With respect to the approximation error between I n and I (γn ), some conclusions have been shown.

(2.9)

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

Table 1 Quadrature weight coefficients

ETR2s-4

5 12 33 98

GBDF-5

0

GAM-5

29 90 147 470 141 446

GAM-3

ETR2s-6 TOM-6

61

ρsj derived from different BVMs with s = 6. 13 12 9 7 1473 470 973 720 669 470 2057 1466

1

1

81 98 − 792 235 511 720 231 470 508 937

54 49 1443 235 433 360 363 235 226 153

5 4 9 7 1161 470 929 720 669 470 2057 1466

11 12 81 98 − 597 235 283 360 231 470 508 937

1 3 33 98 39 235 27 80 147 470 141 446

Lemma 2.1. (cf. [2,6]) If the used BVM (2.3)–(2.5) is consistent of order p, then the quadrature rule (2.9) is convergent of order p. More precisely, it holds that

| I n − I (γn )| ≤ C 1 (γ ) p for n · γ ≤ γs , where C 1 = C 0 || A −1 || with 1-norm || · || and C 0 is a constant independent of n and γ . Remark 2.2. From the quadrature rule (2.9), one can find that the quadrature weight coefficients ρnj depend only on the used BVM. That is, ρnj are independent of the form of the integrand φ(γ ). And also, the approximation order p is only determined by the used BVM. In particular, Brugnano and Trigiante in their monograph [3] have listed five types of BVMs of order up to 10, including the generalized backward differentiation formulae (GBDFs), the generalized Adams methods (GAMs), the extended trapezoidal rules of the first and second classes (ETRs, ETR2s) and the top-order methods (TOMs). Matrices A corresponding to the above BVMs can be all nonsingular for particular choice of k1 and k2 , and the entries of A −1 have a uniform upper bound with respect to s, for s sufficiently large, meanwhile the above BVMs are all A k1 ,k2 -stable which is a generalization of the Absolute stability for initial value methods. Hence the new quadrature rule (2.9) is competitive compared with the classical numerical quadrature formulae such as the mid-point quadrature rule (cf. [23]), the composite trapezoid formula and the composite Simpson formula (cf. [10]). In the following we shall use only the last row of matrix A −1 B e , i.e., ρsj ( j = 0, 1, . . . , s) if set γs = 1. To do this we report six groups of quadrature weight coefficients ρsj ( j = 0, 1, . . . , s) in Table 1, which are derived from the GAM of order 3 and 5 (GAM-3, GAM-5), the ETR2s of order 4 and 6 (ETR2s-4, ETR2s-6), the GBDF of order 5 (GBDF-5) and the TOM of order 6 (TOM-6) with s = 6, respectively. All the parameters related to the BVMs mentioned above can be found in [3]. Note that

1 Dtw u (x, t ) =

γ

w (γ ) 0C D t u (x, t )dγ

(2.10)

0

can be viewed as the special case of the problem (2.2) with φ(γ ) = w (γ ) of order p to (2.10), it follows from Lemma 2.1 that

Dtw u (x, t ) = γ

s 

C γ 0 D t u (x, t ), φ0

γ

ρsj w (γ j ) 0C D t j u (x, t ) + O((γ ) p ).

= 0 and γs = 1. Applying a BVM

(2.11)

j =0

So far, Equation (1.1) is converted to the multi-term fractional differential equation

γ

s 

dj

C γj 0 D t u (x, t ) +

j =0

∂ 4u (x, t ) = f (x, t ) + O((γ ) p ), ∂ x4

(2.12)

where d j = ρsj w (γ j ). Take positive integers M , N. Let h = L / M , τ = T / N be the uniform mesh step sizes in space and time, respectively. Denote xi = ih, i = 0, 1, . . . , M , tn = nτ , n = 0, 1, . . . , N , h = {xi | 0 ≤ i ≤ M }, τ = {tn | 0 ≤ n ≤ N }. Let Vh = { v | v = ( v 0 , v 1 , · · · , v M )} be grid function space on h , and V˚ h = { v | v ∈ Vh , v 0 = v M = 0}. For any v ∈ Vh , we define

1

1

h

h2

δx v i = ( v i − v i −1 ), δx2 v i =

( v i −1 − 2v i + v i +1 ), δx4 v i = δx2 (δx2 v i ).

For u , v ∈ V˚ h , we introduce the discrete inner product

(u , v ) = h

M −1  i =1

ui v i ,

62

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

and norms and seminorms

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

max

1 ≤ i ≤ M −1

|u i |, |u |21 = h

M M −1   (δx u i )2 , |u |22 = h (δx2 u i )2 . i =1

i =1

To achieve second-order accuracy in time we introduce the following lemma which is proposed by Tian et al. [22] and collated by Gao et al. [10]. Lemma 2.3. (cf. [10,22]) Suppose that f (t ) ∈ C 3 (R) (0 < γ < 1), and let

1

γ

−∞ D t f (t ) =

t

d

(1 − γ ) dt

−∞

f (ν )

(t − ν )γ

dν .

Then ∞ 1  (γ )

γ

−∞ D t f (t ) =

λk

τγ

f (t − kτ ) + O (τ 2 )

k =0

uniformly holds in t ∈ R as τ → 0, where

γ  (γ ) (γ )

γ  (γ ) γ (γ ) (γ ) λ0 = 1 + g 0 , λk = 1 + g k − g k −1 , k ≥ 1, 2

and (γ )

g0

(γ )

= 1, g k

2

2

γ + 1 (γ ) = 1− g k −1 , k ≥ 1.

(2.13)

(2.14)

k

(γ )

Lemma 2.4. (cf. [10]) Let {λk } be defined as in (2.13). Then for any positive integer m and real vector ( v 0 , v 1 , . . . , v m ) T ∈ Rm+1 , it holds that

 n  m   (γ ) λk v n−k v n ≥ 0. n =0

k =0

Lemma 2.5. Let

μ = γ

s 

dj ·

j =0

1

τ γj

(γ )

λ0 j ,

(2.15)

then it holds that

μ=

1

O(τ | ln τ |)

, as τ → 0.

Proof. This proof is similar to [23], see also [10]. Since w (γ ) is non-negative continuous function on [0, 1], it follows from the integral mean theorem that

μ = γ

s 

dj ·

j =0

1

1

τ



w (γ ) 1 +

→

(γ j )

γ j λ0

= γ

s 

ρsj w (γ j ) ·

j =0

1

τ γj

1+

γj  2

γ 1 dγ 2 τγ

0

∗



= w (γ ) 1 +

γ∗

1

2

1



γ∗

dγ = w (γ ∗ ) 1 + τγ 2

0

where γ ∗

∈ (0, 1). This inequality implies that μ =

1

1−τ

τ ln τ1

O (τ | ln τ |) as

Also, the following lemmas will be used later on.

,

τ → 0. 2

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

63

Lemma 2.6. (cf. [11]) For v ∈ V˚ h , and δx2 v 0 = δx2 v M = 0. It holds that

(δx4 v , v ) = (δx2 v , δx2 v ) = | v |22 . Lemma 2.7. (cf. [20]) Let v ∈ V˚ h . Then

(δx2 v , v )

=

−| v |21 ,

√

L

|| v || ≤ √ | v |1 , || v ||∞ ≤ 6

L

2

| v |1 , | v |21 ≤

4 h2

|| v ||2 .

Lemma 2.8. (cf. [18]) For v ∈ V˚ h , the following inequality holds,



2 sin(π h/2)

h

| v |1 ≤ | v |2 .

3. New compact difference scheme In the section we are devoted to deriving an higher order finite difference scheme for the initial-boundary value problem (1.1)–(1.4). Let

U in = u (xi , tn ), z(x, t ) =

∂ 4u (x, t ), zni = z(xi , tn ). ∂ x4

(3.1)

Suppose u (x, ·) ∈ C 8 [0, L ]. Using Taylor expansion and the above notations, we have

∂ 4u h2 ∂ 6 u h2 ∂ 2 z ( xi , t n ) + (xi , tn ) + O(h4 ) = zni + (xi , tn ) + O(h4 ) 4 6 6 ∂x 6 ∂ x2 ∂x 

 h2 2 n h2 ∂ 4 z h2 2 n 4 δx z i − (ζ , t ) + O ( h ) = 1 + δ zi + O (h4 ) = zni + n i 6 12 ∂ x4 6 x

δx4 U in =

1

= ( zni−1 + 4zni + zni+1 ) + O(h4 ), 1 ≤ i ≤ M − 1, 6

(3.2)

where ζi ∈ (xi −1 , xi +1 ). For v ∈ V˚ h we introduce the difference operator

⎧ ⎨ 1 (v i −1 + 4v i + v i +1 ), 1 ≤ i ≤ M − 1, Hv i = 6 ⎩ v i , i = 0, M .

Considering (2.12) at the point (xi , tn ), one has

γ

s 

dj

C γj 0 D t u ( xi , t n ) +

j =0

∂ 4u (xi , tn ) = f (xi , tn ) + O((γ ) p ), 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N . ∂ x4

(3.3)

Performing the operator H on both sides of the above equation leads to

γ

s 

dj

j =0

∂ 4u C γj D H u ( x , t ) + H ( xi , t n ) n i 0 t ∂ x4

= H f (xi , tn ) + O((γ ) p ), 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N .

(3.4)

Based on (3.2), it follows that

H

∂ 4u (xi , tn ) = δx4 U in + O(h4 ), 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N . ∂ x4 γ

(3.5) γ

Since the Caputo fractional derivative 0C D t f (t ) is equivalent to the Riemann–Liouville fractional derivative −∞ D t f (t ) with f (t ) = 0 when t ≤ 0, for more details on this treatment see [9,10]. Hence, by the use of Lemma 2.3 yields C γj n 0 Dt U i

n 1  (γ j )

= γj τ

λk

U in−k + O (τ 2 ), 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N .

k =0

Substituting (3.5) and (3.6) into (3.4) yields

(3.6)

64

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

s  dj γ t γ j



j =0

 n  (γ j ) n−k λk HU i + δx4 U in = H f in + R ni ,

(3.7)

k =0

where f in = f (xi , tn ) and

R ni = O (τ 2 )γ

s 

d j + O (h4 + (γ ) p ), 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N .

(3.8)

j =0

In view of

γ

s 

d j = γ

s 

j =0

1

ρsj w (γ j ) →

j =0

w (γ )dγ = c 0 > 0,

(3.9)

0

so there exists a positive constant c such that (see also [23])

γ

s 

d j ≤ c.

(3.10)

j =0

This together with (3.8) leads to

| R ni | ≤ C (τ 2 + h4 + (γ ) p ), 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N ,

(3.11)

where C a positive constant independent of τ , h and γ . Therefore, we can construct the compact difference scheme for the problem (1.1)–(1.4) as follows:

γ

s  dj j =0



τ γj

 n  (γ j ) n−k λk Hu i + δx4 uni = H f in , 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N ,

(3.12)

k =0

u 0i = 0, 0 ≤ i ≤ M ,

(3.13)

uk0

(3.14)

= u 0 (tk ),

ukM

= u L (tk ), 0 ≤ k ≤ N ,

where

δx4 un1 =

1 h2

[ηa (tn ) − 2δx2 un1 + δx2 un2 ], δx4 unM −1 =

1 h2

[δx2 unM −2 − 2δx2 unM −1 + ηb (tn )].

4. Stability and convergence In the section, we analyze the stability and convergence of the compact difference scheme (3.12)–(3.14). Before doing that, a prior lemma is given below. Lemma 4.1. For v ∈ V˚ h , it holds that

||H v ||2 ≤ || v ||2 ,

1 3

|| v ||2 ≤ (H v , v ) ≤ || v ||2 ,

1 3

| v |22 ≤ (δx4 v , H v ) ≤ | v |22 .

(4.1)

Proof. According to the definition of the operator H, the proof of the first inequality in (4.1) is straightforward. Noticing that



Hv i = 1 +

h2 6

δx2

vi,

(4.2)

then

(H v , v ) = ( v , v ) +

h2 6

(δx2 v , v ) = || v ||2 −

Using the inverse estimation | v |21 ≤

4 || v ||2 h2

2

1

3

3

(H v , v ) ≥ || v ||2 − || v ||2 = || v ||2 .

h2 6

| v |21 .

(4.3)

in Lemma 2.7 gives

(4.4)

Combing (4.3) and (4.4), the second inequality in (4.1) is obtained. We now prove the third inequality in (4.1). Using Lemma 2.7 again, one gets

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

(δx4 v , H v ) = (δx4 v , v ) +

h2

(δx4 v , δx2 v ) = | v |22 −

6

h2 6

|δx2 v |21 .

65

(4.5)

By the inverse estimation, the above inequality implies that

2

1

3

3

(δx4 v , H v ) ≥ | v |22 − ||δx2 v ||2 = | v |22 ,

(4.6)

where | v |2 = ||δx2 v || has been used. This proof is completed.

2

We are now ready to prove two main result of this paper. Theorem 4.2. Let {uni | 0 ≤ i ≤ M , 0 ≤ n ≤ N } be the solution of the following difference scheme

γ

s  dj j =0



τ γj

 n  (γ j ) n−k λk Hu i + δx4 uni = H f in , 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N ,

(4.7)

k =0

u 0i = φi , 0 ≤ i ≤ M ,

(4.8)

uk0

(4.9)

=

ukM

= 0, 0 ≤ k ≤ N .

Then for h sufficiently small, one has

τ

m 

|un |21 ≤ 3μτ ||φ||2 +

n =1

1 8

T L 2 max || f k ||2 , 1 ≤ m ≤ N .

(4.10)

1≤k≤ N

Proof. Computing inner product of (4.7) with Hun , one gets

γ

s  dj j =0



τ γj

 n  (γ j ) n−k n λk (Hu , Hu ) + (δx4 un , Hun ) = (H f n , Hun ).

(4.11)

k =0

From Lemma 4.1, it follows that

γ

s  dj j =0



τ γj

 n  1 (γ j ) n−k n λk (Hu , Hu ) ≤ − |un |22 + (H f n , Hun ).

(4.12)

3

k =0

By Lemma 2.6 and Lemma 2.8, for h sufficiently small one has (see also [18])

(δx4 un , un ) = |un |22 ≥ 4|un |21 . Thus,

γ

s  dj j =0



τ γj

(4.13)

 n  4 (γ j ) n−k n λk (Hu , Hu ) ≤ − |un |21 + (H f n , Hun ).

(4.14)

3

k =0

Using Cauchy–Schwarz inequality, Lemma 2.7 and Lemma 4.1 yields

(H f n , Hun ) ≤

6

L2

L

24

||Hun ||2 + 2

||H f n ||2 ≤

6

L2

L

24

||un ||2 + 2

|| f n ||2 ≤ |un |21 +

L2 24

|| f n ||2 .

(4.15)

Substituting (4.15) into (4.14) gets

γ

s  dj j =0



τ γj

 n  1 L2 (γ j ) n−k n λk (Hu , Hu ) ≤ − |un |21 + || f n ||2 . k =0

3

(4.16)

24

Summing up the above inequality for n from 1 to m yields

γ

s  dj j =0

Adding

τ γj



 n m  m m  1 n2 L2  n 2 (γ j ) n−k n λk (Hu , Hu ) ≤ − | u |1 + || f || , 1 ≤ m ≤ N . n =1 k =0

3

n =1

24

n =1

μ(Hu 0 , Hu 0 ) on both sides of (4.17), the above inequality implies that

(4.17)

66

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

γ

s  dj j =0



τ γj

 n m  m m  1 n2 L2  n 2 (γ j ) n−k n λk (Hu , Hu ) ≤ − |u |1 + μ||u 0 ||2 + || f || , 1 ≤ m ≤ N , 3

n =0 k =0

24

n =1

(4.18)

n =1

where we have used ||Hu 0 ||2 ≤ ||u 0 ||2 . With the help of Lemma 2.4, one gets

   M −1 m n n n m  m  M −1       (γ j ) (γ j ) (γ j ) n−k n−k n−k n n n λk (Hu , Hu ) = λk h (Hu i )(Hu i ) = h λk (Hu i )(Hu i ) ≥ 0. n =0 k =0

n =0 k =0

i =1

i =1

n =0 k =0

(4.19) Substituting (4.19) into (4.18) yields

τ

m 

|un |21 ≤ 3μτ ||u 0 ||2 +

L2 8

n =1

τ

m 

|| f n ||2 , 1 ≤ m ≤ N .

(4.20)

n =1

Since mτ ≤ T , it follows from the inequality (4.20) that

τ

m 

|un |21 ≤ 3μτ ||φ||2 +

n =1

1 8

T L 2 max || f k ||2 , 1 ≤ m ≤ N . 1≤k≤ N

(4.21)

2

This proof is completed.

From Lemma 2.5, one can see that μτ is bounded when τ is small enough. Hence, Theorem 4.2 implies that the compact difference scheme (3.12)–(3.14) is stable in the L 1 norm with respect to initial value φ and the inhomogeneous term f when h, τ are sufficiently small. Moreover, noticing the equivalence between the | · |1 norm and the || · ||∞ norm, one can conclude that the compact difference scheme (3.12)–(3.14) is also stable in the L ∞ norm under the same conditions. 8, 3

Theorem 4.3. Let u (x, t ) ∈ C x,t ([0, L ] × [0, T ]) be the solution of the problem (1.1)–(1.4), and {uni | 0 ≤ i ≤ M , 0 ≤ n ≤ N } be the solution of the compact difference scheme (3.12)–(3.14). Denote

eni = u (xi , tn ) − uni , 0 ≤ i ≤ M , 0 ≤ n ≤ N . Then

τ

m 

||en ||∞ ≤

n =1

1√ 8

2L 2 T C (τ 2 + h4 + (γ ) p ), 1 ≤ m ≤ N ,

where C is defined in (3.11). Proof. From the derivation of the scheme (3.12)–(3.14), one can get the error system as follows:

γ

s  dj j =0



τ γj

 n  (γ j ) n−k λk He i + δx4 eni = R ni , 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N ,

(4.22)

k =0

e 0i = 0, 0 ≤ i ≤ M ,

(4.23)

ek0

(4.24)

= ekM

= 0, 0 ≤ k ≤ N .

Combing (3.11) and Theorem 4.2 yields

τ

m 

|en |21 ≤

n =1

1 8

L 3 T [C (τ 2 + h4 + (γ ) p )]2 , 1 ≤ m ≤ N .

(4.25)

Using Cauchy–Schwarz inequality and mτ ≤ T one has



m 

2

n

τ ||e ||∞

≤

m 

n =1



τ

2

n =1

m  n =1



n 2

||e ||∞ ≤ T τ

m  n =1

||en ||2∞ .

(4.26)

By the use of Lemma 2.7, it follows from (4.25) that

τ

m  n =1

||en ||2∞ ≤

1 4

Lτ

m  n =1

|en |21 ≤

1 32

L 4 T [C (τ 2 + h4 + (γ ) p )]2 .

(4.27)

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

67

Substituting (4.27) into (4.26) gives



m 

2 n

τ ||e ||∞

1

≤

32

n =1

L 4 T 2 [C (τ 2 + h4 + (γ ) p )]2 .

(4.28)

Namely,

τ

m 

||en ||∞ ≤

n =1

1√ 8

This proof is completed.

2L 2 T C (τ 2 + h4 + (γ ) p ), 1 ≤ m ≤ N .

(4.29)

2

In view of Theorem 4.3 above, one can conclude that the compact difference scheme (3.12)–(3.14) is convergent in the L ∞ norm with the convergence order O (τ 2 + h4 + (γ ) p ). Remark 4.4. Replacing the compact operator H in (3.12) with the identity operator, we can obtain the second-order difference scheme in space for the same initial-boundary value problem (1.1)–(1.4) as follows:

γ

s  dj j =0



τ γj

 n  (γ j ) n−k λk u i + δx4 uni = f in , 1 ≤ i ≤ M − 1, 1 ≤ n ≤ N ,

(4.30)

k =0

u 0i = 0, 0 ≤ i ≤ M ,

(4.31)

uk0

(4.32)

= u 0 (tk ),

ukM

= u L (tk ), 0 ≤ k ≤ N .

By using the arguments similar to those used in Theorem 4.2 and Theorem 4.3, we can prove that the second-order difference scheme (4.30)–(4.32) is stable with respect to initial value φ and the inhomogeneous term f , and convergent with order O (τ 2 + h2 + (γ ) p ) under the L ∞ norm. 5. Numerical results In this section, some numerical results are reported to illustrate our theoretical statements. The maximum error between the exact solution and the difference solution is defined by

E (τ , h, γ ) =

max

0≤i ≤ M ,0≤k≤ N

|U ik − uki |.

The convergence order in time is defined by



Ord1 = log2

E (τ , h, γ )

E (τ /2, h, γ )

for sufficiently small γ and h, and the convergence order in space is defined by



Ord2 = log2

E (τ , h, γ )

E (τ , h/2, γ )

for sufficiently small γ and

τ . All computations are carried out in MATLAB.

Example 5.1. (cf. [10]) Consider the following problem

∂ 4u (x, t ) = f (x, t ), 0 < x < π , 0 < t ≤ 0.5, ∂ x4 u (x, 0) = 0, 0 ≤ x ≤ π ,

Dtw u (x, t ) +

∂ 2u ∂ 2u u (0, t ) = u (π , t ) = 0, (0, t ) = 2 (π , t ) = 0, 0 ≤ t ≤ 0.5, 2 ∂x ∂x

(5.1) (5.2) (5.3)

where w (γ ) = (4 − γ ) and f (x, t ) = 8[6(t 3 − t 2 )/ ln(t ) + t 3 ] sin(x). It is easy to verify that the exact solution of the above problem is u (x, t ) = 8t 3 sin(x). In what follows, the suggested compact difference scheme (3.12)–(3.14) associated with different BVMs will be used to solve this example. The BVMs involved including the GAM-3, ETR2s-4, GBDF-5, GAM-5, ETR2s-6 and TOM-6 mentioned in Section 2.

68

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

Table 2 Maximum errors and convergence orders in time variable for different BVMs. N

ETR2s-4

5 10 20 40 80

GAM-5

TOM-6

Error

Ord1

Error

Ord1

Error

Ord1

2.7085e-2 7.1402e-3 1.8283e-3 4.6240e-4 1.1625e-4

– 1.9234 1.9655 1.9833 1.9919

2.7085e-2 7.1402e-3 1.8283e-3 4.6223e-4 1.1620e-4

– 1.9234 1.9655 1.9838 1.9920

2.7085e-2 7.1402e-3 1.8283e-3 4.6222e-4 1.1619e-4

– 1.9234 1.9655 1.9838 1.9921

Table 3 Maximum errors and convergence orders in space variable for different BVMs. M

ETR2s-4

5 10 20 40

GAM-5

TOM-6

Error

Ord2

Error

Ord2

Error

Ord2

2.1602e-5 1.3249e-6 8.1647e-8 5.3587e-9

– 4.0272 4.0203 3.9294

2.1602e-5 1.3249e-6 8.1643e-8 5.3560e-9

– 4.0272 4.0204 3.9301

2.1602e-5 1.3249e-6 8.1643e-8 5.3551e-9

– 4.0272 4.0204 3.9304

Table 4 Maximum errors and convergence orders in distributed-order variable. GAM-3

ETR2s-4

s

M

N

Error

Ord3

s

M

N

Error

Ord3

6 12 24 48

6 10 17 29

6 17 48 136

1.9251e-2 2.5344e-3 3.2230e-4 4.0519e-5

– 2.9252 2.9752 2.9917

6 12 24 48

6 12 24 48

6 24 96 384

1.9173e-2 1.2752e-3 8.0845e-5 5.0706e-6

– 3.9103 3.9794 3.9949

Table 5 Maximum errors and convergence orders in distributed-order variable. s

M

6 12 24 48

N

6 14 34 81

GBDF-5

6 34 192 1086

GAM-5

Error

Ord3

Error

Ord3

1.9171e-2 6.3886e-4 2.0252e-5 6.3410e-7

– 4.9073 4.9794 4.9972

1.9174e-2 6.3893e-4 2.0254e-5 6.3415e-7

– 4.9073 4.9794 4.9972

We first test the convergence orders in time and space variables. The maximum errors and convergence orders in time “Ord1” for different BVMs with M = 600 and s = 100 are shown in Table 2. Similarly, the maximum errors and convergence orders in space “Ord2” with N = 50000 and s = 200 are listed in Table 3. Based on these data, we can conclude that the suggested compact difference scheme (3.12)–(3.14) has accuracy of order 2 in time and accuracy of order 4 in space with respect to different BVMs. From Theorem 4.3 above, we know that the suggested compact difference scheme (3.12)–(3.14) is convergent with the order O (τ 2 + h4 + (γ ) p ) when a BVM of order p is used. That is, the maximum error has the form of E (τ , h, γ ) = C (τ 2 + h4 + (γ ) p ). It implies that when the step sizes in time, space and distributed-order variables are reduced by a factor of 2 p /2 , 2 p /4 and 21 respectively, one has

E (2− p /2 τ , 2− p /4 h, 2−1 γ ) = C ((2− p /2 τ )2 + (2− p /4 h)4 + (2−1 γ ) p )

= 2− p E (τ , h, γ ). As a result,



log2

E (τ , h, γ ) E (2− p /2 τ , 2− p /4 h, 2−1 γ )

= p.

It shows that we can verify the convergence order in distributed-order variable by using the way. The maximum errors and convergence orders in distributed-order variable “Ord3” for different BVMs with an optimal step size ratio are reported in Tables 4–6, respectively. We can observe that the compact difference scheme (3.12)–(3.14) has accuracy of order p in distributed-order variable when a BVM of order p is used. It confirmed the correctness of our theoretical results.

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

69

Table 6 Maximum errors and convergence orders in distributed-order variable. s

6 12 24 48

M

N

6 17 48 136

ETR2s-6

6 48 384 3072

TOM-6

Error

Ord3

Error

Ord3

1.9174e-2 3.2044e-4 5.0695e-6 7.9226e-8

– 5.9030 5.9821 5.9997

1.9174e-2 3.2044e-4 5.0695e-6 7.9287e-8

– 5.9030 5.9821 5.9986

Table 7 Comparison with the G-C scheme for M = 100, N = 10000 and s = 10. Algorithm

Error

CPU

Algorithm

Error

CPU

C-S ETR2s-4 GBDF-5

5.0130e-7 2.0808e-7 3.1366e-7

1413s 1362s 1301s

GAM-5 ETR2s-6 TOM-6

8.3610e-8 4.1938e-9 4.1068e-9

1396s 1329s 1405s

Finally, the suggested scheme is compared with the method obtained by using the composite Simpson formula to approximate the integral term in the distributed-order derivative. For simplicity, the latter method is denoted by C-S scheme. It can be easily proved that the C-S scheme has accuracy order of 2, 4 and 4 in time, space and distributed-order variables, respectively. In contrast to the C-S scheme, the biggest advantage of our method is the high order accuracy in the distributed-order variable. To reflect this, we fix M = 100 and N = 10000. The comparison results is given in Table 7. From which one can see that the six schemes cost almost the same CPU time (second), while the suggested compact difference scheme (3.12)–(3.14) can achieve more higher accuracy than the C-S scheme. 6. Conclusion In the current work, we derived a class of high order compact difference schemes for solving the fourth-order time fractional sub-diffusion of distributed order. The derived method is based on a new quadrature rule which originated by boundary value method. Using the discrete energy method, we proved the suggested scheme is stable with respect the initial value and the inhomogeneous term, and convergent in L ∞ norm with the convergence order O (τ 2 + h4 + (γ ) p ) when a boundary value method of order p was used, where τ , h and γ are the step sizes in time, space and distributed-order variables, respectively. Numerical results verified the effectiveness and high accuracy of our method. It is evident that this idea can be applied to solve other related problems. Acknowledgements We would like to express our sincere thanks to the referee and the editors for their valuable comments and suggestions, which have let to great improvements over the original manuscript. This work was supported by the National Natural Science Foundation of China (Grant No. 11571128). References [1] L. Aceto, C. Magherini, P. Novati, Fractional convolution quadrature based on generalized Adams methods, Calcolo 51 (2014) 441–463. [2] L. Brugnano, D. Trigiante, Convergence and stability of boundary value methods for ordinary differential equations, J. Comput. Appl. Math. 66 (1996) 97–109. [3] L. Brugnano, D. Trigiante, Solving Differential Equations by Multistep Initial and Boundary Value Methods, Gordan and Breach, Amsterdam, 1998. [4] C. Çelik, M. Duman, Crank–Nicolson method for the fractional diffusion equation with the Riesz fractional derivative, J. Comput. Phys. 231 (2012) 1743–1750. [5] A. Chechkin, R. Gorenflo, I. Sokolov, V. Gonchar, Distributed order time fractional diffusion equation, Fract. Calc. Appl. Anal. 6 (2003) 259–280. [6] H. Chen, C. Zhang, Boundary value methods for Volterra integral and integro-differential equations, Appl. Math. Comput. 218 (2011) 2619–2630. [7] K. Diethelm, N. Ford, Numerical solution methods for distributed order differential equations, Fract. Calc. Appl. Anal. 4 (2001) 531–542. [8] K. Diethelm, N. Ford, Numerical analysis for distributed-order differential equations, J. Comput. Appl. Math. 225 (2009) 96–104. [9] G. Gao, 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–948. [10] G. Gao, H. Sun, Z. Sun, Some high-order difference schemes for the distributed-order differential equations, J. Comput. Phys. 298 (2015) 337–359. [11] X. Hu, L. Zhang, On finite difference methods for fourth-order fractional diffusion-wave and subdiffusion systems, Appl. Math. Comput. 218 (2012) 5019–5034. [12] Z. Jiao, Y. Chen, I. Podlubny, Distributed-Order Dynamic Systems: Stability, Simulation, Applications and Perspectives, Springer, London, 2012. [13] A. Kochubei, Distributed order calculus and equations of ultraslow diffusion, J. Math. Anal. Appl. 340 (2008) 252–281. [14] X. Li, B. Wu, A numerical method for solving distributed order diffusion equations, Appl. Math. Lett. 53 (2016) 92–99. [15] M. Meerschaert, E. Nane, P. Vellaisamy, Distributed-order fractional diffusions on bounded domains, J. Math. Anal. Appl. 379 (2011) 216–228. [16] M. Morgado, M. Rebelo, Numerical approximation of distributed order reaction–diffusion equations, J. Comput. Appl. Math. 275 (2015) 216–227. [17] M. Naber, Distributed order fractional sub-diffusion, Fractals 12 (2004) 23–32. [18] K. Omrani, F. Abidi, T. Achouri, N. Khiari, A new conservative finite difference scheme for the Rosenau equation, Appl. Math. Comput. 201 (2008) 35–43.

70

[19] [20] [21] [22]

M. Ran, C. Zhang / Applied Numerical Mathematics 129 (2018) 58–70

I. Sokolov, A. Chechkin, J. Klafter, Distributed-order fractional kinetics, Acta Phys. Pol. B 35 (2004) 1323. Z. Sun, Numerical Methods of Partial Differential Equations, second ed., Science Press, Beijing, 2009 (in Chinese). Z. Sun, X. Wu, A fully discrete difference scheme for a diffusion-wave system, Appl. Numer. Math. 56 (2006) 193–209. Y. Tian, H. Zhou, W. Deng, A class of second order difference approximations for solving space fractional diffusion equations, Math. Comput. 84 (2015) 1703–1727. [23] H. Ye, F. Liu, V. Anh, Compact difference scheme for distributed-order time-fractional diffusion-wave equation on bounded domains, J. Comput. Phys. 298 (2015) 652–660.