An efficient split-step method for distributed-order space-fractional reaction-diffusion equations with time-dependent boundary conditions

Applied Numerical Mathematics 147 (2020) 142–160

Contents lists available at ScienceDirect

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

An efficient split-step method for distributed-order space-fractional reaction-diffusion equations with time-dependent boundary conditions Kamran Kazmi a,∗ , Abdul Q.M. Khaliq b a

Department of Mathematics, University of Wisconsin Oshkosh, Oshkosh, WI 54901, USA Department of Mathematical Sciences and Center for Computational Sciences, Middle Tennessee State University, Murfreesboro, TN 37132, USA b

a r t i c l e

i n f o

Article history: Received 11 June 2019 Received in revised form 17 August 2019 Accepted 20 August 2019 Available online 27 August 2019 Keywords: Fractional Laplacian Distributed-order space-fractional reaction diffusion equations Time-dependent boundary conditions Matrix transfer technique Predictor-corrector method

a b s t r a c t A split-step predictor-corrector method, based on Adams-Moulton formula, is presented for the solution of nonlinear distributed-order space-fractional reaction-diffusion equations with time-dependent boundary conditions. The distributed-order space-fractional equation is first discretized to a multi-term space-fractional equation by applying a quadrature rule. Multi-term space-fractional equation is spatially discretized by matrix transfer technique and a linearly implicit split-step predictor-corrector method is applied for time-stepping to avoid solving nonlinear systems at each time step. The method is shown to be unconditionally stable and second order convergent. Numerical experiments are performed to confirm the stability and second order convergence of the method. The split-step predictor-corrector method is also compared with a second order IMEX predictor-corrector method. The IMEX predictor-corrector method is found to incur oscillatory behavior when implemented on problems with nonsmooth initial conditions or with mismatch between the initial and boundary conditions. Our method, which is an L-stable method, produces reliable and oscillation free results for such problems. © 2019 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction Fractional partial differential equations with time or space fractional derivatives has become very popular in the recent years. They are found to model various physical processes better than the corresponding integer order differential equations due to the nonlocal nature of fractional derivatives. Many problems in physical sciences, engineering and finance can be modeled by fractional partial differential equations. Some examples are anomalous transport [35], anomalous diffusion [34, 40], highly heterogeneous aquifer [3], materials and mechanics [6,43], wave propagation and turbulence [36], and financial derivatives [12,53]. In general, a single or multi-term fractional equation is used to model these phenomenons. Several numerical methods have been developed for solving problems involving single term fractional operators. See [8,13,49,21] for space-fractional differential equations and [15,19,31,42,54,56] for time-fractional differential equations. There are many anomalous diffusion processes that can not be characterized by single or multi-term fractional equations as their diffusion exponents vary with time. Such processes are more suitably described by distributed order fractional equations. In [14],

*

Corresponding author. E-mail address: [email protected] (K. Kazmi).

https://doi.org/10.1016/j.apnum.2019.08.019 0168-9274/© 2019 IMACS. Published by Elsevier B.V. All rights reserved.

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

143

Chechkin et al. used distributed order time and space fractional diffusion equations to analyze retarding subdiffusion and acceleration superdiffusion. Lévy flights in nonhomogeneous media are studied using distributed order fractional equation in [48]. For further applications of distributed order fractional equations, see [5,10,11,23,32,38,47]. Several authors have presented numerical methods for solving distributed-order fractional differential equations. A general form of distributed-order time-fractional ordinary differential equation is solved in [24,33]. Numerical methods for solving distributed-order time-fractional diffusion equations are presented in [27,37,55]. Numerical methods for distributedorder space-fractional diffusion equations are provided in [1,7,18,30]. For numerical methods for solving distributed-order space-fractional advection-diffusion equations, see [29,52,57]. Most of the articles in literature have reported numerical results for linear distributed-order space or time fractional partial differential equations under homogeneous boundary conditions. In this paper, we study a nonlinear distributed-order space-fractional differential equation under nonhomogeneous time-dependent boundary conditions. In particular, we present a numerical method for the solution of a distributed-order space-fractional reaction-diffusion equation

2 u t ( x, t ) = −

α

κ (α )(−) 2 u (x, t )dα + f (x, t , u (x, t )), (x, t ) ∈  × (0, T ],

(1.1)

1

with the initial condition

u (x, 0) = u 0 (x),

x ∈ ,

(1.2)

and nonhomogeneous time-dependent boundary conditions of one of the following types:

where

u (x, t ) = h(t ),

(x, t ) ∈ ∂  × (0, T ],

(1.3a)

∇ u (x, t ) · n = h(t ),

(x, t ) ∈ ∂  × (0, T ],

(1.3b)

β u (x, t ) + ∇ u (x, t ) · n = h(t ),

(x, t ) ∈ ∂  × (0, T ],

κ (α ) is a nonnegative weight function which satisfies κ (α ) ≡ 0 and

(1.3c)

2 1

κ (α )dα < ∞, (−)α /2 is the d-dimensional

space-fractional Laplacian of order α for 1 < α ≤ 2, n is the outward unit normal, and  ⊂ Rd (d = 1, 2) is a bounded domain. We assume that f , h, and u 0 are sufficiently smooth functions. Also, the space-fractional Laplacian (−)α /2 is expressed using the spectral representation (see [21,22]). For other definitions of the space-fractional Laplacian, see [39,44, 45]. Definition 1. [22] Let {ϕn } be the complete set of the orthonormal eigenfunctions corresponding to the eigenvalues λn2 of the Laplacian (−) on a bounded domain  with homogeneous boundary conditions. Then α

(−) φ = 2

⎧ ⎪ ⎨

α = 2m, m = 0, 1, 2, · · · ,

(−)m φ, α −m

(−) 2 (−) φ, m − 1 < α2 < m, m = 1, 2, · · · , ⎪ ⎩ ∞ α α < 0. n=1 λn φ, ϕn ϕn , m

α

The action of (−) 2 can be extended to functions which do not satisfy homogeneous boundary conditions through the formula (1.4). The formula is stated for one dimensional case but it can be extended to two or three dimensions. Proposition 2. [22] Let ϕn be an eigenfunction corresponding to the eigenvalue λn2 of the one dimensional Laplacian − on the interval a ≤ x ≤ b and φ be a sufficiently smooth function. Then α

α

ϕn , (−) 2 φ − (−) 2 ϕn , φ = λnα −2 [φ(b)ϕn (b) − φ(a)ϕn (a) − φ (b)ϕn (b) + φ (a)ϕn (a)].

(1.4)

The purpose of this paper is to extend the approach in [25] to the distributed-order space-fractional reaction diffusion problems in one and two space dimensions. The integral term of the distributed-order differential equation brings extra challenges for stability and convergence and makes the full dense matrix highly ill conditioned, especially in multi space dimensions which needs careful adaptation of reliable algorithms. Furthermore, the computational cost of the numerical solution for the distributed-order fractional operator is much higher than the single order fractional operator, requiring efficient and accurate algorithms. As far as we know, the literature is lacking such algorithms for distributed-order spacefractional partial differential equations with nonhomogeneous time dependent boundary conditions. We formulate and analyze a second order L-stable method for the solution of the problem (1.1)-(1.3). First, the equation (1.1) is discretized to a multi-term space-fractional equation by applying the midpoint quadrature rule. Multi-term space-fractional equation is then spatially discretized by matrix transfer technique (MTT). Finally, a linearly implicit splitstep predictor-corrector method, based on Adams-Moulton formula, is applied for time-stepping to avoid solving nonlinear

144

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

systems at each time step. The unconditional stability and second order convergence of the method are theoretically established. The split-step predictor-corrector method is also compared with a second order IMEX predictor-corrector method. The IMEX predictor-corrector method produces unwanted oscillations in numerical solutions when applied to problems with nonsmooth initial conditions or with mismatch between the initial and boundary conditions. The split-step predictorcorrector method, which is an L-stable method, produces reliable and oscillation free solutions for such problems. For IMEX numerical schemes for solving partial differential equations, see [4,9,20,28,41,51]. MTT was introduced by Ilic´ et al. [21] for the space-fractional diffusion equation with homogeneous boundary conditions and later extended for nonhomogeneous boundary conditions in [22]. The basic idea of MTT is to approximate the fractional Laplacian (−)α /2 by the discrete matrix representation A α /2 , where the matrix A is the discrete representation of the standard Laplacian operator. This approach is later followed by many researchers for solving space-fractional diffusion equations. A fourth order method based on MTT is presented by Ding and Zhang [17]. A scheme based on rational approximation to the fractional Laplacian operator is developed by Aceto and Novati [2]. A space fractional telegraph equation is solved by using diagonal Padé approximants by Chen et al. [16] and a space-fractional reaction-diffusion equation is solved by (1, 1)-Padé and (0, 2)-Padé approximants by Khaliq et al. [26]. A finite difference numerical method based on MTT for solving a distributed-order space-fractional reaction-diffusion equation is presented in [7]. The remainder of the paper is organized as follows. In Section 2, we derive the analytical solution of a linear distributedorder space-fractional reaction-diffusion equation. The development of the split-step predictor-corrector method is discussed in Section 3. Section 4 provides the stability and convergence results for the numerical method. The numerical experiments are performed in Section 5 and the concluding remarks are presented in Section 6. 2. Analytic solution In this section, we derive the analytic solution of a linear distributed-order space-fractional reaction diffusion equation. Since nonlinear distributed-order space-fractional equations often can not be analytically solved, the analytic solution of this linear test problem will be compared with the numerical solution computed by the predictor-corrector method to demonstrate the effectiveness of our method. We consider a one-dimensional problem of the form

∂u =− ∂t

2



∂2 κ (α ) − 2 ∂x

α /2 u dα + f (x) g (t ),

0 < x < 1, t > 0,

(2.1)

1

u (0, t ) = h1 (t ),

t > 0,

(2.2)

u (1, t ) = h2 (t ),

t > 0,

(2.3)

u (x, 0) = u 0 (x),

0 < x < 1,

where 1 < α ≤ 2 and

κ (α ) is a nonnegative weight function, which satisfies κ (α ) ≡ 0 and

2 1

(2.4)

κ (α )dα < ∞. The operator √

2 − ∂∂x2 with homogeneous Dirichlet boundary conditions is self-adjoint and has eigenfunctions ϕn = 2 sin(nπ x) and correα 2 2 sponding eigenvalues λn2 = (nπ )2 for n = 1, 2, · · · . For a fixed α , the fractional Laplacian (− ∂∂x2 ) 2 as a function of − ∂∂x2 is

also self-adjoint. Using equation (1.4), we can extend its definition to a function φ that do not vanish at x = 0 and x = 1 as follows:

ϕn , (−

 ∂2 α ) 2 φ = λnα ϕn , φ + λnα −2 φ(1)ϕn (1) − φ(0)ϕn (0) . ∂ x2

(2.5)

Applying the Laplace transform to equations (2.1)–(2.3), we obtain



2 su¯ (x, s) − u 0 (x) = −

κ (α ) −

∂2 ∂ x2

α /2 u¯ (x, s) dα + f (x) g¯ (s),

(2.6)

1

u¯ (0, s) = h¯ 1 (s),

(2.7)

u¯ (1, s) = h¯ 2 (s), (2.8)  ∞ −st where u¯ (x, s) = 0 e u (x, t ) dt is the Laplace transform of u (x, t ). Similarly g¯ (s), h¯ 1 (s), and h¯ 2 (s) are the Laplace trans1 forms of g (t ), h1 (t ), and h2 (t ) respectively. With u , v  = 0 u (x) v (x) dx as the inner product, the inner product of equation (2.6) with the eigenfunction

ϕn results in 2

s ϕn , u¯  − ϕn , u 0  = −



κ (α ) ϕn , − 1

∂2 ∂ x2

α /2 u¯  dα + g¯ (s) ϕn , f .

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

145

Using equation (2.5), we obtain

2 s ϕn , u¯  − ϕn , u 0  = − ϕn , u¯ 

 2

κ (α )λnα dα − h¯ 2 (s)ϕn (1) − h¯ 1 (s)ϕn (0)

1

Let c¯n (s) = ϕn (x), u¯ (x, s), u 0,n = ϕn (x), u 0 (x), f n = ϕn (x), f (x), and kn = comes

s + kn

kn

−

λn2 (s

κ (α )λnα dα . Then, the above equation be-

λn

or

u 0,n

1

1

 kn ¯ ¯ 1 (s)ϕ (0) + g¯ (s) f n h ( s ) ϕ ( 1 ) − h 2 n n 2

(s + kn )¯cn (s) = u 0,n −

c¯n (s) =

2

κ (α )λnα −2 dα + g¯ (s) ϕn , f .



h¯ 2 (s)ϕn (1) − h¯ 1 (s)ϕn (0) + f n

+ kn )

g¯ (s)

.

s + kn

Now, using the splitting

kn

λn2 (s

+ kn )

=

1

s

−

λn2

λn2 (s

√

+ kn )

, √

ϕn (0) = 2 λn and ϕn (1) = (−1)n 2 λn , we obtain √ √ √ 

u 0,n 2 2 2s g¯ (s) c¯n (s) = + (−1)n+1 h¯ 2 (s) + h¯ 1 (s) + h¯ 2 (s)(−1)n − h¯ 1 (s) + f n . s + kn λn λn λn (s + kn ) s + kn

with the fact that

Thus,

u¯ (x, s) =

∞  n =1

=

∞  n =1

+

c¯n (s)ϕn (x)

√

2u 0,n

s + kn

∞  n =1

sin(nπ x) + xh¯ 2 (s) + (1 − x)h¯ 1 (s)



2s

h¯ 2 (s)(−1)n − h¯ 1 (s) sin(nπ x) +

λn (s + kn )

n =1

√

2

λn

√

n =1

where we used the facts that ∞ 

∞ 

(−1)n+1 sin(nπ x) = x and

∞ 

√

n =1

2

λn

2 g¯ (s) f n

s + kn

sin(nπ x),

(2.9)

sin(nπ x) = 1 − x.

Finally, applying the inverse Laplace transform, we obtain the solution of the problem (2.1)-(2.4) given by

u (x, t ) =

∞ √ 

2u 0,n sin(nπ x)e −kn t + xh2 (t ) + (1 − x)h1 (t )

n =1

+

 ∞  2 sin(nπ x) λn

n =1

+

t





e −kn (t −τ ) (−1)n H 2 (τ ) − H 1 (τ ) dτ

0

t

∞ √ 

2 f n sin(nπ x)

n =1

e −kn (t −τ ) g (τ )dτ ,

(2.10)

0

where H i (t ) = L−1 {sh¯ i (s)} for i

= 1, 2, and g (t ) = L−1 { g¯ (s)}. √ In Example 1 of Section 5, we take κ (α ) = α , u 0 (x) = −(x2 − x − 1), f (x) = 2 sin(π x), g (t ) = et and h1 (t ) = h2 (t ) = 1.

Then

u 0,n = ϕn (x), u 0 (x)

√ =− 2

and



2

(nπ

)3

+

1 nπ



 (−1)n − 1 ,

n ≥ 1,

146

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

f n = ϕn (x), f (x)



=

1 for n = 1, 0 for n = 1.

Also h¯ 1 (s) = h¯ 2 (s) =

s s + kn

and g¯ (s) = s−1 1 . Thus

1 s

h¯ i (s) =

1 s + kn

i = 1, 2,

,

and

g¯ (s) s + kn

= =

1

(s + kn )(s − 1) An s + kn

+

Bn s−1

,

where A n = −1/(kn + 1) and B n = 1/(kn + 1). Substituting these in equation (2.9), we obtain

u¯ (x, s) = −4

∞  [(−1)n − 1] sin(nπ x)

(nπ )3

n =1

s + kn

+

1 s

√

+



2

k1 + 1

1 s−1

−

1



s + k1

sin(π x).

Applying the inverse Laplace transform, we obtain the solution of Example 1 of Section 5, given by

√

u (x, t ) = 1 +

2

k1 + 1

(et − e −k1 t ) sin(π x) − 4

∞  [(−1)n − 1] n =1

(nπ )3

e −kn t sin(nπ x),

(2.11)

.

(2.12)

where

2

α (nπ )α dα = −

kn =

nπ (nπ − 1) − nπ log(nπ )(2nπ − 1)

[log(nπ )]2

1

3. Numerical method In this section, we derive the numerical method for the distributed-order space-fractional reaction-diffusion differential problem (1.1)-(1.3). We demonstrate the technique on the one-dimensional problem

∂u =− ∂t



2

κ (α ) −

∂2 ∂ x2

α /2 u dα + f (x, t , u ),

0 < x < L, 0 < t ≤ T ,

(3.1)

1

with nonhomogeneous Robin boundary conditions

u x (0, t ) − β1 u (0, t ) = h1 (t ),

0
u x ( L , t ) + β2 u ( L , t ) = h2 (t ),

0
(3.2)

and initial condition

u (x, 0) = u 0 (x),

0 < x < L,

2

(3.3)

where κ (α ) is a nonnegative weight function for 1 < α ≤ 2 satisfying κ (α ) ≡ 0 and 1 κ (α )dα < ∞, and β1 and β2 are non-negative. To approximate the integral in equation (3.1), let P > 0 be a fixed integer and 1 = ρ0 < ρ1 < · · · < ρ P = 2 be the equally spaced grid points discretizing the interval [1, 2]. Applying the midpoint quadrature rule with ρ = 1/ P and αi = (ρi −1 + ρi )/2, we obtain

2 −



∂2 κ (α ) − 2 ∂x

α /2 u dα = −

1

P 1

P

i =1



∂2 κ (αi ) − 2 ∂x

α2i

u + O ((ρ )2 ).

Thus, the distributed-order equation (3.1) is approximated by the following multi-term space-fractional equation

P ∂2 2 ∂u 1 =− κ (αi ) − 2 u + f (x, t , u ), ∂t P ∂x αi

i =1

0 < x < L, 0 < t ≤ T .

(3.4)

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

147

3.1. Spatial discretization via matrix transfer technique In this subsection, we apply the matrix transfer technique, introduced by Ilic´ et al. [21,22], to spatially discretize the multi-term fractional equation (3.4). First, consider the standard space reaction-diffusion initial boundary value problem

∂ 2u ∂u = − 2 + f (x, t , u ), ∂t ∂x u x (0, t ) − β1 u (0, t ) = h1 (t ), 0 < t ≤ T ,

0 < x < L, 0 < t ≤ T ,

(3.5) (3.6)

u x ( L , t ) + β2 u ( L , t ) = h2 (t ),

0
(3.7)

u (x, 0) = u 0 (x),

0 < x < L.

(3.8)

We use the method of lines (MOL) to derive a semi-discretized approximation of the problem (3.5)-(3.8). For a given positive integer M, let h = L / M be the spatial step size and define the spatial grid points xm = mh, m = 0, 1, . . . , M. Let um = u (xm , t ) and f m, = f (xm , t , um ) for m = 0, 1, . . . , M. The differential equation (3.5) is spatially discretized at each xm , m = 0, . . . , M using the second order central difference approximation

∂2 um+1 − 2um + um−1 u (xm , t ) = + O (h2 ) ∂ x2 h2 of the spatial derivatives. The derivative in the boundary conditions is approximated by

∂ u m +1 − u m −1 + O (h2 ) u (xm , t ) = ∂x 2h for m = 0, M. Using these approximations for the spatial derivatives, we obtain

du 0 dt

=

1 h2

2

([2 + 2hβ1 ]u 0 − 2u 1 ) + h1 (t ) + f 0 , h

(3.9)

for m = 0. For 1 ≤ m ≤ M − 1, we have

dum dt

=

1 h2

(−um−1 + 2um − um+1 ) + f m .

(3.10)

When m = M, we have

du M dt

=

1 h2

2

(−2u M −1 + [2 + 2hβ2 ]u M ) − h2 (t ) + f M . h

(3.11)

Also, from the initial condition (3.8), we have

um (0) = u 0 (xm ),

0 ≤ m ≤ M.

(3.12)

The equations (3.9)-(3.12) give a semi-discrete approximation of the standard space reaction-diffusion initial boundary value problem (3.5)-(3.8), which can be written in the matrix form as

du

2

2

h

h

= Au + h1 (t )e 0 − h2 (t )e M + f (t , u ),

dt u (0) = u 0 ,

(3.13)

where u = (u 0 , u 1 , . . . , u M ) T , f (t , u ) = ( f 0 , f 1 , . . . , f M ) T , u 0 = (u 0 (x0 ), u 0 (x1 ), . . . , u 0 (x M )) T , the ( M + 1)-dimensional veci ) for m = 0, . . . , M are defined by tors em = (em



i em =

1 for m = i , 0 for m = i ,

(3.14)

and the matrix A is the tridiagonal matrix of size ( M + 1) × ( M + 1) given by

⎛

2 + 2hβ1 ⎜ −1 ⎜

1 ⎜ A= 2⎜ h ⎜

⎝

−2 2

..

.

⎞ −1 .. . −1

..

.

2 −1 −2 2 + 2hβ2

⎟ ⎟ ⎟ ⎟. ⎟ ⎠

(3.15)

148

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

Next, we consider the equation (3.4) rewritten in the split form as

−1

P ∂2 ∂2 2 ∂u 1 − 2 u + f (x, t , u ). =− κ (αi ) − 2 ∂t P ∂x ∂x αi

(3.16)

i =1



2



The matrix representation of the Laplacian with homogeneous Robin boundary conditions is given by m − ∂∂x2 = A, where A is defined by equation (3.15). If the function u satisfies the nonhomogeneous boundary conditions (3.6)-(3.7), then

m −



∂ 2u ∂ x2

2

2

h

h

= Au + h1 (t )e 0 − h2 (t )e M . 

  2  αi −1 αi 2 ∂ − ∂ x2 = A 2 −1 for each αi , Assuming the fractional Laplacian under homogeneous boundary conditions satisfies m i = 1, 2, · · · , P , a semi-discrete approximation of the multi-term space-fractional differential equation (3.16) is given by

du dt

=− =−

P 1

P

i =1

P 1

P



αi

2

2

h

h



κ (αi ) A 2 −1 Au + h1 (t )e0 − h2 (t )e M + f (t , u ) P 1

αi

κ (αi ) A u − 2

P

i =1

κ (αi ) A

αi 2

−1



2 h

i =1

h1 (t )e 0 −

2 h

h2 (t )e M

+ f (t , u ).

Denoting

Aα =

P 1

P

αi

κ (αi ) A 2 ,

i =1

and

b(t ) = −

P 1

P

αi

κ (αi ) A 2 −1



2 h

i =1

h1 (t )e 0 −

2 h

h2 (t )e M

,

we can write the above system as

du

= − A α u + f (t , u ) + b(t ), dt u (0) = u 0 .

(3.17)

In a similar manner, we can derive the semi-discrete approximation of the one-dimensional distributed-order spacefractional reaction-diffusion equation (3.1) with nonhomogeneous Dirichlet boundary conditions

u (0, t ) = h1 (t ) and u ( L , t ) = h2 (t ). In this case, the matrix A is the ( M − 1) × ( M − 1) tridiagonal matrix given by A =

b(t ) =

P 1

P

αi

κ (αi ) A 2 −1



i =1

1

1

h

h

h (t )e 1 + 2 1

1 tridiag{−1, 2, −1} h2

and



h (t )e M −1 , 2 2

where the ( M − 1)-dimensional vectors e 1 and e M −1 are defined by (3.14). The semi-discrete approximation of the two-dimensional distributed-order space-fractional reaction-diffusion problem (1.1)-(1.3) can also be written in the form of the system (3.17). For the Dirichlet boundary conditions (1.3a), the matrix A is block tridiagonal given by A = h12 tridiag{− I , B , − I }, where B = tridiag{−1, 4, −1} and I is the identity matrix. When Robin boundary conditions (1.3c) are used, the matrix A is defined by

⎛

B0 ⎜ C1 ⎜

1 ⎜ . A = 2 ⎜ .. h ⎜ ⎝ 0 0 where

C0 B1

0 C1

..

..

. ··· 0

.

C M −1

···

··· ··· .. . B M −1 CM

0 0

⎞

⎟ ⎟ ⎟ , 0 ⎟ ⎟ ⎠ C M −1 BM

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

⎛

4 + 4hβ −1

⎜ ⎜ ⎜ Bi = ⎜ ⎜ ⎝

−2 4 + 2hβ ..

149

⎞ ⎟ −1 ⎟ ⎟ .. .. ⎟, . . ⎟ ⎠ −1 4 + 2hβ −1 −2 4 + 4hβ

.

for i = 0, M, and

⎛

4 + 2hβ ⎜ −1 ⎜

−2 4

⎜ Bi = ⎜ ⎜ ⎝

..

.

⎞ −1 .. . −1

..

.

4 −2

−1 4 + 2hβ

⎟ ⎟ ⎟ ⎟, ⎟ ⎠

for 1 ≤ i ≤ M − 1. The matrix C i = − I for i = 1, 2, . . . , M − 1, and C 0 = C M = −2I , where I is the identity matrix. Remark 3. The matrix A α can be constructed from the eigenvalues and eigenvectors of the matrix A. If A = H  H −1 , where  is the diagonal matrix containing the eigenvalues of the matrix A and H is the matrix of corresponding eigenvectors, then the fractional powers of A can be computed as A

Aα =

=

P 1

P

2

αi

= H  2 H −1 for i = 1, 2, · · · , P . Consequently,

αi

κ (αi ) A 2

i =1

P 1

P

αi

αi

κ (αi ) H  2 H −1

i =1

1

H  α H −1 , P αi P where α = i =1 κ (αi ) 2 .

=

3.2. Split-step predictor-corrector method Here, we present a split-step predictor-corrector time stepping method for the numerical solution of the semi-discretized problem (3.17). This method was discussed in [25] for nonlinear space-fractional reaction-diffusion equations. First, write the system (3.17) as

du

= F (t , u ), dt u (0) = u 0 ,

(3.18)

F (t , u ) = − A α u + f (t , u ) + b(t ).

(3.19)

where

Next, with the step size τ = T / N where N is a given positive integer, define the temporal grid points tn = nτ for n = 0, 1, 2, . . . , N and let un = u (tn ), bn = b(tn ), f n = f (tn , un ) and F n = F (tn , un ). Now, applying the split form of the second order Adams-Moulton formula (see [50,51]) to the ODE system (3.18), we obtain

u n +1 − u n = τ where φ = 0, predictor

u˜

n +1

1 2

and F˜



1 2

n+1

 − un = τ

Fn +

1 2



n +1 − φ F n+1 + φ τ F˜ ,

= F (tn+1 , u˜ n+1 ). Here u˜ n+1 is the predicted value of the solution obtained from a separate implicit 1 2





1 n +1 . + φ Fn + − φ F˜ 2

Using equation (3.19) in the above split-step predictor-corrector method, we obtain

150

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

u˜

n +1

 n

−u =τ

u n +1 − u n = τ

1 2



+φ



n

n

− Aα u + f + b

n



+

1 2

−φ



− A α u˜

n +1

n +1 + ˜f + b n +1

 ,



   1 − A α u n + f n + bn + − φ − A α u n +1 + f n +1 + b n +1

1 2

2



+ φ τ − A α u˜

n +1

+ ˜f

n +1

 + b n +1 ,

for n = 0, 1, . . . , N − 1. This fully implicit predictor-corrector method requires solving two nonlinear systems of equations at each time step. This can be avoided by treating the term f (t , u ) explicitly in both the predictor and the corrector. Moreover, n+1 we use the most current value available of f (t , u ). Thus, we replace f n+1 by f n in the predictor but replace f n+1 by ˜f in the corrector as this is the most current value available at this stage. This yields the following split-step scheme

T u˜





n +1



1



 n

n



1

n

1



n +1



, = I −τ + φ Aα u + τ f + τ +φ b + −φ b 2 2 2    1 τ n ˜ n +1  τ n T u n +1 = I − τ A α u n + f + f b + b n +1 , − φ τ A α u˜ n+1 + 2

2

2



(3.20)



for n = 0, 1, . . . , N − 1, where the coefficient matrix T = I + τ 12 − φ A α . The linearly implicit predictor-corrector scheme (3.20) reduces the computational cost significantly as it only requires solving two linear systems of equations with the same coefficient matrix at each time step. 3.3. Implementation

P

αi

2 Recall that the matrix A α = 1P i =1 κ (αi ) A . Here, A is the matrix representation of the standard Laplacian under ρ +ρ homogeneous boundary conditions and αi = i 2 i−1 for i = 1, 2, · · · , P , where 1 = ρ0 < ρ1 < · · · < ρ P = 2 is the equally spaced discretization of the interval [1, 2] with ρ = 1/ P . We construct the matrix A α from the eigenvalues and eigenvectors of the matrix A. If A = H  H −1 , where H and  are the matrices of eigenvectors and eigenvalues of the matrix

αi

αi

= H  2 H −1 for i = 1, 2, · · · , P . Consequently, the αi 2 i =1 κ (αi ) . The pseudocode implementing the split-step

A respectively, then the fractional powers of A can be computed as A

P

2

matrix A α is constructed as A α = 1P H α H −1 , where α = predictor-corrector method for solving nonlinear distributed-order space-fractional reaction-diffusion equation is given in Algorithm 1. Algorithm 1: Split-step predictor-corrector scheme. α = 0 for i = 1 to P do

αi =

ρi + ρi−1 2

 α =  α + κ (α i ) Aα =

1 P

αi 2

H α H −1

Precompute the inverse of the matrix T = I + τ for n = 0 to N − 1 do Compute RHS vector of predictor:





g˜ = I − τ Predict:

u˜

n+1

=

1 2

1 2

 − φ Aα





  1 1 + φ A α un + τ f n + τ + φ bn + − φ bn+1 2

2

T −1 g˜

Compute RHS vector of corrector:

 g= I− Correct:

1 2



τ A α un +

τ 2



n+1 f n + ˜f − φ τ A α u˜ n+1 +

τ 2

bn + bn+1



un+1 = T −1 g

3.4. IMEX predictor-corrector method In this subsection, we briefly describe another second-order method for solving the semi-discretized problem (3.17), which we use for comparison with the proposed split-step predictor-corrector method. IMEX predictor-corrector method

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

151

was introduced in [28] to solve integer order nonlinear parabolic equations. With the time stepsize temporal grid points tn = nτ for n = 0, 1, 2, . . . , N, the predictor is an IMEX method given by

τ = T / N and the

  τ n − A α u n − A α u n +1 + τ f n + b + b n +1 ,

τ

u n +1 − u n =

2

2

and the corrector is the Crank-Nicolson method given by

 τ   τ n − A α u n − A α u n +1 + f + f n +1 + b n + b n +1 .

τ

u n +1 − u n =

2

2

2

With M = I + τ2 A α , the above IMEX predictor-corrector method becomes

  τ  τ n = I − A α un + τ f n + b + b n +1 , 2 2   τ  τ  n ˜ n +1  τ  n Mun+1 = I − A α un + f + f b + b n +1 , + M u˜

n +1

2

2

(3.21)

2

for n = 0, 1, . . . , N − 1. 4. Stability and accuracy This section studies the stability and accuracy of the split-step predictor-corrector method (3.20). First, we analyze the stability of the method. A numerical method is stable if the spectral radius of the matrix representing the method does not exceed one. We use the Gerschgorin’s theorem (see [46]) to bound the spectrum of the matrix representing the proposed scheme (3.20). According to Gerschgorin’s theorem, each eigenvalue λ of a matrix A = ai , j n×n lies inside or on the boundary of at least one of the circles |λ − aii | =  j =i |ai , j |, i = 1, 2, . . . , n. Theorem 4. The split-step predictor-corrector method (3.20) is unconditionally stable for φ ≤ 14 . Proof. The predictor-corrector method (3.20) can be represented as

u n +1 = S u n + h n , where the matrix S is given by



 S = I +τ

1 2

− φ Aα

− 1 !

I−

1 2



τ Aα − φτ Aα I + τ



1 2

− φ Aα

−1

I −τ

1 2



"

+ φ Aα

n+1

and the vector hn includes the terms involving f n , ˜f , bn , and bn+1 . For the sake of simplicity, we present the proof for d = 1 with Robin boundary conditions (3.2). In this case, the matrix A α is defined by

Aα =

P 1

P

αi

κ (αi ) A 2 ,

i =1

P

αi

2 where the matrix A is given by the equation (3.15). Let λα be an eigenvalue of the matrix A α . Then λα = 1P i =1 κ (αi )λ , where λ is an eigenvalue of the matrix A. Since the off-diagonal elements of the matrix A are one-signed, all the eigenvalues of A are real. By Gerschgorin’s theorem (see [46]),

# # # # #λ − 2 + 2hβi # ≤ 2 # # 2 2 h

or

h

# # # # #λ − 2 # ≤ 2 , # h2 # h2

for i = 1, 2. Hence it follows that either

2β i h

≤λ≤

4 h2

+

2β i h

,

for i = 1, 2, or

0≤λ≤

4 h2

.

Since βi ≥ 0 (i = 1, 2), we must have

0≤λ≤

4 h2

+

2β h

,

152

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

P

αi

2 of the matrix where β = max{β1 , β2 }. Since κ (α ) is a nonnegative weight function, the eigenvalue λα = 1P i =1 κ (αi )λ A α is also nonnegative. The method (3.20) is unconditionally stable if the moduli of the eigenvalues of the matrix S do not exceed one, that is,

# # # 1 − 1 τ λ − φ τ λ 1 + τ  1 − φ  λ −1 1 − τ  1 + φ  λ # α α α α # # 2 2 2   # # ≤ 1. # # 1 + τ 1 − φ λα 2

After some algebraic manipulation, the stability condition for the split-step method (3.20) becomes

# # # 1 − 2φ τ λ + φ 2 + φ − 1  τ 2 λ2 # α # # α 4   # # ≤ 1. # 1 + (1 − 2φ)τ λα + φ 2 − φ + 1 τ 2 λ2α # 4

Since the eigenvalues λα of the matrix A α are nonnegative, the above stability condition is satisfied when φ ≤ 14 .

2

Next, we consider the accuracy of the scheme. We have the following result that shows that the predictor-corrector method has second order accuracy in time. Theorem 5. The split-step predictor-corrector method (3.20) is second order accurate in time. Proof. Substituting the exact solution u (t ) of (3.17) in the predictor formula of the scheme (3.20), we obtain





T u (tn+1 ) = I − τ

 +τ

1





+ φ A α u (tn ) + τ f (tn , u (tn ))



 1 1 + φ b(tn ) + − φ b(tn+1 ) + L p (u (tn ); τ ), 2

2

2

where L p (u (tn ); τ ) is the local truncation error of the predictor given by



L p (u (tn ); τ ) =

1 u (tn ) + 2



1 2



 1 − φ A α u (tn ) − − φ b (tn ) τ 2 + O (τ 3 ). 2

Now, subtracting the predictor formula from the above expression and using the localizing assumption, we obtain





I +τ

1 2



− φ Aα

Thus, it follows that

u (tn+1 ) − u˜

n +1

 =



u (tn+1 ) − u˜

1 u (tn ) + 2



1 2

n +1



= L p (u (tn ); τ ).



 1 − φ A α u (tn ) − − φ b (tn ) τ 2 + O (τ 3 ). 2

(4.1)

Similarly, for the corrector formula of the scheme (3.20), we have





I +τ

  τ n +1 − φ A α (u (tn+1 ) − un+1 ) = f (tn+1 , u (tn+1 )) − f (tn+1 , u˜ ) 2 2   − φ τ A α u (tn+1 ) − u˜ n+1 + L c (u (tn ); τ ), 1

where L c (u (tn ); τ ) is the local truncation error of the corrector given by

L c (u (tn ); τ ) = −

1 12

u (tn )τ 3 + O (τ 4 ).

Now

f (tn+1 , u (tn+1 )) − f (tn+1 , u˜

n +1

  ) = J (tn+1 , η¯ n+1 ) u (tn+1 ) − u˜ n+1 ,

n+1 where J (tn+1 , η¯ n+1 ) is the Jacobian matrix of f with respect to u evaluated at η¯ n+1 = θ u˜ + (1 − θ)u (tn+1 ), 0 < θ < 1. Thus, we obtain





I +τ

1 2

 

    1 n +1 − φ A α u (tn+1 ) − un+1 = τ + L c (u (tn ); τ ). J (tn+1 , η¯ n+1 ) − φ A α u (tn+1 ) − u˜ 2

On substituting equation (4.1) into the above equation, we obtain

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

u (tn+1 ) − un+1 =





1

J (tn+1 , η¯ n+1 ) − φ A α

2

1 1 u (tn ) + ( − φ) A α u (tn ) 2 2

 1 −( − φ)b (tn ) − u (tn ) τ 3 + O (τ 4 ). 1



2

153

12

2

Finally, we discuss the L-stability of the proposed scheme (3.20). Finite oscillations may arise in the numerical solution of the distributed-order space-fractional reaction-diffusion equation if the eigenvalues of the amplification matrix of the numerical method approach −1 as τ λα becomes large. L-stable methods have the property that the eigenvalues of the amplification matrix approach zero as τ λα becomes large due to which the numerical solution is oscillation free. To determine the L-stability, the numerical method is applied to the scalar problem

du dτ

= −λu ,

Re(λ) > 0,

resulting in a numerical approximation of the form un+1 = R ( z)un , where R is the rational function of z = τ λ called the stability function of the method. Definition 6. A numerical method is said to be L-stable if | R ( z)| < 1 for all z with Re( z) > 0 and | R ( z)| → 0 as Re( z) → ∞, where R ( z) is the stability function of the method. Although, the IMEX predictor-corrector method (3.21) is unconditionally stable, it is not L-stable. However, the split-step predictor-corrector method (3.20) is L-stable. Theorem 7. The split-step predictor-corrector method (3.20) is L-stable for φ = − 12 ±

√

2 . 2

The proof of this theorem is omitted as it is similar to the proof given in [25]. 5. Numerical results In this section, we consider several test problems to illustrate the effectiveness of the method proposed in this paper. We present numerical examples with nonhomogeneous boundary conditions. Problems with both smooth and nonsmooth initial conditions as well as with mismatch between initial and boundary conditions are considered. The proposed split-step predictor-corrector method is also compared with the second order IMEX predictor corrector method for accuracy and robustness. We denote the split-step predictor-corrector method (3.20) by SS-PC and the IMEX predictor-corrector method (3.21) by IMEX-PC. The discrete L 2 norm of a vector x = (x1 , x2 , . . . , x M ) T with respect to step size h is defined as

 ||x||2 = h

M 

1/2 | xi |

2

.

i =1

The convergence of a numerical method is demonstrated by computing the numerical order of convergence. When exact solution is available, the L 2 error is calculated as ||u − u h,τ ||2 and the numerical convergence rate is defined by

Rate = log2

||u − uh,τ ||2 , ||u − u h , τ ||2 2 2

where u is the exact solution and u h,τ is the numerical solution computed using space step size h and time step size τ . When the exact solution is not available, the L 2 error is computed as ||uh,τ − u h , τ ||2 and the numerical convergence rate 2

is defined as

Rate = log2

||uh,τ − u h , τ ||2 2 2

||u h , τ − u h , τ ||2 2 2

2

.

4 4

Remark 8. The split-step predictor-corrector method (3.20) is L-stable for φ = − 12 ± ical results are more accurate for φ = − 12 +

√

2 2

≈ 0.207 than for φ = − 12 −

√

2 2

√

2 . 2

In the examples tested, the numer-

≈ −1.207. Consequently, we will only report

the numerical results for the split-step predictor-corrector method (3.20) with φ = − 12 +

√

2 . 2

154

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

Fig. 1. Surface plots with h = 1/160,

τ = 1/80 and P = 800.

5.1. Example 1 We consider the linear one dimensional distributed order space-fractional reaction-diffusion equation

∂u =− ∂t



2

κ (α ) −

∂2 ∂ x2

α /2 u dα +

√

2et sin(π x),

0 < x < 1 , t > 0,

1

subject to nonhomogeneous Dirichlet boundary conditions

u (0, t ) = 1

and

u (1, t ) = 1,

and initial condition

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

0 < x < 1,

κ (α ) = α . The exact solution of this problem is derived in Section 2 (see equation (2.11)) and is given by √ ∞  2 [(−1)n − 1] −kn t u (x, t ) = 1 + (et − e −k1 t ) sin(π x) − 4 e sin(nπ x), k1 + 1 (nπ )3 2

n =1

where kn = 1 κ (α )(nπ )α dα is given by the equation (2.12). Fig. 1 shows the surface plots of the exact solution and the numerical solution computed by the two predictor-corrector methods using h = 1/160, τ = 1/80 and P = 800. Table 1 presents L 2 errors and the numerical rate of convergence of the two methods at t = 1. A second order convergence is observed for both methods.

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

155

Table 1 L 2 error and rate of convergence at t = 1.

τ

h

1/20 1/40 1/80 1/160 1/320

1/10 1/20 1/40 1/80 1/160

P

100 200 400 800 1600

SS-PC

IMEX-PC

L 2 error

Rate

L 2 error

Rate

5.16e-04 1.30e-04 3.26e-05 8.16e-06 2.04e-06

– 1.9903 1.9947 1.9972 1.9986

4.50e-04 1.07e-04 2.64e-05 6.56e-06 1.64e-06

– 2.0724 2.0218 2.0066 2.0021

Fig. 2. Surface plots with h = 1/160,

τ = 1/80 and P = 800.

5.2. Example 2 In this example, we solve the nonlinear one dimensional distributed-order space-fractional reaction-diffusion equation

∂u =− ∂t

2



κ (α ) −

∂2 ∂ x2

α /2 u dα + u 2 ,

0 < x < 1 , t > 0,

1

where

κ (α ) = α . The initial condition is the tent function given by ⎧ if 0.3 < x ≤ 0.5, ⎨ 25x − 7.5, u (x, 0) = −25x + 17.5, if 0.5 < x ≤ 0.7, ⎩ 0, otherwise,

and time-dependent Robin boundary conditions are taken to be

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

1 t+1

,

u x (1, t ) + u (1, t ) =

1 t2

+1

, t > 0.

Fig. 2 shows the surface plots of the numerical solution computed by the two methods using h = 1/160, τ = 1/80 and P = 800. The solution computed by the IMEX predictor-corrector method oscillates near the points where the initial condition is not smooth. The L-stable split-step method, on the other hand, produces oscillation free solution. L 2 errors and numerical rate of convergence of the two predictor-corrector schemes at t = 1 are presented in Table 2. The splitstep method once again demonstrates second order convergence, but the convergence rate of the IMEX predictor-corrector method appears to be less than two due to the oscillations in the computed numerical solution. 5.3. Example 3 We consider the nonlinear one dimensional distributed-order space-fractional equation

∂u =− ∂t

2



κ (α ) − 1

∂2 ∂ x2

α /2 u dα + u (1 − u ),

0 < x < 1 , t > 0,

156

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

Table 2 L 2 error and rate of convergence at t = 1.

τ

h

1/40 1/80 1/160 1/320 1/640

1/20 1/40 1/80 1/160 1/320

P

200 400 800 1600 3200

SS-PC

IMEX-PC

L 2 error

Rate

L 2 error

Rate

7.24e-02 1.70e-02 3.81e-03 7.83e-04 1.60e-04

– 2.0868 2.1591 2.2835 2.2875

8.48e-02 2.45e-02 7.02e-03 2.10e-03 6.40e-04

– 1.7920 1.8033 1.7418 1.7125

Fig. 3. Surface plots with h = 1/160,

τ = 1/80 and P = 800.

Table 3 L 2 error and rate of convergence at t = 1. h

1/40 1/80 1/160 1/320 1/640

where

τ 1/20 1/40 1/80 1/160 1/320

P

200 400 800 1600 3200

SS-PC

IMEX-PC

L 2 error

Rate

L 2 error

Rate

4.28e-03 1.28e-03 3.73e-04 9.78e-05 2.00e-05

– 1.7346 1.7841 1.9310 2.2935

3.11e-02 1.86e-02 1.12e-02 6.67e-03 3.98e-03

– 0.7396 0.7403 0.7424 0.7450

κ (α ) = − cos(πα /2). The initial condition is defined by the discontinuous step function ⎧ ⎨ 0, 0 ≤ x < 0.25, u (x, 0) = 1, 0.25 ≤ x < 0.75, ⎩ 0, 0.75 ≤ x ≤ 1,

and the time-dependent boundary conditions are given by

u x (0, t ) − u (0, t ) = et ,

u x (1, t ) + u (1, t ) = et , t > 0.

Fig. 3 gives the surface plots of the numerical solutions computed using h = 1/160, τ = 1/80 and P = 800. The solution produced by IMEX predictor-corrector method oscillates near the points of discontinuity of the initial condition, whereas the solution computed by the split-step predictor-corrector method is smooth and oscillation free. L 2 -errors and the numerical rate of convergence of the two methods are presented in Table 3. 5.4. Example 4 The final example of this section is a two dimensional distributed-order space-fractional equation

∂u =− ∂t

2

α

κ (α )(−) 2 udα + f (u ) 1

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

Fig. 4. Solution profiles at t = 1.5 with h = 0.05,

157

τ = 0.075, P = 80 when f (u ) = 0.

subject to time-dependent Dirichlet boundary conditions

u (0, y , t ) = e −t ,

u (1, y , t ) = e −t ,

u (x, 0, t ) = e −t ,

u (x, 1, t ) = e −t ,

and the initial condition

u (x, y , 0) = sin(π y ). In this example, κ (α ) = − cos(πα /2) and the spatial domain is  = [0, 1]2 . If f (u ) = 0, the exact solution of the problem can be derived in a manner similar to Example 1 and is given by

u (x, y , t ) = e −t + 2

∞  [1 − (−1)n ] n =1

nπ

e −kn1 t sin(nπ x) sin(π y )

∞  ∞  [1 − (−1)n ][1 − (−1)m ] [e −t − knm e −knm t ] sin(nπ x) sin(mπ y ), −4 1 − knm nmπ 2

2

n =1 m =1

α dα with λ2 = (nπ )2 + (mπ )2 . where knm = 1 κ (α )λnm nm

The solution profiles at t = 1.5 for f (u ) = 0 and f (u ) = u (1 − u ) are presented in Figs. 4 and 5 respectively. Due to the mismatch between the initial and the boundary conditions, the solutions computed using the IMEX method suffer from oscillations near the boundary of the spatial domain. However, the L-stable split-step predictor-corrector method produces smooth solutions in both cases. Tables 4 and 5 give L 2 -errors and the numerical rate of convergence of the two methods for f (u ) = 0 and f (u ) = u (1 − u ) respectively.

158

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

Fig. 5. Solution profiles at t = 1.5 with h = 0.05,

τ = 0.075, P = 80 when f (u ) = u (1 − u ).

Table 4 L 2 error and rate of convergence at t = 1.5 when f (u ) = 0. h

0.2 0.1 0.05 0.025

τ 0.3 0.15 0.075 0.0375

P

SS-PC

20 40 80 160

IMEX-PC

L 2 error

Rate

L 2 error

Rate

3.03e-04 6.82e-05 1.61e-05 3.97e-06

– 2.1524 2.0836 2.0180

1.42e-02 5.97e-03 2.41e-03 9.74e-04

– 1.2525 1.3073 1.3085

Table 5 L 2 error and rate of convergence at t = 1.5 when f (u ) = u (1 − u ). h

0.1 0.05 0.025 0.0125

τ 0.15 0.075 0.0375 0.01875

P

40 80 160 320

SS-PC

IMEX-PC

L 2 error

Rate

L 2 error

Rate

1.03e-03 2.73e-04 6.95e-05 1.74e-05

– 1.9101 1.9751 1.9940

1.10e-02 7.88e-03 3.13e-03 1.24e-03

– 0.4863 1.3318 1.3315

6. Conclusion In this paper, we considered the numerical approximation of nonlinear distributed-order space-fractional reactiondiffusion equations with time-dependent boundary conditions. Matrix transfer technique is employed for spatial discretization and a linearly implicit split-step predictor-corrector method based on Adams-Moulton formula is used for time-stepping to avoid solving nonlinear systems at each time step. Stability and second order convergence of the method are theoretically established. The order of convergence of the method is also demonstrated computationally. The proposed method is also compared with a second order IMEX predictor-corrector method. Although, the IMEX predictor-corrector method is unconditionally stable, it produced unwanted oscillations in the numerical solution when the initial data is nonsmooth or when there is a mismatch between the initial and the boundary conditions. The L-stable split-step predictor-corrector method, on the other hand, produced oscillation free and accurate solutions for these problems. Acknowledgements The authors would like to thank the referees for their very helpful comments and suggestions which greatly improved the quality of this paper. References [1] M. Abbaszadeh, Error estimate of second-order finite difference scheme for solving the Riesz space distributed-order diffusion equation, Appl. Math. Lett. 88 (2019) 179–185. [2] L. Aceto, P. Novati, Rational approximation to the fractional Laplacian operator in reaction-diffusion problems, SIAM J. Sci. Comput. 39 (1) (2017) A214–A228.

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

159

[3] E.E. Adams, L.W. Gelhar, Field study of dispersion in a heterogeneous aquifer: 2. Spatial moments analysis, Water Resour. Res. 28 (1992) 3293–3307. [4] U.M. Ascher, S.J. Ruuth, R.J. Spiteri, Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations, Appl. Numer. Math. 25 (1997) 151–167. [5] T.M. Atanackovic, S. Pilipovic, D. Zorica, Time distributed-order diffusion-wave equation. I. Volterra-type equation, Proc. R. Soc. A 465 (2009) 1893–1917. [6] D. Baleanu, A.K. Golmankhaneh, R. Nigmatullin, Fractional Newtonian mechanics, Cent. Eur. J. Phys. 8 (2010) 120–125. [7] T.A. Biala, Second-order predictor-corrector schemes for nonlinear distributed-order space-fractional differential equations with non-smooth initial data, Int. J. Comput. Math. (2018), https://doi.org/10.1080/00207160.2018.1539480. [8] K. Burrage, N. Hale, D. Kay, An efficient implicit FEM scheme for fractional-in-space reaction-diffusion equations, SIAM J. Sci. Comput. 34 (4) (2012) A2145–A2172. [9] W. Cao, F.H. Zeng, Z. Zhang, G.E. Karniadakis, Implicit-explicit difference schemes for nonlinear fractional differential equations with nonsmooth solutions, SIAM J. Sci. Comput. 38 (2016) A3070–A3093. [10] M. Caputo, Distributed order differential equations modelling dielectric induction and diffusion, Fract. Calc. Appl. Anal. 4 (2001) 421–442. [11] M. Caputo, Diffusion with space memory modelled with distributed order space fractional differential equations, Ann. Geophys. 46 (2003) 223–234. [12] A. Cartea, D. del Castillo-Negrete, Fractional diffusion models of option prices in markets with jumps, Physica A 374 (2) (2007) 749–763. [13] C. Celik, M. Duman, Crank-Nicolson method for the fractional diffusion equation with the Riesz fractional derivative, J. Comput. Phys. 231 (4) (2012) 1743–1750. [14] A. Chechkin, R. Gorenflo, I. Sokolov, Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations, Phys. Rev. E 66 (2002) 046129. [15] C.M. Chen, F. Liu, V. Anh, I. Turner, Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation, SIAM J. Sci. Comput. 32 (2010) 1740–1760. [16] S. Chen, X. Jiang, F. Liu, I. Turner, High order unconditionally stable difference schemes for the Riesz space-fractional telegraph equation, J. Comput. Appl. Math. 278 (2015) 119–129. [17] H.F. Ding, Y.X. Zhang, New numerical methods for the Riesz space fractional partial differential equations, Comput. Math. Appl. 63 (7) (2012) 1135–1146. [18] W. Fan, F. Liu, A numerical method for solving the two-dimensional distributed order space-fractional diffusion equation on an irregular convex domain, Appl. Math. Lett. 77 (2018) 114–121. [19] G. Gao, Z. Sun, A compact finite difference scheme for the fractional sub-diffusion equations, J. Comput. Phys. 230 (2011) 586–595. [20] W. Hundsdorfer, S.J. Ruuth, IMEX extensions of linear multistep methods with general monotonicity and boundedness properties, J. Comput. Phys. 225 (2007) 2016–2042. ´ F. Liu, I. Turner, V. Anh, Numerical approximation of a fractional-in-space diffusion equation I, Fract. Calc. Appl. Anal. 8 (2005) 323–341. [21] M. Ilic, ´ F. Liu, I. Turner, V. Anh, Numerical approximation of a fractional-in-space diffusion equation (II) with nonhomogeneous boundary conditions, [22] M. Ilic, Fract. Calc. Appl. Anal. 9 (2006) 333–349. [23] Z. Jiao, Y. Chen, I. Podlubny, Distributed-Order Dynamic Systems: Stability, Simulation, Applications and Perspectives, Springer, 2012. [24] J.T. Katsikadelis, Numerical solution of distributed order fractional differential equations, J. Comput. Phys. 259 (2014) 11–22. [25] K. Kazmi, A. Khaliq, A split-step predictor-corrector method for space-fractional reaction-diffusion equations with nonhomogeneous boundary conditions, Commun. Appl. Math. Comput. (2019), https://doi.org/10.1007/s42967-019-00030-z. [26] A.Q.M. Khaliq, T.A. Biala, S.S. Alzahrani, K.M. Faruti, Linearly implicit predictor-corrector methods for space-fractional reaction-diffusion equations with non-smooth initial data, Comput. Appl. Math. 75 (2018) 2629–2657. [27] X.Y. Li, B.Y. Wu, A numerical method for solving distributed order diffusion equations, Appl. Math. Lett. 53 (2016) 92–99. [28] D. Li, C. Zhang, W. Wang, Y. Zhang, Implicit-explicit predictor-corrector schemes for nonlinear parabolic differential equations, Appl. Math. Model. 35 (2011) 2711–2722. [29] J. Li, F. Liu, L. Feng, I. Turner, A novel finite volume method for the Riesz space distributed-order advection-diffusion equation, Appl. Math. Model. 46 (2017) 536–553. [30] J. Li, F. Liu, L. Feng, I. Turner, A novel finite volume method for the Riesz space distributed-order diffusion equation, Comput. Math. Appl. 74 (2017) 772–783. [31] B. Li, H. Luo, X. Xie, A time-spectral algorithm for fractional wave problems, J. Sci. Comput. 77 (2018) 1164–1184. [32] C. Lorenzo, T. Hartley, Variable order and distributed order fractional operators, Nonlinear Dyn. 29 (2002) 57–98. [33] S. Mashayekhi, M. Razzaghi, Numerical solution of distributed order fractional differential equations by hybrid functions, J. Comput. Phys. 315 (2016) 169–181. [34] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339 (2000) 1–77. [35] R. Metzler, J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A, Math. Gen. 37 (2004) R161–R208. [36] R. Metzler, T.F. Nonnenmacher, Space- and time-fractional diffusion and wave equations, fractional Fokker-Planck equations, and physical motivation, Chem. Phys. 284 (2002) 67–90. [37] M.L. Morgado, M. Rebelo, Numerical approximation of distributed order reaction-diffusion equations, J. Comput. Appl. Math. 275 (2015) 216–227. [38] M. Naber, Distributed order fractional subdiffusion, Fractals 12 (2004) 23–32. [39] E.D. Nezza, G. Palatucci, E. Valdinoci, Hitchhikers guide to the fractional Sobolev spaces, Bull. Sci. Math. 136 (5) (2012) 521–573. [40] R. Nigmatulin, The realization of the generalized transfer equation in a medium with fractal geometry, Phys. Status Solidi B 133 (1986) 425–430. [41] K.M. Owolab, K.C. Patidar, Higher-order time-stepping methods for time-dependent reaction-diffusion equations arising in biology, Appl. Math. Comput. 240 (2014) 30–50. [42] K.M. Owolab, E. Pindza, Modeling and simulation of nonlinear dynamical system in the frame of nonlocal and non-singular derivatives, Chaos Solitons Fractals 127 (2019) 146–157. [43] I. Podlubny, Fractional Differential Equations in Mathematics, Science and Engineering, vol. 198, Academic Press, 1999. [44] C. Pozrikidis, The Fractional Laplacian, CRC Press, 2016. [45] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives: Theory and Applications, Gordon and Breach Science Publishers, New York, 1987. [46] G.D. Smith, Numerical Solution of Partial Differential Equations: Finite Difference Methods, Clarendon Press, Oxford, 1985. [47] I. Sokolov, A. Chechkin, J. Klafter, Distributed-order fractional kinetics, Acta Phys. Pol. B 35 (2004) 1323–1341. [48] T. Srokowski, Lévy flights in nonhomogeneous media: distributed-order fractional equation approach, Phys. Rev. E 78 (2008) 031135. [49] C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comput. Phys. 213 (1) (2006) 205–213. [50] D.A. Voss, M.J. Casper, Efficient split linear multistep methods for stiff ODEs, SIAM J. Sci. Stat. Comput. 10 (1989) 990–999. [51] D.A. Voss, A.Q.M. Khaliq, A linearly implicit predictor-corrector method for reaction diffusion equations, Comput. Math. Appl. 38 (1999) 207–216. [52] X. Wang, F. Liu, X. Chen, Novel second-order accurate implicit numerical methods for the Riesz space distributed-order advection-dispersion equations, Adv. Math. Phys. 2015 (2015) 590435.

160

K. Kazmi, A.Q.M. Khaliq / Applied Numerical Mathematics 147 (2020) 142–160

[53] W. Wyss, The fractional Black-Scholes equations, Fract. Calc. Appl. Anal. 3 (1) (2000) 51–61. [54] Y. Yang, Y. Yan, N. Ford, Some time stepping methods for fractional diffusion problems with nonsmooth data, Comput. Methods Appl. Math. 18 (2018) 129–146. [55] 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. [56] S.B. Yuste, Weighted average finite difference methods for fractional diffusion equations, J. Comput. Phys. 216 (2006) 264–274. [57] H. Zhang, F. Liu, X. Jiang, F. Zeng, I. Turner, A Crank-Nicolson ADI Galerkin-Legendre spectral method for the two-dimensional Riesz space distributedorder advection-diffusion equation, Comput. Math. Appl. 76 (2018) 2460–2476.