A gWSGL numerical scheme for generalized fractional sub-diffusion problems

Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Research paper

A gWSGL numerical scheme for generalized fractional sub-diffusion problems Xuhao Li, Patricia J.Y. Wong∗ School of Electrical and Electronic Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore

a r t i c l e

i n f o

Article history: Received 28 April 2019 Revised 29 July 2019 Accepted 3 September 2019 Available online 9 September 2019 Keywords: Generalized fractional derivative Sub-diffusion problem Generalized weighted shifted Grünwald-Letnikov approximation Numerical scheme

a b s t r a c t In this paper, an efficient numerical scheme is constructed for a generalized fractional subdiffusion problem using a newly proposed generalized weighted shifted Grünwald-Letnikov (gWSGL) approximation for the generalized fractional derivative. The solvability, stability and convergence of the numerical scheme are analyzed using the discrete energy method. It is proven that the temporal convergence order is 2 and this is the best result to date. Simulation is further carried out to demonstrate the accuracy of the proposed numerical scheme. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Fractional problems have gathered increasing interests during the past four decades since many physical phenomena can be modeled much more accurately using fractional operators, see [14–18]. The classical fractional problem, e.g. fractional sub-diffusion and diffusion-wave problems discussed in [14–18], normally involves Caputo or Riemann-Liouville fractional derivative (or integral) which can be considered as a kind of generalization of integer order derivative (or repeated integral). Recently, Agrawal [3] did a survey on some generalizations of fractional operators [9,11] and presented a new generalized fractional operator involving a scale function z(t) and a weight function w(t ). The scale function z(t) can stretch or contract the time domain, this may be necessary to accurately capture the phenomenon of interest over a desired range of time, which could be within micro/nano-seconds, or over several decades. On the other hand, the weight function w(t ) allows an event to be weighed differently at different time points, for example the modelling of the memory of a child may require a heavier weight at current time point, whereas the same for an older person may require more weight on the past. As shown in [3], one can obtain well known fractional derivatives in the literature by selecting specific scale and weight functions, as such the generalized fractional operator in [3] indeed unifies many fractional operators proposed in the literature. In this sense, this generalized fractional operator is able to relate to different physical phenomena for different choices of z(t) and w(t ), just like how the different phenomena are related to different types of fractional derivatives. As an example, when z(t ) = t and w(t ) ≡ 1, the generalized fractional operator gives the classical Caputo fractional derivative, which is used to describe anomalous diffusion [16]. Though this new generalization of Caputo or Riemann-Liouville fractional derivative or integral brings convenience in expression, it is very challenging to present analytical solutions of equations involving generalized fractional operators except ∗

Corresponding author. E-mail addresses: [email protected] (X. Li), [email protected] (P.J.Y. Wong).

https://doi.org/10.1016/j.cnsns.2019.104991 1007-5704/© 2019 Elsevier B.V. All rights reserved.

2

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

for a small class of problems. There is some pioneer work [10,24] contributed to this topic. In fact, Kim et al. [10] have proved the existence and uniqueness of a class of initial value problem under appropriate assumptions, moreover the analytical solution is given by applying operational method. Different from the method used in [10], in [24] a generalized fractional sub-diffusion equation with zero boundary conditions is analyzed in a rather simple and straightforward way by using the method of separation of variables. However, for most generalized fractional problems, in general it is difficult if not impossible to find the explicit analytical solutions. It is therefore necessary to propose efficient numerical schemes for such problems in practice. There has been some numerical investigation on generalized fractional problems such as fractional Burgers equation [22], fractional sub-diffusion equation [24], fractional advection-diffusion equation [25], generalized van der Pol equation [21] as well as oscillator equation [23]. To get these generalized fractional equations, one can simply replace the original fractional operator by the new generalized fractional operator (as in [22,24,25]) or by an even more general operator involving a kernel function (as in [4,21,23]). The main difficulty in deriving the numerical scheme for a generalized fractional equation lies in the approximation of the generalized fractional operator. To tackle this, Xu et al. [24] have proposed a useful approximation that generalizes the classical L1 discretization [6] in the sense that it coincides with L1 when z(t ) = t and w(t ) ≡ 1. Thereafter, this generalized L1 approximation has been applied successfully to problems in [22,25]. As for generalized van der Pol equation [21] and oscillator problem [23], the authors follow a similar procedure as in [24] to approximate the generalized fractional derivative, that is, using linear approximation for the first derivative of the integrand in each subinterval and then compute the integral. However, this is not easy in practice if the kernel function is complex. In this work, we shall consider a generalized fractional sub-diffusion equation of order α , where 0 < α < 1. As demonstrated in [2,16], from a physical point of view, such an equation basically describes a slow diffusion process as compared to ordinary diffusion (α = 1 ). In a slow diffusion process, the mean square displacement of particles from the original starting site is not linear in time but instead follows a generalized Fick’s second law. Indeed, the sub-diffusion motion is characterized by the following asymptotic behavior of the mean square displacement (assume that z(t) → ∞ as t → ∞)

x(z(t )) ∼ z(t )α ,

t→∞

where α (0 < α < 1) is the anomalous diffusion exponent. Let v˜ (x, z(t )) be the probability density function that describes this sub-diffusion motion. Similar to [26], we have C α 0 Dz ˜

v(x, z(t )) =

∂ 2 v˜ (x, z(t )) , ∂ x2

t>0

(1.1)

˜ (x, z(t )) is the classical Caputo fractional derivative defined by where C0 Dα zv C α 0 Dz ˜

v(x, z(t )) =

1 (1 − α )

 z(t ) 0

1

(z(t ) − z˜)α

∂ v˜ (x, z˜) dz˜. ∂ z˜

Let z˜ = z(s ) and v(x, t ) = v˜ (x, z(t )). It is easy to check that C α 0 Dz ˜

v(x, z(t )) = C0 Dtα;[z(t ),1] v(x, t )

(see (1.4) for the definition of C0 Dtα;[z(t ),1] v(x, t )). Hence, (1.1) yields a generalized fractional sub-diffusion equation C α 0 Dt ;[z(t ),1]

v(x, t ) =

∂ 2 v(x, t ) , ∂ x2

t > 0.

(1.2)

In this paper, we shall investigate a more general formulation than (1.2) that is first studied in [24]

C

α

v(x, t ) = avxx (x, t ) + g(x, t ), (x, t ) ∈ (0, 1 ) × (0, 1 ) v(x, 0 ) = v0 (x ), x ∈ ( 0, 1 ) v(0, t ) = g1 (t ), v(1, t ) = g2 (t ), t ∈ ( 0, 1 ) 0 Dt ;[z(t ),w(t )]

(1.3)

where 0 < α < 1 is the fractional order, a ( > 0) is the diffusion constant, v0 (x ), g1 (t ), g2 (t ),g(x, t) are continuous functions in corresponding domain, z(t) is the scale function and w(t ) is the weight function. The generalized Caputo fractional derivative C Dα v(x, t ) is defined by [3] 0 t ;[z (t ),w(t )] C α 0 Dt ;[z(t ),w(t )]

v(x, t ) =

[w(t )]−1 (1 − α )



0

t

1 [z(t ) − z(s )]α

∂ (w(s )v(x, s ) ) ds. ∂s

(1.4)

We shall first develop a new approximation for the generalized fractional derivative, namely the generalized weighted shifted Grünwald-Letnikov (gWSGL) approximation. Using the gWSGL approximation, we shall construct an efficient numerical scheme for the above generalized fractional sub-diffusion problem. Let u(x, t ) = w(t )v(x, t ) − w(0 )v0 (x ). It is easy to see that (1.3) is equivalent to the following generalized fractional problem with the same scale function z(t) and unity weight function

C

α

(x, t ) = auxx (x, t ) + f (x, t ), u(x, 0 ) = 0, x ∈ ( 0, 1 ) u(0, t ) = φ1 (t ), u(1, t ) = φ2 (t ), 0 Dt ;[z(t ),1] u

(x, t ) ∈ (0, 1 ) × (0, 1 ) (1.5) t ∈ ( 0, 1 )

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

3

where φ1 (t ) = w(t )g1 (t ) − w(0 )v0 (0 ), φ2 (t ) = w(t )g2 (t ) − w(0 )v0 (1 ) and f (x, t ) = aw(0 )v

0 (x ) + w(t )g(x, t ). Once the numerical solution of (1.5) is obtained, the numerical solution of the original problem (1.3) is easily computed from the relation u(x, t ) = w(t )V (x, t ) − w(0 )V0 (x ). The paper is organized as follows. In Section 2, we shall first develop a new generalized weighted shifted Grünwald-Letnikov (gWSGL) approximation for the generalized fractional derivative. Then, a numerical scheme for the equivalent generalized fractional problem (1.5) is constructed. In Section 3, we shall establish a priori estimate of the solution of (1.5) with zero boundary conditions through energy methodology. This a priori estimate is then used to analyze the stability and convergence in Section 4. In section 5, four numerical experiments are presented to demonstrate the efficiency and accuracy of the proposed numerical scheme. Finally, we give a brief conclusion in Section 6. 2. gWSGL numerical scheme In this section, we shall first develop a generalized weighted shifted Grünwald-Letnikov (gWSGL) approximation for the generalized fractional derivative. This new approximation is then used to derive a numerical scheme for problem (1.5). Throughout, the scale function z(t) is assumed to be strictly increasing. Let Z : t → z(t ) be a mapping from [0,1] to [z(0), z(1)]. To begin, let

P : 0 = x0 < x1 < · · · < xM = 1 be a uniform partition of the space interval [0,1] with step size h =

1 M.

Next, let

P : 0 = t0 < t1 < · · · < tN = 1 (0 ) be a mesh of the time interval [0,1] such that τz ≡ z(tk ) − z(tk−1 ), 1 ≤ k ≤ N is a constant. Obviously, τz = z(1 )−z . N For any given function y(x, t), we denote y(xj , tn ) by ynj and y(xj , t) by yj for fixed t. Let C,r ([0, 1] × [0, 1]) be a space j

of functions {y(x, t)} such that the derivatives Dix Dt y(x, t ), i = 0, 1, · · · , , j = 0, 1, · · · , r are continuous in [0, 1] × [0, 1]. Moreover, denote

Uh = {y = (y0 , y1 , · · · , yM ) and C α +m (IR ) =

    φ (t ) 

∞ −∞

| y0 = yM = 0}

     (1 + |ξ |α+m )φˆ (ξ )dξ < ∞ ,

∞ where φˆ (ξ ) = −∞ φ (t )e−it ξ dt is the Fourier transform of φ (t). For any y, yˆ ∈ Uh , we define

δx y j+ 12 =

1 y j+1 − y j , h

δx2 y j =

1 h



and

(y, yˆ ) = h

M−1

y j yˆ j ,

(y, yˆ )1 = h

j=1

Obviously, y =

M−1

 δx y j+ 12 − δx y j− 12 ,



j=0

(y, y ) is a norm and |y|1 =

δx y j+ 12



(2.1)

 δx yˆ j+ 12 .

(2.2)

(y, y )1 is a kind of semi-norm. We also define y ∞ = max1≤ j≤M−1 |y j |.

2.1. gWSGL approximation To proceed to the development of the gWSGL approximation, we first introduce the definition of the generalized RiemannLiouville fractional derivative as in [3]. Definition 2.1 [3]. For 0 < α < 1 and fixed x ∈ (0, 1), the generalized Riemann-Liouville fractional derivative of u(x, t) with scale function z(t) and weight function w(t ) is defined by RL α 0 Dt ;[z(t ),w(t )] u

(x, t ) =

[w(t )]−1 1 D (1 − α ) z (t ) t



t 0



w(s )z (s )u(x, s ) ds . [z(t ) − z(s )]α

(2.3)

Next, we shall present some lemmas. The first one relates a subspace of C (2+m ) to C α +m . Lemma 2.1. [7,27] Let 0 < α < 1. If φ (t ) ∈ { (t ) ∈ C (2+m ) (IR ) | ,

, · · · , (m+2 ) ∈ L1 (IR )}, then we have φ (t ) ∈ C α +m (IR ). The next lemma gives the typical weighted shifted Grünwald-Letnikov (WSGL) approximation for the classical RiemannLiouville fractional derivative, which is defined by RL α −∞ Dt

φ (t ) =

RL α −∞ Dt ;[z(t )=t ,w(t )≡1]

1 φ (t ) = D (1 − α ) t



t

−∞

 φ (s ) ds . (t − s )α

(2.4)

4

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

Lemma 2.2. [5,19] Let 0 < α < 1 and φ (t ) ∈ L1 (IR ) ∩ C α +2 (IR ). The following approximation holds RL α −∞ Dt

α − 2q α 2p − α α A φ (t ) + A φ (t ) + O(τ 2 ), 2( p − q ) τ ,p 2( p − q ) τ ,q

φ (t ) =

(2.5)

where p, q ( = p) are integers, τ is the step size,

Aατ,p φ (t ) = τ −α

∞

rα ,k

k=0



φ (t − (k − p)τ ),



rα ,0 = 1 and rα ,k = 1 − α +1 rα ,k−1 for k ≥ 1. k We are now ready to construct the gWSGL approximation for the generalized Riemann-Liouville fractional derivative defined in (2.3). We write

u(x, t ) = u(x, Z −1 (z )) ≡ u¯ (x, z ), and extend the function u¯ as follows:

⎧ ⎪ ⎨0, u¯ (x, z ), u¯ e (x, z ) = ⎪ ⎩ψ (z ),

z < z (0 ) z (0 ) ≤ z ≤ z (1 ) z ( 1 ) < z ≤ 2z ( 1 ) z > 2z ( 1 )

0,

where ψ (z) is a smooth function that can be determined later. Theorem 2.1 (gWSGL approximation). Let 0 < α < 1. Suppose that u(·, Z −1 (z )) ≡ u¯ (·, z ) ∈ C (4 ) ([z(0 ), z(1 )] ) and ∂ k u¯ |z=z(0 ) = ∂z 0, k = 0, 1, 2, 3, 4. Then, the generalized Riemann-Liouville fractional derivative with unity weight function at (xj , tn ) has the following approximation k

RL α 0 Dt ;[z(t ),1] u

( x j , tn ) =

α − 2q α n 2 p − α α n G u + G u + O(τz2 ), 2( p − q ) τz ,p j 2( p − q ) τz ,q j

(2.6)

(0 ) where p, q ( = p) are integers, τz = z(1 )−z is the step size, N

Gατz ,p unj = τz−α



n+ p

rα ,k unj −k+ p,

k=0

(2.7)



rα ,0 = 1 and rα ,k = 1 − α +1 rα ,k−1 for k ≥ 1. k Proof. In view of the assumptions, we can choose a function ψ (z) such that u¯ e (·, z ) ∈ C (4 ) (IR ). In fact, we can apply Hermite interpolation to get a ψ (z) that ensures the continuity at z = z(1 ) and z = 2z(1 ). Obviously, u¯ e (·, z ) ∈ L1 (IR ). Also, using Lemma 2.1, we get u¯ e (·, z ) ∈ C α +2 (IR ). Therefore, u¯ e (·, z ) ∈ L1 (IR ) ∩ C α +2 (IR ). (0 ) Note that P

: z(0 ) = z(t0 ) < z(t1 ) < · · · < z(tN ) = z(1 ) is the uniform mesh of [z(0), z(1)] with step size τz = z(1 )−z . N Hence, denoting z(tn ) by zn , we apply Lemma 2.2 to get RL α¯ e −∞ Dz u

( x j , zn ) =

α − 2q α 2p − α α A u¯ e (x j , zn ) + A u¯ e (x j , zn ) + O(τz2 ), 2( p − q ) τz ,p 2( p − q ) τz ,q

(2.8)

where

Aατz ,p u¯ e (x j , zn ) = τz−α

n+ p

rα ,k u¯ (x j , zn−k+ p ) = τz−α

k=0

n+ p

rα ,k u(x j , tn−k+ p ) = Gατz ,p unj .

(2.9)

k=0

Since u¯ e (x, z ) = 0 for z < z(0), for z ≥ z(0), we have, noting the definition (2.4), RL α¯ e −∞ Dz u

(x j , z ) =

RL α ¯ z0 Dz u

(x j , z ) =

1 d (1 − α ) dz



z z0



u¯ (x j , y ) dy . ( z − y )α

(2.10)

Using the transformation y = z(s ) and noting that z = z(t ) and u¯ (x, z(t )) = u(x, t ), we further get RL α ¯ z0 Dz u

(x j , z ) =

1 d (1 − α ) dt



z(t )

z0





u¯ (x j , z(s )) dt dz(s ) [z(t ) − z(s )]α dz(t )



u (x j , s ) 1 z (s )ds

α z (t ) 0 [z (t ) − z (s )]  t  u (x j , s ) 1 1 d

= z (s )ds (1 − α ) z (t ) dt 0 [z(t ) − z(s )]α =

1 d (1 − α ) dt

=

RL α 0 Dt ;[z(t ),1] u

t

( x j , t ),

(2.11)

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

5

where we have used the definition (2.3) in the last equality. It follows from (2.10) and (2.11) that RL α¯ e −∞ Dz u

( x j , zn ) =

RL α ¯ z0 Dz u

(x j , z ) =

RL α 0 Dt ;[z(t ),1] u

( x j , tn ).

(2.12) 

The approximation (2.6) now follows directly from (2.12), (2.8) and (2.9).

 k  Remark 2.1. The conditions in Theorem 2.1, namely, u(·, Z −1 (z )) ≡ u˜ (·, z ) ∈ C (4 ) [z(0 ), z(1 )] and ∂ k u˜ ∂z L1 (IR )

= 0, k =

z=z (0 ) ∩ C α +2 (IR ), which is sufficient enough for us to

0, 1, 2, 3, 4, can be replaced by the much weaker condition u˜e (·, z ) ∈ apply Lemma 2.2 in the proof of the theorem. The much weaker condition ensures that the gWSGL approximation (2.6) is applicable to many more functions.

Remark 2.2. In practice, since u is computed in consecutive time steps t1 , t2 , , it is necessary to take ( p, q ) = (0, −1 ) or (−1, 0 ) in Theorem 2.1. In this case, the gWSGL approximation (2.6) yields RL α 0 Dt ;[z(t ),1] u

(x j , tn ) = τz−α

n

λα,k unj −k + O(τz2 ),

(2.13)

k=0

where

α+2

λα,0 =

2

r α ,0 =

α+2 2

α+2

λα,k =

,

rα ,k −

2

α

r , k ≥ 1. 2 α ,k−1

(2.14)

Next, we shall derive the relation between the generalized Caputo fractional derivative (1.4) and the generalized Riemann-Liouville fractional derivative (2.3). This relation is very important in constructing the numerical scheme. Theorem 2.2. Let 0 < α < 1. Suppose u(x, t) is sufficiently smooth as a function of t. We have the relation RL α 0 Dt ;[z(t ),w(t )] u

(x, t ) =

[w(t )]−1 w(0 )u(x, 0 ) + (1 − α ) [z(t ) − z(0 )]α

Proof. Using the definition (2.3), we write RL α 0 Dt ;[z(t ),w(t )] u

(x, t ) = =

[w(t )]−1 1 D (1 − α ) z (t ) t

−

t 0



[z(t ) − z(s )]1−α w(s )u(x, s )d 1−α

0



t 0



w(s )u(x, s )d

(2.15)



[z(t ) − z(s )]1−α 1−α

 .

(2.16)



[z(t ) − z(s )]1−α t  + 1−α 0  t [z(t ) − z(0 )]1−α = w(0 )u(x, 0 ) + 1−α 0 = − w(s )u(x, s )

(x, t ).

w(s )z (s )u(x, s ) ds [z(t ) − z(s )]α

t

 

[w(t )] 1 D − (1 − α ) z (t ) t −1

Applying integration by parts, we have





C α 0 Dt ;[z(t ),w(t )] u



∂ (w(s )u(x, s )) [z(t ) − z(s )]1−α ds ∂s 1−α 0 ∂ (w(s )u(x, s )) [z(t ) − z(s )]1−α ds, ∂s 1−α t

which, upon differentiating with respect to t, yields

 





[z(t ) − z(s )]1−α 1−α 0  t w(0 )u(x, 0 )

1

(t ) = z ( t ) + z α [z(t ) − z(0 )]α 0 [z (t ) − z (s )]

Dt −

t

w(s )u(x, s )d

∂ (w(s )u(x, s )) ds. ∂s

(2.17)

Now, substituting (2.17) into (2.16) and noting the definition of the generalized Caputo fractional derivative (1.4), we obtain (2.15) immediately.  Remark 2.3. If u(x, 0 ) = 0, then Theorem 2.2 leads to the equality of the generalized Riemann-Liouville and the generalized Caputo fractional derivatives. Since u(x, 0 ) = u¯ (x, z(0 )), we see that the conditions of Theorem 2.1 actually include the condition u(x, 0 ) = 0. Hence, noting Remark 2.3, the gWSGL approximation (2.6) is also applicable to the generalized Caputo fractional derivative. Combining with Remark 2.2, we have the following corollary. Corollary 2.1 (gWSGL approximation). Let 0 < α < 1. Suppose that u(·, Z −1 (z )) ≡ u¯ (·, z ) ∈ C (4 ) ([z(0 ), z(1 )] ) and ∂ k u¯ |z=z(0 ) = ∂z 0, k = 0, 1, 2, 3, 4. Then, the generalized Caputo fractional derivative with unity weight function at (xj , tn ) has the following approximation k

C α 0 Dt ;[z(t ),1] u

(x j , tn ) = τz−α

n k=0

λα,k unj −k + O(τz2 ),

(2.18)

6

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

where λα ,k is given in (2.14). Remark 2.4. In view of Remarks 2.1 and 2.3, the conditions in Corollary 2.1 can be replaced by the weaker conditions u(x, 0 ) = 0 and u˜e (·, z ) ∈ L1 (IR ) ∩ C α +2 (IR ). 2.2. Numerical scheme We are now ready to derive a numerical scheme for the equivalent problem (1.5) using the gWSGL approximation and finite difference method. First, we discretize (1.5) at (xj , tn ) to get C α 0 Dt ;[z(t ),1] u

(x j , tn ) = auxx (x j , tn ) + f (x j , tn ).

(2.19)

Noting Corollary 2.1, we apply the gWSGL approximation (2.18) to the generalized Caputo fractional derivative on the left side of (2.19). Further, the classical finite difference method gives for u(x, · ) ∈ C(4) [0, 1],

uxx (x j , tn ) =

unj+1 − 2unj + unj−1

+ O(h2 ) = δx2 unj + O(h2 ).

h2

(2.20)

Substituting (2.18) and (2.20) into (2.19), we obtain the following scheme after omitting the small terms

⎧ n ⎪ ⎪ λα,k unj −k = aδx2 unj + f jn , ⎨τz−α

1≤ j ≤M−1 (2.21)

k=0

⎪ u0 = 0, 0≤ j≤M ⎪ ⎩ nj u0 = φ1 (tn ), unM = φ2 (tn ), where τz =

z (1 )−z (0 ) N

1≤n≤N

and λα ,k is given in (2.14).

Remark 2.5. The proposed numerical scheme (2.21) has a unique solution at each time level t = tn , 1 ≤ n ≤ N since the coefficient matrix

⎛

−2 − ahτ α λα ,0 z ⎜ 1 2

⎜

a⎜ − 2⎜ h ⎜

⎝

⎞

1 −2 − ahτ α λα ,0 z .. . 2

1 .. . 1

..

⎟ ⎟ ⎟ ⎟ ⎟ ⎠

.

−2 − ahτ α λα ,0 z 1 2

1 −2 − ahτ α λα ,0 z 2

(M−1 )×(M−1 )

is easily seen to be strictly diagonally dominant for any α ∈ (0, 1). To implement the numerical scheme (2.21), we first solve for u1j , 1 ≤ j ≤ M − 1; thereafter, we can get u2j , 1 ≤ j ≤ M − 1. Continuing in this manner, we can obtain the numerical solution unj , 1 ≤ j ≤ M − 1, 1 ≤ n ≤ N of the problem (1.5). The numerical solution {vnj } of the original problem (1.3) is then computed from the relation u(x, t ) = w(t )v(x, t ) − w(0 )v0 (x ). 3. A priori estimate In this section, we shall establish a priori estimate of the solution of the numerical scheme (2.21) with un0 = unM = 0 in ‘average norm’. This a priori estimate will be used in the stability and convergence analysis in the next section. It is noted that a priori estimates in ‘average norm’ have also been sighted in [7,8]. We first present two lemmas which are essential in getting the a priori estimate. Lemma 3.1 [20]. For any positive integer i and real sequence {yn }, we have



i n n=0



λα,k yn−k yn ≥ 0,

(3.1)

k=0

where λα ,k is given in (2.14). Lemma 3.2. [1,6,7] For u ∈ Uh , we have

1

u ≤ √ |u|1 6

and

1 2

u ∞ ≤ |u|1 .

(3.2)

The main result of this section is the following. Theorem 3.1. (A priori estimate) Let {unj } denote the solution of the following numerical scheme

⎧ n ⎪ ⎪ λα,k unj −k = aδx2 unj + f jn , ⎨τz−α k=0

⎪ u0 = 0, ⎪ ⎩ nj

u0 = 0,

0≤ j≤M unM = 0,

1≤ j ≤M−1 (3.3)

1 ≤ n ≤ N.

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

We have

τz

i

u ∞ ≤ n

  

!

[z(1 ) − z(0 )]

i τz

24a2

n=1

f n 2 + τ

1 −α z

n=1

λα,0 2a

7

" ( u0 , u0 ) , 1 ≤ i ≤ N

(3.4)

where λα ,0 = α +2 2 . Proof. First, we multiply both sides of the first equation in (3.3) by hunj , then sum j from 1 to M − 1, and sum n from 1 to i (where i ∈ {1, 2, , N} is fixed). Noting the notation in (2.2), we obtain

τz−α

i n

i  i  λα,k un−k , un = d δx2 un , un + ( f n , un ),

n=1 k=0

n=1

(3.5)

n=1

or equivalently

τz−α

i n

i  i   λα,k un−k , un = d δx2 un , un + ( f n , un ) + τz−α λα,0 u0 , u0 .

n=0 k=0

n=1

(3.6)

n=1

Applying Lemma 3.1 to the term in the left side of (3.6) gives

τ

−α z

i n



λα,k u

n−k

,u

n



=h

n=0 k=0



M−1

τ

−α z

i n



λ

n−k n α ,k u j u j

≥ 0.

(3.7)

n=0 k=0

j=1





Since un0 = unM = 0, obviously we have un ∈ Uh . Therefore, by simple calculation, it is found that δx2 un , un = −|un |21 . Using this result and (3.7), it follows from (3.6) that

d

i

|un |21 ≤

n=1

η,

i

 ( f n , un ) + τz−α λα,0 u0 , u0 .

Next, using the inequality xy ≤

( f n , un ) ≤ Setting η =

(3.8)

n=1

1 , 6d

η 2

f n 2 +

1 2 2 (x

+ y2 ) and the first inequality in (3.2) of Lemma 3.2, we find for any positive constant

1 η 1 1 n 2 · |u |1 . un 2 ≤ f n 2 + 2η 2 2η 6

the above relation becomes

( f n , un ) ≤

1 d f n 2 + |un |21 . 12a 2

(3.9)

Using (3.9) in (3.8) and then applying the second inequality in (3.2) of Lemma 3.2, it is found that

τz

i

un 2∞ ≤

n=1

i τz

24a2

f n 2 + τz1−α

λα,0

n=1

2a

( u0 , u0 ),

Finally, applying Cauchy-Schwarz inequality gives



τz

i

2

un ∞



≤

τz

n=1

i n=1



1

τz !

≤ [z ( 1 ) − z ( 0 ) ]

i

1 ≤ i ≤ N.

(3.10)



un 2∞

n=1 i τz

24a2

n=1

f + τ n 2

1 −α z

λα,0 2a

" (u , u ) , 0

0

(3.11)

where the fact iτz ≤ Nτz = z(1 ) − z(0 ) and (3.10) are used in the last inequality. This completes the proof. Remark 3.1. From the proof of Theorem 3.1, we notice that the a priori estimate (3.4) holds for it is not necessary to have u0 = 0, i.e., zero initial condition.

un

∈ Uh and for any

 u0 .

Here,

4. Stability and convergence Using the a priori estimate established in Section 3, we shall analyze the stability and convergence of the proposed numerical scheme (2.21) rigorously in terms of ‘average norm’. It is noted that stability and convergence analyses in terms of ‘average norm’ have also been sighted in [7,8].

8

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

4.1. Stability analysis Let {unj } be the solution obtained from the numerical scheme (2.21) and {ynj } be the solution of the corresponding perturbed system with non-homogeneous term and initial condition

⎧ n ⎪ ⎪ λα,k ynj −k = aδx2 ynj + f¯jn , ⎨τz−α

1≤ j ≤M−1 (4.1)

k=0

⎪ y0 = θ j , 0≤ j≤M ⎪ ⎩ nj y0 = φ1 (tn ), ynM = φ2 (tn ),

1≤n≤N

where θ0 = θM = 0. We claim that the numerical scheme (2.21) is robust to perturbation of non-homogeneous and initial data. To be specific, we have the following theorem. Theorem 4.1. (Stability) Let  nj = unj − ynj . Then, { nj } satisfies the following inequality

τz

i

 ∞ ≤ n

  

!

[z ( 1 ) − z ( 0 ) ]

n=1

" i τz 1−α λα ,0 n 2 2 ˆ f + τz θ , 2

24a

2a

n=1

1≤i≤N

(4.2)

where fˆ = f − f¯. Proof. From (2.21) and (4.1), it is obvious that { nj } is the solution of the following system

⎧ n ⎪ ⎪ λα,k  nj −k = aδx2  nj + fˆjn , ⎨τz−α

1≤ j ≤M−1 (4.3)

k=0

⎪  0 = −θ j , 0≤ j≤M ⎪ ⎩ nj 0 = 0, Mn = 0, 1 ≤ n ≤ N. Noting that  n ∈ Uh and Remark 3.1, we apply Theorem 3.1 to get

τz

i

 ∞ ≤ n

  

!

[z ( 1 ) − z ( 0 ) ]

n=1

" i τz 1−α λα ,0 n 2 2 ˆ f + τz θ , 2

24a

2a

n=1

1 ≤ i ≤ N.

(4.4)

This shows that the numerical scheme (2.21) is robust to perturbation of non-homogeneous and initial data.



4.2. Convergence analysis We shall obtain the local truncation error of the numerical scheme (2.21) which is vital in the convergence analysis. Let U(x, t) be the exact solution of the problem (1.5). The local truncation error will be incurred if we replace unj in (2.21) by U nj = U (x j , tn ), in fact we have

⎧ n ⎪ ⎪ λα,kU jn−k = aδx2U jn + f jn + Tjn , ⎨τz−α

1≤ j ≤M−1 (4.5)

k=0

⎪ U 0 = 0, 0≤ j≤M ⎪ ⎩ jn n U0 = φ1 (tn ), UM = φ2 (tn ), where

T jn

1≤n≤N

denotes the local truncation error at (xj , tn ).

k Lemma 4.1. Suppose that U (x, Z −1 (z )) = U¯ (x, z ) ∈ C 4,4 ([0, 1] × [z(0 ), z(1 )] ) and ∂ k U¯ |z=z(0 ) = 0, k = 0, 1, 2, 3, 4. Then, the lo∂z cal truncation error T jn satisfies

T jn = O(τz2 + h2 ),

1 ≤ j ≤ M − 1, 1 ≤ n ≤ N

(4.6)

(0 ) where h is the spatial step size and τz = z(1 )−z . N

Proof. From the first equation in (4.5) and the fact that U(x, t) is the exact solution of (1.5), we get

T jn =

τz−α

n

λα,kU jn−k − aδx2U jn − f jn

k=0

=

τz−α

n k=0

λα,kU jn−k − aδx2U jn − f jn −

#C

α

0 Dt ;[z(t ),1]U

(x j , tn ) − dUxx (x j , tn ) − f (x j , tn )

$

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

! =

τz−α

n

9

"

# $ λα,kU jn−k − C0 Dtα;[z(t ),1]U (x j , tn ) + d Uxx (x j , tn ) − δx2U jn

k=0

= O ( τ ) + O ( h2 ), 2 z



where we have used (2.18) and (2.20) in the last equality. The proof is complete.

{unj }

be the solution obtained from the numerical scheme (2.21) and let Let convergence of the numerical scheme (2.21) is established in the next result.

enj

= U ( x j , tn )

− unj

be the error at (xj , tn ). The

Theorem 4.2 (Convergence). Suppose that U (x, Z −1 (z )) = U¯ (x, z ) ∈ C 4,4 ([0, 1] × [z(0 ), z(1)]) and 0, 1, 2, 3, 4. Then,

τz

i

en ∞ = O(τz2 + h2 ),

1≤i≤N

∂ k U¯ | = 0, ∂ zk z=z(0 )

k=

(4.7)

n=1

(0 ) where h is the spatial step size and τz = z(1 )−z . N

Proof. From the numerical scheme (2.21) and the system (4.5), obviously we have

⎧ n ⎪ ⎪ λα,k enj −k = aδx2 enj + Tjn , ⎨τz−α k=0

⎪ e0 = 0, ⎪ ⎩ nj

0≤ j≤M enM = 0,

e0 = 0,

1≤ j ≤M−1 (4.8)

1 ≤ n ≤ N.

Noting that en ∈ Uh , we apply Theorem 3.1 to the system (4.8) to get, for 1 ≤ i ≤ N,

τz

i

%

e ∞ ≤

[z ( 1 ) − z ( 0 ) ]

n

n=1

& [z ( 1 ) − z ( 0 ) ]

≤

& ≤

i τz

24a2

T n 2

n=1

iτz max T n 2 24a2 1≤n≤N

[z ( 1 ) − z ( 0 ) ] max T n 2 , 1≤n≤N 24a2 2

(4.9)

where the fact iτz ≤ z(1 ) − z(0 ) has been used in the last inequality. Next, using (4.6) of Lemma 4.1 and (M − 1 )h < Mh = 1, it is found that

T n 2 = h

M−1

 n 2 Tj

#

= h(M − 1 ) O(τz2 + h2 )

$2

#

$2

= O(τz2 + h2 ) .

(4.10)

j=1

Now, substituting (4.10) into (4.9) gives (4.7) immediately.



Remark 4.1. When i = N, the relation (4.7) establishes the convergence of the numerical scheme (2.21) in ‘average norm’ as N 1 n e ∞ = O(τz2 + h2 ). N

(4.11)

n=1

(0 ) This is clear since τz = z(1 )−z . N

Remark 4.2. Note that the relations (4.7) and (4.11) indicate second order convergence O(τz2 ) with respect to the uniform mesh of [z(0), z(1)]. To investigate the convergence with respect to the mesh of the time interval [0,1], let us denote τ = max1≤n≤N |tn − tn−1 |, the maximum temporal step size of [0,1]. Suppose 0 ≤ z (t) ≤ C for t ∈ [0, 1]. Then, applying mean value theorem gives

τz = z(tn ) − z(tn−1 ) = z (ξ )(tn − tn−1 ) for some ξ ∈ (tn−1 , tn ), and so

τz ≤ C |tn − tn−1 | ≤ C τ , which, together with (4.7) and (4.11), yields

τz

i n=1

e n ∞ = O ( τ 2 + h2 ),

1≤i≤N

(4.12)

10

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

and N 1 n e ∞ = O ( τ 2 + h2 ) N

(4.13)

n=1

respectively. In particular, the relation (4.13) establishes the convergence of the numerical scheme (2.21) in ‘average norm’ with respect to the maximum step size τ of the time interval [0,1]. Remark 4.3. (a) We claim that the temporal convergence order of 2 obtained in (4.7) (or (4.11)) improves the work of [24] where the temporal convergence order is observed to be (2 − α ). Even though our method is based on non-uniform mesh of the time domain while the method in [24] is based on uniform mesh, the improvement can be demonstrated in the next section (Example 5.3) where we compute the maximum absolute error at t = tN . (b) In [24], the convergence of the method is established by using the Lax-Richtmyer Theorem which states that if the numerical scheme is well-posed and consistent, then the scheme is convergent if and only if it is stable. The convergence order in [24] is not proved explicitly but is observed numerically to be (2 − α ). In comparison, the proof of the convergence of our method is more straightforward and the convergence order is explicitly given in (4.7) (or (4.11)). −1 4,4 Remark 4.4.  In view of Remark 2.4, the conditions in Theorem 4.2, namely, U (x, Z (z )) = U˜ (x, z ) ∈ C ([0, 1] × [z(0 ), z(1 )] ) k  and ∂ k U˜  = 0, k = 0, 1, 2, 3, 4, can be replaced by the weaker conditions U(x, · ) ∈ C(4) [0, 1] and U˜e (·, z ) ∈ L1 (IR ) ∩

∂z

C α +2 (IR ).

z=z (0 )

5. Simulation In this section, we shall present four examples. The main objective is to demonstrate the efficiency of the proposed numerical scheme (2.21), and to confirm the theoretical temporal convergence of O(τz2 ) in (4.11) – which is immediate from Theorem 4.2. Moreover, some numerical experiments are conducted when the conditions of Theorem 4.2 are not entirely satisfied, just to see the influence of the conditions on the convergence of the numerical scheme. Finally, the effects of different scale functions z(t), weight functions w(t ) and diffusion constants a on the numerical solutions are also investigated. We implement all the numerical experiments using Matlab R2014b on a MacBook Pro-with macOS Mojave system (Version 10.14.5), 2.6 GHz Intel Core i5 Processor, 8GB 1600 MHz DDR3 Memory and Intel Iris 1536MB Graphics. Let {U nj } be the exact solution of the sub-diffusion problem and {unj } be the numerical solution obtained from the scheme (2.21) or some other schemes in the literature. In the following examples, we shall compute the average norm and maximum norm errors as follows

Eav (h, N ) = τz

N n=1

max

1≤ j≤M−1

|U jn − unj |

and

E∞ (h, N ) = max

1≤ j≤M−1

|U jN − uNj |.

The temporal convergence orders in terms of average norm and maximum norm are computed respectively by



Tav = log2

Eav (h, N ) Eav (h, 2N )





and

T∞ = log2



E∞ (h, N ) . E∞ (h, 2N )

(5.1)

Note that Theorem 4.2 (or (4.11)) shows the convergence in terms of average norm, nonetheless it would also be interesting to see the convergence in terms of maximum norm. In some examples where the exact solution of the sub-diffusion problem exists but is unknown or extremely difficult to get, we shall compute the ‘exact’ solution by using (2.21) with sufficiently small spatial step size h and sufficiently large number of time intervals N. In fact, we shall use h = 10100 and N = 10 0 0 0 to pre-compute the ‘exact’ solution in Examples 5.2 and 5.3. Example 5.1. Consider the sub-diffusion problem

⎧   ( γβ + 1 ) γ −βα γ ⎪ ⎪ ⎨C0 Dtα;[z(t ),1] u(x, t ) = uxx (x, t ) + ex t − t , ( γβ + 1 − α ) ⎪ u(x, 0 ) = 0, x ∈ ( 0, 1 ) ⎪ ⎩ u ( 0, t ) = t γ , u(1, t ) = et γ , t ∈ ( 0, 1 )

(x, t ) ∈ (0, 1 ) × (0, 1 ) (5.2)

where z(t ) = t β , γ ≥ βα > 0 and 0 < α < 1. The exact solution is U (x, t ) = exp(x )t γ . Case 5.1. Let γ = 6 + α and β = 0.5, 1.5. It is easy to verify that the conditions in Theorem 4.2 are satisfied. Applying the scheme (2.21) with fixed h = 10100 , in Table 1 we present the error Eav (h, N ), the temporal convergence order Tav and the total CPU time taken T(sec) (in seconds) for different α ’s. Clearly, the proposed numerical scheme (2.21) gives second order accuracy which is consistent with the theoretical result of O(τz2 ).

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

11

Table 1 (Example 5.1) Temporal convergence order Tav . z(t) t

0.5

t1.5

N

α = 0.2

Tav

T(sec)

α = 0.5

Tav

T(sec)

α = 0.8

Tav

T(sec)

40 80 160 40 80 160

1.08E-04 2.70E-05 6.73E-06 2.91E-05 7.28E-06 1.82E-06

2.01 2.00 2.00 2.00

0.15 0.16

5.89E-04 1.49E-04 3.74E-05 1.24E-04 3.10E-05 7.76E-06

1.99 1.99 2.00 2.00

0.18 0.18

1.83E-03 4.67E-04 1.18E-04 3.29E-04 8.27E-05 2.07E-05

1.97 1.98 1.99 2.00

0.14 0.16

Table 2 (Example 5.1) Tav when conditions of Theorem 4.2 are not satisfied.

γ

N

α = 0.2

Tav

T(sec)

α = 0.5

Tav

T(sec)

α = 0.8

Tav

T(sec)

α+1

40 80 160 40 80 160

6.64E-05 2.33E-05 8.08E-06 6.21E-06 1.62E-06 4.25E-07

1.51 1.53 1.94 1.93

0.18 0.12

1.69E-04 6.15E-05 2.22E-05 2.17E-05 5.56E-06 1.43E-06

1.46 1.47 1.96 1.96

0.17 0.15

3.49E-04 1.34E-04 5.15E-05 4.57E-05 1.16E-05 2.94E-06

1.37 1.38 1.98 1.98

0.15 0.14

α+2

Table 3 (Example 5.2) Temporal convergence order T∞ . z(t)

N

α = 0.2

T∞

T(sec)

α = 0.5

T∞

T(sec)

α = 0.8

T∞

T(sec)

t0.5

50 100 200 50 100 200 50 100 200 50 100 200

6.42E-06 1.62E-06 4.06E-07 7.98E-07 1.99E-07 4.98E-08 6.30E-08 1.66E-08 4.34E-09 2.43E-07 6.03E-08 1.50E-08

1.99 1.99 2.00 2.00 1.92 1.94 2.01 2.01

0.48 0.45 0.44 0.43

2.31E-05 5.83E-06 1.46E-06 1.56E-06 3.88E-07 9.68E-08 2.78E-07 6.75E-08 1.65E-08 7.09E-07 1.75E-07 4.34E-08

1.99 1.99 2.00 2.00 2.04 2.03 2.02 2.01

0.49 0.43 0.42 0.45

5.05E-05 1.27E-05 3.19E-06 1.14E-06 2.82E-07 7.03E-08 9.40E-07 2.32E-07 5.78E-08 3.46E-07 8.27E-08 2.02E-08

1.99 1.99 2.01 2.00 2.02 2.01 2.06 2.03

0.45 0.43 0.44 0.48

t

t1.5

e2t − 1

Case 5.2. We shall investigate some cases where the exact solution U (x, t ) = exp(x )t γ may not be as smooth as assumed in Theorem 4.2 which is possible in engineering problems. Let β = 1.5 and fix h = 10100 in the implementation of the scheme (2.21). Consider the cases (i) γ = α + 1, and (ii) γ = α + 2. (i) When γ = α + 1, we have U¯ (·, z ) ∈ C ([z(0 ), z(1 )] ) for α ∈ (0, 0.5), and U¯ (·, z ) ∈ C (1 ) ([z(0 ), z(1 )] ) for α ∈ [0.5, 1). (ii) When γ = α + 2, we have U¯ (·, z ) ∈ C (1 ) ([z(0 ), z(1 )] ) for α ∈ (0, 1). Obviously, both cases somewhat violate the conditions of Theorem 4.2. The errors, temporal convergence orders and total CPU times taken of these two cases are displayed in Table 2. It is observed from Table 2 that the temporal convergence order is indeed lower when γ = α + 1; but the temporal convergence order is still close to 2 (as in Theorem 4.2) when γ = α + 2. The next example is modified from examples in [12,13] by replacing the classical Caputo fractional derivative by the generalized Caputo fractional derivative. Example 5.2. Consider the sub-diffusion problem

⎧ ' ⎪ ⎨C0 Dtα;[z(t ),1] u(x, t ) = uxx (x, t ) + u(x, 0 ) = 0, ⎪ ⎩ u ( 0, t ) = 0,

x ∈ ( 0, 1 ) u ( 1, t ) = 0,

(

2 + 4π 2 t 2 sin 2π x, (3 − α )t 2−α

(x, t ) ∈ (0, 1 ) × (0, 1 ) (5.3)

t ∈ ( 0, 1 )

where z(t ) = t 0.5 , t, t 1.5 , exp(2t ) − 1 and 0 < α < 1. Here, we fix h = 10100 in the scheme (2.21) and let N vary. The maximum absolute error at t = tN ,E∞ (h, N), the temporal convergence order T∞ and the total CPU time taken T(sec) are presented in Table 3; while the average norm error Eav (h, N ), the temporal convergence order Tav and the total CPU time taken T(sec) are displayed in Table 4. It is obvious from the tables that the numerical scheme (2.21) can achieve second order accuracy in temporal dimension for various fractional orders and scale functions.

12

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991 Table 4 (Example 5.2) Temporal convergence order Tav . z(t) t

0.5

t

t1.5

e2t − 1

N

α = 0.2

Tav

T(sec)

α = 0.5

Tav

T(sec)

α = 0.8

Tav

T(sec)

50 100 200 50 100 200 50 100 200 50 100 200

2.33E-06 5.83E-07 1.46E-07 9.86E-07 2.47E-07 6.17E-08 2.71E-07 7.25E-08 1.94E-08 1.05E-06 2.68E-07 6.76E-08

2.00 2.00 2.00 2.00 1.90 1.90 1.97 1.99

0.50 0.40 0.40 0.38

9.30E-06 2.33E-06 5.83E-07 2.97E-06 7.47E-07 1.87E-07 1.95E-06 5.71E-07 1.59E-07 6.58E-06 1.68E-06 4.27E-07

2.00 2.00 1.99 1.99 1.77 1.85 1.97 1.98

0.45 0.40 0.37 0.45

2.26E-05 5.67E-06 1.42E-06 4.92E-06 1.25E-06 3.15E-07 1.30E-05 3.89E-06 1.09E-06 2.68E-05 7.01E-06 1.77E-06

1.99 2.00 1.98 1.98 1.74 1.84 1.94 1.98

0.39 0.43 0.41 0.45

Table 5 (Example 5.3) Temporal convergence orders T∞ and Tav .

w(t ) ≡ 1, z(t ) = t

w(t ) ≡ 1, z(t ) = t 0.5

w(t ) = et , z(t ) = t

w(t ) = et , z(t ) = t 0.5

N

E∞ (h, N)

T∞

T(sec)

Eav (h, N )

Tav

T(sec)

50 100 200 50 100 200 50 100 200 50 100 200

1.06E-06 2.64E-07 6.58E-08 4.82E-05 1.22E-05 3.05E-06 1.83E-05 4.64E-06 1.17E-06 9.92E-05 2.60E-05 6.67E-06

2.01 2.00 1.99 1.99 1.98 1.99 1.93 1.96

0.49 0.55 0.48 0.48

5.11E-06 1.29E-06 3.27E-07 2.07E-05 5.19E-06 1.30E-06 5.35E-05 1.45E-05 3.77E-06 3.29E-05 8.32E-06 2.09E-06

1.98 1.98 1.99 2.00 1.88 1.95 1.98 1.99

0.46 0.50 0.43 0.41

Our third example is from [24] where the authors only conducted numerical experiment for z(t ) = t and w(t ) ≡ 1 (classical Caputo fractional derivative). Here, we shall illustrate the improvement in the temporal convergence order of the proposed scheme (2.21) for different scale functions and weight functions, including the special case of z(t ) = t and w(t ) ≡ 1. Example 5.3 [24]. Consider the sub-diffusion problem

⎧ ⎪ ⎨C0 Dt0;[.85 z(t ),w(t )] u (x, t ) = uxx (x, t ) +

2

(x2 − x )t 1.15 + π 2 sin(π x ) − 2t 2 , (2.15 ) u(x, 0 ) = sin(π x ), x ∈ ( 0, 1 ) ⎪ ⎩ u ( 0, t ) = 0, u ( 1, t ) = 0, t ∈ ( 0, 1 ).

(x, t ) ∈ (0, 1 ) × (0, 1 ) (5.4)

When z(t ) = t and w(t ) ≡ 1, the exact solution is U (x, t ) = sin(π x ) + x(x − 1 )t 2 . We fix h = 10100 in the implementation of the scheme (2.21). In Table 5, we present the errors E∞ (h, N ), Eav (h, N ), the temporal convergence orders T∞ , Tav , as well as the total CPU time taken T(sec) for four sets of w(t ) and z(t). Note that the numerical scheme (2.21) is applied directly when w(t ) ≡ 1. If w(t ) ≡ 1, we shall first convert the original problem to the equivalent form (1.5) and then use the scheme (2.21) to compute the numerical solution of (1.5). The numerical solution of the original problem is then obtained via a relation that relates the solutions of (1.3) and (1.5) (see Introduction). From the experimental results, it is clear that the temporal convergence order of 2 is achieved in all cases, which is consistent with the theoretical result. Next, we shall compare the proposed numerical scheme (2.21) with the method in [24]. With h = 10100 fixed, we present the maximum absolute error at t = tN ,E∞ (h, N), the convergence order T∞ and the total CPU time taken T(sec) in Table 6. It is observed that the proposed scheme (2.21) achieves second order accuracy in the temporal dimension as well as smaller errors for all the cases except when w(t ) ≡ 1, z(t ) = t 0.5 . However, since our scheme achieves T∞ ≈ 2 while T∞ ≈ 2 − α = 1.15 in [24], our method will eventually give a smaller error as N becomes larger. The last example is first used in [24] to illustrate the effects of the scale function z(t) and weight function w(t ) on the numerical solution obtained by the method in [24]. Example 5.4 [24]. Consider the sub-diffusion problem

⎧ ⎪ ⎨C0 Dtα;[z(t ),w(t )] u(x, t ) = duxx (x, t ) + (x − 1 )(t − 1 ) , (x, t ) ∈ (0, 1 ) × (0, 1 ) 2 + sin(xt ) 3 x ∈ ( 0, 1 ) ⎪ ⎩u(x, 0 ) = x(1 − x ) exp(sin(3.5π x )), u ( 0, t ) = 0, u ( 1, t ) = 0, t ∈ ( 0, 1 )

(5.5)

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991

13

Fig. 1. (Example 5.5) Numerical solution of (a) w(t ) = exp(t ), z(t ) = t, a = 1; (b) w(t ) = exp(t ), z(t ) = t 5 , a = 1; (c) w(t ) = exp(5t ), z(t ) = t, a = 1; (d) w(t ) = exp(t ), z(t ) = t, a = 0.5. Table 6 (Example 5.3) Comparison between scheme (2.21) and scheme in [24].

w(t ) ≡ 1, z(t ) = t

w(t ) ≡ 1, z(t ) = t 0.5

w(t ) = et , z(t ) = t

w(t ) = et , z(t ) = t 0.5

N

Scheme (2.21)

T∞

T(sec)

Scheme in [24]

T∞

T(sec)

50 100 200 50 100 200 50 100 200 50 100 200

1.06E-06 2.64E-07 6.58E-08 4.82E-05 1.22E-05 3.05E-06 1.83E-05 4.64E-06 1.17E-06 9.92E-05 2.60E-05 6.67E-06

2.01 2.00 1.99 1.99 1.98 1.99 1.93 1.96

0.49 0.55 0.48 0.48

2.25E-04 9.72E-05 3.98E-05 7.21E-06 3.14E-06 1.29E-06 2.56E-04 1.12E-04 4.62E-05 2.68E-04 1.16E-04 4.76E-05

1.21 1.29 1.20 1.29 1.19 1.28 1.20 1.29

0.75 0.87 0.74 0.74

where d > 0 is the diffusion constant and 0 < α < 1. In this example, we shall present illustrative figures to investigate the effects of the scale function z(t), weight function w(t ) and diffusion constant d on the numerical solution computed from the scheme (2.21). To have proper comparisons, the 1 parameters α = 0.97, h = 100 and N = 10 0 0 are fixed for all cases. In Fig. 1, we plot the numerical solutions for the following four cases: (a) (b) (c) (d)

Base case: w(t ) = exp(t ), z(t ) = t, a = 1 Change z(t): w(t ) = exp(t ), z(t ) = t 5 , a = 1 Change w(t ): w(t ) = exp(5t ), z(t ) = t, a = 1 Change d: w(t ) = exp(t ), z(t ) = t, a = 0.5

14

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991 1

3

1

0.8

0.00 0.01

5

0.8

0.01

0.01

5

0.7 0 .0

0.6

0.0

15

0.

02

0.5

0.01

25

0.0

x

1

0.6

25

0.0

0.02

0.5

0.003

0.01 0.0 1

0.9

3

03 0.0

0.01 0.01 5 0.02 0.02 5

0.7

x

0.003

0.003

00 0.

0.9

0.4

0.01 5

0.3

0.3

2

0.0

0.1 0 0

0.003 0.1 0.2

5

0.01

0.2 0.01

0.015

0.1

0.01

0.003 0.3

0.4

0.5

0.6

t

0 0.7

0.8

0.9

1

02 0. 15 0 . 00.01 0.003 0.1 25

25 0.0 0.02

0.01

0.0

0.2

0

0.0 03

0.4

0.003

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

t

Fig. 2. (Example 5.5) Contour of (a) (left) and contour of (c) (right) in Fig. 1 with level list [0.0 0 03, 0.01, 0.015, 0.02, 0.025].

The following observations are noted. (i) In Fig. 1(b), it is obvious that the scale function z(t ) = t 5 stretches the time domain at the beginning which results in a coarser mesh compared to Fig. 1(a). (ii) Fig. 1(c) depicts a larger weight function w(t ) = exp(5t ) compared to Fig. 1(a). Typically, the solution tends to be smaller with larger weights, although this is hard to see from Fig. 1(c). Therefore, for cases (a) and (c), we plot the contours with a given level list [0.0 0 03 0.01 0.015 0.02 0.025] in Fig. 2. From the contours in Fig. 2, it is easily seen that the numerical solutions of cases (a) and (c) are very different. In fact, since the weight function w(t ) = exp(5t ) is larger than w(t ) = exp(t ), the contour of case (c) is much closer to t = 0 than case (a) for the same level. This indicates that the solution of case (c) tends to be smaller at the same mesh point than that of case (a). (iii) Fig. 1(d) depicts a smaller diffusion constant a = 0.5 compared to Fig. 1(a). Indeed, it is observed from Fig. 1(d) that the diffusion process of case (d) has become slower than that of case (a) where a = 1 is used. It is noted in [3] that the scale function z(t) can stretch or contract the domain, the weight function w(t ) can weigh differently at different time points, and the diffusion constant d measures the rate of the diffusion process. From the above observations (i)–(iii), we conclude that the numerical solutions computed from the scheme (2.21) do illustrate correctly the effects of the scale function z(t), the weight function w(t ) and the diffusion constant d. This further illustrates the efficacy of the numerical scheme (2.21). 6. Conclusion In this paper, we consider the numerical treatment of a generalized fractional sub-diffusion problem. A new approximation for the generalized fractional derivative, namely the generalized weighted shifted Grünwald-Letnikov (gWSGL) approximation has been developed. Using the gWSGL approximation, we construct an efficient numerical scheme for the generalized fractional sub-diffusion problem. The solvability, stability and convergence have been established rigorously. It is shown that second order accuracy is achieved in the temporal dimension and this is the best result to date. To demonstrate the efficiency and accuracy of the proposed method, we have carried out simulation for four examples. The numerical experiments indicate that the convergence of our method is consistent with the theoretical result. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

Agarwal RP. Difference equations and inequalities. Marcel Dekker Inc. New York; 1992. Agrawal OP. A general solution for a fourth-order fractional diffusionwave equation defined in a bounded domain. Comp Struct 2001;79:1497–501. Agrawal OP. Some generalized fractional calculus operators and their applications in integral equations. Fract Calc Appl Anal 2012;15:700–11. Agrawal OP. Generalized variational problems and Euler-Lagrange equations. Comput Math Appl 2010;59:1852–64. Gao G, Sun ZZ. Two alternating direction implicit difference schemes for two-dimensional distributed-order fractional diffusion equations. J Sci Comput 2016;66:1281–312. Hu X, Zhang L. On finite difference methods for fourth-order fractional diffusion-wave and sub-diffusion systems. Appl Math Comput 2012;218:5019–34. Ji CC, Sun ZZ. A high-order compact finite difference scheme for the fractional sub-diffusion equation. J Sci Comput 2015;64:959–85. Ji CC, Sun ZZ. The high-order compact numerical algorithms for the two-dimensional fractional sub-diffusion equation. Appl Math Comput 2015;269:775–91. Katugampola UN. A new approach to generalized fractional derivatives. Bull Math Anal Appl 2014;6:1–15. Kim M-H, Ri G-C, O H-C. Operational method for solving multi-term fractional differential equations with the generalized fractional derivatives. Fract Calc Appl Anal 2014;17:79–95. Kiryakova V. A brief story about the operators of the generalized fractional calculus. Fract Calc Appl Anal 2008;11:203–20. Li X, Wong PJY. A higher order non-polynomial spline method for fractional sub-diffusion problems. J Comput Phys 2017;328:46–65.

X. Li and P.J.Y. Wong / Commun Nonlinear Sci Numer Simulat 82 (2020) 104991 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]

15

Li X, Wong PJY. Non-polynomial spline approach in two-dimensional fractional sub-diffusion problems. Appl Math Comput 2019;357:222–42. Lubich C. Discretized fractional calculus. SIAM J Math Anal 1986;17:704–19. Mainardi F, Paradisi P. Fractional diffusive waves. J Comput Acoust 2001;9:1417–36. Metzler R, Klafter J. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys Rep 20 0 0;339:1–77. Oldham KB, Spanier J. The fractional calculus: theory and applications of differentiation and integration to arbitrary order. Academic Press, New York; 1974. Povstenko Y. Linear fractional diffusion-Wave equation for scientists and engineers. Birkhäuser, New York; 2015. Tian W, Zhou H, Deng W. A class of second order difference approximations for solving space fractional diffusion equations. Math Comp 2015;84:1703–27. Wang Z, Vong S. Compact difference schemes for the modified anomalous fractional sub-diffusion equation and the fractional diffusion-wave equation. J Comput Phys 2014;277:1–15. Xu Y, Agrawal OP. Models and numerical schemes for generalized van der Pol equations. Commun Nonlinear Sci Numer Simul 2013;18:3575–89. Xu Y, Agrawal OP. Numerical solutions and analysis of diffusion for new generalized fractional Burgers equation. Fract Calc Appl Anal 2013;16:709–36. Xu Y, Agrawal OP. Models and numerical solutions of generalized oscillator equations. J Vib Acoust 2014;136:050903–050903–7. Xu Y, He Z, Agrawal OP. Numerical and analytical solutions of new generalized fractional diffusion equation. Comput Math Appl 2013;66:2019–29. Xu Y, He Z, Xu Q. Numerical solutions of fractional advection-diffusion equations with a kind of new generalized fractional derivative. Int J Comput Math 2014;91:588–600. Zhang Y, Sun ZZ. Alternating direction implicit schemes for the two-dimensional fractional sub-diffusion equation. J Comput Phys 2011;230:8713–28. Zorich VA. Mathematical analysis II. Springer, Berlin; 2004.