- Email: [email protected]
Numerical approximation of nonlinear fractional parabolic differential equations with Caputo–Fabrizio derivative in Riemann–Liouville sense
Chaos, Solitons and Fractals 99 (2017) 171–179
Contents lists available at ScienceDirect
Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Numerical approximation of nonlinear fractional parabolic differential equations with Caputo–Fabrizio derivative in Riemann–Liouville sense Kolade M. Owolabi∗, Abdon Atangana Institute for Groundwater Studies, Faculty of Natural and Agricultural Sciences, University of the Free State, Bloemfontein 9300, South Africa
a r t i c l e
i n f o
Article history: Received 2 November 2016 Revised 4 April 2017 Accepted 4 April 2017 Available online 11 April 2017 MSC: 34A34 35A05 35K57 65L05 65M06 93C10
a b s t r a c t This paper considers the Caputo–Fabrizio derivative in Riemann–Liouville sense for the spatial discretization fractional derivative. We formulate two notable exponential time differencing schemes based on the Adams–Bashforth and the Runge–Kutta methods to advance the fractional derivatives in time. Our approach is tested on a number of fractional parabolic differential equations that are of current and recurring interest, and which cover pitfalls and address points and queries that may naturally arise. The effectiveness and suitability of the proposed techniques are justified via numerical experiments in one and higher dimensions. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Caputo–Fabrizio derivative Exponential decay-law Exponential time differencing Finite difference method Fractional nonlinear PDEs Numerical simulations Riemann–Liouville definition
1. Introduction
In this paper, we consider the general fractional-in-space reaction-diffusion equation of the form
The theory of fractional calculus or fractional derivative has a long standing history with applications in various scientific disciplines like chemistry, physics, signal and image processing, aerodynamics, economics, polymer rheology, electrodynamics, economics, control theory, biophysics, and blood flow phenomena [11,21,36,39]. Even though the discussion on derivatives of non-integer order could be dated back as far as the development of the classical theory of integer order differential calculus, Many researchers have investigated and still developing fractional calculus concepts in various ways as well as the numerical solutions of fractional-order differential equations using different numerical techniques [12,14,35].
∂ u(x, t ) = K ()α u(x, t ) + f (u(x, t )), t > 0, ∂t
∗
Corresponding author. E-mail addresses: [email protected] (K.M. Owolabi), [email protected] (A. Atangana). http://dx.doi.org/10.1016/j.chaos.2017.04.008 0960-0779/© 2017 Elsevier Ltd. All rights reserved.
(1.1)
subject to the initial condition u(x, 0 ) = u0 (x ), x ∈ [a, b] and the homogeneous Dirichlet or Neumann boundary conditions. Here, κ > 0 is the diffusion tensor or conductivity, α is called the laplacian operator with order α in superdiffusive domain 1 < α ≤ 2. The function f(u(x,t)) is nonlinear term that accounts for the local kinetics. Several numerical methods have been used to discretize the Laplacian operator of Eq. (1.1) [8,9,13,35,37,41]. Knowing well that the fractional differential operator is nonlocal in nature, which has posed a lot numerical and computational challenges other than the issue of stiffness encountered in the standard reaction-diffusion equations. In this paper, we utilize the Caputo–Fabrizio fractional derivative in Riemann–Liouville sense to discretize the space fractional derivative. We formulate some exponential time difference schemes based on Adams–Bashforth and Runge–Kutta schemes to advance in time.
172
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
The paper is organized as follows. In Section 2, we briefly review the main definitions concerning fractional differential equations and fractional derivatives. Section 3 presents the discretization scheme for the spatial fractional derivative of Caputo–Fabrizio in Riemnn-Liouville sense with the exponential decay law in the superdiffusive context. We also formulate some exponential time differencing schemes here for the temporal discretization. Numerical experiments in one, two and three dimensions are performed in Section 4. We conclude the paper with Section 5. 2. Some basic definitions of fractional differentiation There are many mathematical definitions about fractional derivative [7,21,36]. In relation with the Laplace transform, we show in this section some definitions of the Riemann–Liouville (R–L) fractional differentiation and their Fourier transform properties. Definition 2.1. The left and right R–L fractional derivatives of order α , (n − 1 < α < n ) of the function u(x) on a closed interval [a, b] are given as [36]
for 1 ≤ α ≤ 2 [26,36]. Definition 2.2. Assume that f is a differentiable function, Let α be a real number bounded in the interval 0 ≤ α ≤ 2, then the derivative with order α with power law is defined as [8,36] RL α 0 Dx
{ f ( x )} =
1 d (1 − α ) dx
α u (x ) =
1 dn (n − α ) dxn
x a
u (ξ ) dξ , (x − ξ )α−n+1
x
0
f (x )(x − τ )−α dτ .
(2.10)
Definition 2.3. Assume that f is a differentiable function, Let α be a real number bounded in the interval 0 ≤ α ≤ 2, then the derivative with order α with exponential decay law is defined as [8,36] AD α 0 Dx
{ f ( x )} =
B (α ) d 1 − α dx
x
0
f (x ) exp −α
t −x dx. 1−α
The Laplace transform of (2.11) is given as
L{
AD α 0 Dt
B (α ) d ( f (t ))}( p) = L 1 − α dt
b
t
(2.11)
t −x f (x ) exp −α dx ( p), 1−α
B(α ) pL{ f (t )}( p) = . 1 − α p + 1−αα
1. left R–L fractional derivative: a Dx
(2.12)
(2.2) 3. Numerical method of approximations
2. right R–L fractional derivative: α
b Db u
(x ) =
(−1 )n dn (n − α ) dxn
If α = n, then a Dα x u (x ) =
b
x
u (ξ ) dξ . (ξ − x )α−n+1
dn u (x ) dxn
(2.3)
d and x Dα u(x ) = (−1 )n dx n u ( x ). b n
If α > 0, u ∈ C0∞ (), ∈ R then, the Fourier transforms of the left and right R–L fractional derivatives are found satisfying
F (−∞ Dαx u(x )) = (iω )α uˆ (ω ),
and F (x Dα∞ u(x )) = (−iω )α uˆ (ω ), (2.4)
where the Fourier transform of u is denoted by uˆ (ω ),
uˆ (ω ) =
R
e−ωx u(x )dx.
The shifted Grünwald operator [25]
A h¯ ,ρ α u(x ) =
∞ 1 α gk u ( x − ( k − ρ ) h ¯ ), h ¯α
(2.5)
k=0
approximates the R–L fractional derivative uniformly with firstorder accuracy, that is,
Aαh¯ ,ρ u(x ) = −∞ Dαx u(x ) + O (h ),
(2.6)
In this section, we first consider methods of approximation of the spatial fractional derivative with exponential decay law using the Caputo–Fabrizio derivative in Riemann–Liouville sense in the superdiffusive interval 1 < α ≤ 2. We then formulate some second order exponential time differencing schemes based on Adams– Bashforth and Runge–Kutta methods to advance in time. 3.1. Approximation of space fractional derivative with exponential decay law Owing to the fact that most physical problems cannot be modelled by the power law, Caputo and Fabrizio [1,3,14] suggested a more vibrant approach of differentiation using the exponential decay law as an alternative to the power law kernel. This new development has caught many scholars attention, see [2,5,6,10,15,17–19]. Here, we briefly present the discretization scheme for the spatial fractional derivative of Caputo–Fabrizio in Riemann–Liouville sense with the exponential decay law in the superdiffusive region 1 < α ≤ 2. Definition 3.1. Let f be a function (not necessarily differentiable), then the Caputo–Fabrizio fractional derivative in Riemann–Liouville sense of order bounded in 1 < α ≤ 2 is defined as
where ρ is a positive integer and g(kα ) = (−1 )k (α ). Obviously, the k coefficients g(kα ) in Eq. (2.5) are the same as the coefficients of the power series function (1 − z )α , where
( 1 − z )α =
∞
∞ (−1 ) ( )z = g(kα ) zk , k α k
k
k=0
k=0
α + 1 k
g(kα−1) , k = 1, 2, . . . , .
g(1α ) = −α < 0,
(2.8)
⎫ ⎪ ⎬
1 ≥ g(2α ) ≥ g(3α ) ≥ · · · ≥ 0, ⎪ ∞ ( α ) ∞ ( α ) ⎭ = 0, < 0, m ≥ 1, k=0 gk k=0 gk
α d2 √ (1 − α ) π dx2 1 < α ≤ 2.
x 0
f (y ) exp −
α2 (x −y )2 dy, (1 − α ) (3.13)
In what follows, we suggest numerical approximation of this derivation, to achieve this, we let
The coefficients in (2.5) are also found satisfying the properties
g(0α ) = 1,
{ f (t )} =
(2.7)
for all |z| ≤ 1, which can be recursively evaluated using
g(0α ) = 1, g(kα ) = 1 −
CF R α 0 Dx
(2.9)
α d2 U (x ) = √ (1 − α ) π dx2
x 0
α2 2 f (y ) exp − (x − y ) dy, (1 − α ) (3.14)
where CF R α 0 Dx
{ f (t )} =
d 2 U (x ) . dx2
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
By replacing the spatial second derivative with a central finite difference approximation, we have
d 2U (x ) U (xi+1 ) − 2U (xi ) + U (xi−1 ) CF R α = . 0 Dt { f (x )} = dx2 2 ( x ) 2
(3.15)
With the above, Eq. (3.15) becomes
CF R α 0 Dt
1 { f ( x )} = 2(x )2
α U (x j+1 ) = √ (1 − α ) π
x j+1 0
f (τ ) exp −
2
α
(1 − α )
+
k=0
(x j+1 − τ ) dτ ,
2
j x k+1 α f (xk+1 ) √ (1 − α ) π k=0 xk α 2 × exp − (x j+1 − τ )2 dτ + G1 , (1 − α ) xk+1 α = √ f (xk+1 ) (1 − α ) π xk α 2 2 × exp − (x j+1 − τ ) dτ + G1 , (1 − α )
=
j f (xk+1 ) α ,1 H j,k + G1 , 2
Hαj,k,2 Hαj,k,3
Hαj,k,1 = erfc −α
x j+1 − xk+1 1−α
G =
x j+1 − xk , 1−α
and
G1 =
j x k+1 α ( f (τ ) − f (xk+1 )) √ (1 − α ) π k=0 xk α 2 2 × exp − (x j+1 − τ ) dτ . (1 − α )
j−1 f (xk+1 ) α ,2 H j,k + G2 2
(3.17)
k=0
where
Hαj,k,2 = erfc −α
x j − xk+1 1−α
− erfc −α
x j − xk , 1−α
and
G2 =
j−1 x k+1 α ( f (τ ) − f (x j )) √ (1 − α ) π k=0 xk α 2 2 × exp − (x j − τ ) d τ . (1 − α )
j−1 f (xk+1 ) α ,3 H j,k + G3 , 2
(3.18)
k=0
with
Hαj,k,3 = erfc −α
x j−1 − xk+1 1−α
− erfc −α
G3 =
α √ (1 − α ) π
× exp −
j−1 k=1
α
tk+1
tk
(1 − α )
( f (τ ) − f (x j−1 ))
2
x j−1 − xk , 1−α
and
(x j−1 − τ )2 dτ .
j x k+1 α ( f (τ ) − f (xk+1 )) √ (1 − α ) π k=0 xk j−1 x k+1 α 2 2 × exp − (x j+1 − τ ) dτ −2 ( f (τ ) − f (x j )) (1 − α ) xk k=0 j−1 t k+1 α 2 2 × exp − (x j − τ ) d τ + ( f (τ ) − f (x j−1 )) (1 − α ) t k=1 k α 2 × exp − (x j−1 − τ )2 dτ , (3.20) (1 − α )
j x 2 k+1 α − ( (1−αα ) ) (x j+1−τ )2 |G1 | = ( f (τ ) − f (xk+1 ))e dτ √ (1 − α ) π k=0 xk j x k+1 ( f (τ ) − f (x α k+1 ))(τ − xk+1 ) = √ (τ − xk+1 ) (1 − α ) π k=0 xk 2 − ( (1−αα ) ) (x j+1 −τ )2 × e dτ j x 2 α 2 k+1 − ( x − τ ) ( (1−α ) ) j+1 = f (y )(τ − xk+1 )e dτ , k=0 xk τ < y < x j+1
j x k+1 α x max f (x ) √ (1 − α ) π 0≤x≤x j+1 xk k=0 2 α × exp − (x j−1 − τ )2 dτ (1 − α ) j x −x x −x x j+1 j+1 k+1 k ≤ max f (x ) erfc −α −erfc −α , 2 0≤x≤x j+1 1−α 1−α k=0
x ( j + 1 )x ⇒ |G1 | ≤ max f (x ) efrc −α . (3.21) 2 0≤x≤x j+1 1−α
≤
Also,
U (x j−1 ) =
the erfc{·} denotes the error function. Conveniently, we now evaluate the first component of the reminder G as
Also, by adopting a similar technique, we get
U (x j ) =
+ G,
(3.19)
(3.16)
− erfc −α
k=0
k=0
where
k=0
x j+1 − xk+1 x j+1 − xk − erfc −α , 1−α 1−α x −x x −x j j k+1 k = erfc −α − erfc −α , 1−α 1−α x −x x −x j−1 j−1 k+1 k = erfc −α − erfc −α , 1−α 1−α
Hαj,k,1 = erfc −α
=
j j−1 f (xk+1 ) α ,1 f (xk+1 ) α ,2 H j,k − 2 H j,k 2 2
f (xk+1 ) Hαj,k,3 2 j
So also,
173
Without loss of generality, similar expressions can be obtained for the second and third terms as:
|G2 | ≤ and
|G3 | ≤
f (x )efrc −α ( j )x , 2 0≤x≤x j+1 1−α
x
max
f (x )efrc −α ( j − 1 )x . 2 0≤x≤x j+1 1−α
x
max
174
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
3.2. Exponential time differencing schemes (ETD) Let B = D (Lα ) ⊂ X, for 0 < α ≤ 1 with L represents a closed linear operator. The space B is a Banach space, Xα with norm υ B = Lα υ . We assume the following property on F holds: Let Y be open subset of R × Xα and F: Y → X holds: For every (t, u) ∈ Y there exists B ⊂ Y and constants c ≤ 0, 0 < θ ≤ 1 such that F (u1 , t1 ) − F (u2 , t2 ) c(|u1 − u2 |α + |t1 − t2 |θ ) for all (ui , ti ) ∈ B. From theory of standard semigroup [4,40] says that for every initial data (u0 , t0 ) ∈ Y, the initial value problem possess a unique local solution of the form
u ∈ C ([t0 , t ) : X ) ∩ C 1 ((t0 , t1 ) : X )
(3.22)
where t1 = t1 (u0 , t0 ) > t0 . It is assumed that the solution satisfies the integral formulation of our initial value problem which we reformulated into the recurrence form
u(tn+1 ) = e−Lh u(tn ) +
tn+1
tn
e−L(tn+1 −s ) F (u(s ), s )ds.
(3.23)
By putting s = tn + τ h with tn = nh(0 ≤ h ≤ h0 ; 0 ≤ n ≤ N ) and τ ∈ [0, 1] the recurrence relation is reduced to the form
u(tn+1 ) = e−Lh u(tn ) + h
1 0
e−L(1−τ F (u(tn + τ h ), tn + τ h )dτ , (3.24)
which form the basis for deriving the exponential time differencing methods. Eq. (3.24) is exact, method and various order ETD schemes come be derived from it depending on how one approximates the integral [12,16,27,28]. We focus primarily in this work on the second order ETD schemes, which adapt a linear approximation of the nonlinear function
F (u(tn + s ), tn + s ) ≈ F (u(tn ), tn )
+s
F (u(tn+1 , tn+1 )) − F (u(tn ), tn ) , h
un+1 = un eLh + Fn ((1 + hL )eLh − 1 − 2hL )/hLα + Fn−1 (−eLh + 1 − hL )/hLα .
(3.28)
Fractional second-order ETD method of Runge–Kutta type, denoted for brevity as ETD2RK:
un+1 = an + (F (an , tn + h ) − Fn )(eLh − 1 − hL )/hLα
(3.29)
where
an = un eLh + Fn (eLh − 1 )/L. To justify the accuracy and efficiency of these schemes, we consider a family of second order standard integrating factor (IF) methods whose formulations are also based on the Adams– Bashforth and Runge–Kutta schemes: Integrating factor method based on second-order Adams– Bashforth scheme (IF2AB):
un+1 = un eLh +
3h Lh h Fn e − Fn−1 e2Lh . 2 2
(3.30)
Integrating factor method based on second-order Runge–Kutta method (IF2RK)
un+1 = un eLh +
h (Fn eLh + F ((un + hFn )eLh , tn + h )). 2
(3.31)
Details of derivation of these schemes can be found in [16]. For any of these schemes to be asymptotically stable, then the ratio un+1 /un ≤ 1. For instance, to show that scheme (3.28) is asymptotically stable for 1 < α ≤ 2, we consider the test equation of the form
du(t ) = Lu(t ) + λu(t ). dt
(3.32)
Next, one writes Eq. (3.28) in the form
s = τ h ∈ [0, h],
Un+1 {((Lh + 1 )eLh − 1 − 2Lh )Fn + (−eLh + Lh + 1 )Fn−1 } = eLh + un Lα hun (3.33)
to yield the semidiscrete formula
The fractional second-order Adams–Bashforth, denoted as the ETD2ADAMS:
un+1 = e−Lh un + L−1 I − e−Lh F (un )
L−2 + hL − I + eLh [F (un+1 ) − F (un )]. h
On substituting Fn = λun into the equation above yields
(3.25)
By employing a locally second order approximation of un+1 :
u∗ = e−Lh un + L−1 (I − e−Lh )F (un ),
(3.26)
Un+1 {((Lh + 1 )eLh −1 −2Lh )λun + (−eLh + Lh + 1 )Fn−1 } = eLh + , un Lα hun (3.34)
the semidiscrete formula then results to
taking α x = λh, y = Lh and r = un+1 /un in the last expression we have
un+1 = e−Lh un + L−1 I − e−Lh F (un )
y2 r 2 − (y2 ey + [(y + 1 )ey −2y −1]x )r + (ey −y −1 )x = 0
L−2 + h L − I + e L h [ F ( u ∗ ) − F ( u n )] , h u∗ = e−Lh un + L−1 (I − e−Lh )F (un ).
therefore,
(3.27)
We assume here that the spatial operator L has been discretized to form a matrix through Caputo–Fabrizio fraction derivative in Riemann–Liouville sense or any of the classical techniques such as spectral, finite elements and finite difference methods [9,13,20,24,26,31,32,35], to mention a few. For the temporal discretization, we focus on two important second order ETD schemes whose formulations are based on the Adams–Bashforth and Runge–Kutta methods. This type of schemes have been applied in computational electrodynamics [34,38]. The following schemes can be used to efficiently integrate the fractional-in-space reaction-diffusion equation of superdiffusive order α for 1 < α ≤ 2 in time:
r = ey +
ey − 1 x+ y
( ey − y − 1 ) y2
Eq. (3.28) is asymptotically stable if
r = ey +
ey − 1 x+ y
( ey − y − 1 ) y2
(3.35)
x2 .
(3.36)
x2 ≤ 1.
Similarly, Eq. (3.29) can be written in the form
un+1 = eLh un +
{((Lh −2 )eLh + Lh + 2 )Fn + 2(eLh −Lh −1 )F (an , tn + h ) + h/2 )} un Lα h
.
(3.37)
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
175
Table 1 Numerical results of problem (4.40) for various schemes obtained at some different instances of fractional index α with final time t = 1 β = 0.75, D = 10 and h = 0.125.
α
ETD2AB
CPUT
ETD2RK
CPUT
IF2AB
CPUT
IF2RK
CPUT
1.1 1.3 1.5 1.7 1.9
2.4341e−05 3.0823e−05 3.2590e−05 4.0334e−05 4.5059e−05
0.174 0.174 0.172 0.183 0.173
5.0120e−08 1.2757e−07 1.7443e−07 1.9539e−07 1.9716e−07
0.184 0.175 0.173 0.180 0.185
3.5670e−05 4.2082e−05 4.8526e−05 5.4999e−05 6.1492e−05
0.181 0.184 0.173 0.180 0.176
4.2278e−05 5.0045e−05 5.7776e−05 6.5477e−05 7.3156e−05
0.191 0.169 0.172 0.173 0.170
Fig. 1. Stability regions of the ETD2 and ETD2RK schemes in the x − y plane.
Again, we let for simplicity α = 2, r = un+1 /un , x = λh and y = Lh to have
r = ey +
+
2(ey − y − 1 )ey/2 + (y − 2 )ey + y + 2 x y2
2(ey − y − 1 )(ey/2 − 1 ) 2 x . y3
(3.38)
r=e +
+
2(ey − y − 1 )ey/2 + (y − 2 )ey + y + 2 x y2
L∞ = max |ui − Ui |,
2(ey − y − 1 )(ey/2 − 1 ) 2 x ≤ 1. y3
i
(3.39)
Fig. 1 shows the stability regions of the schemes (3.28) and (3.29) in the x − y plane. The stability regions of schemes (3.30) and (3.31) follow similar process, see [16,33] for details. 4. Numerical experiments and discussion Applicability and reliability of numerical methods discussed in above section are examined here by considering some illustrative examples which are still of current and recurring interest in the fields of applied mathematics and physics to cover pitfalls and naturally arise questions. 4.1. Fractional diffusion-like models Under this section, we first justify the performance of the novel schemes by considering the numerical experiments of the space fractional one-dimensional problem with the homogeneous Dirichlet boundary conditions [13,35] given by
∂ u(x, t ) ∂ α u(x, t ) =D + f (u(x, t )), 1 < α ≤ 2, t > 0. ∂t ∂ xα
3 u(x, t ) = t α sin (2π x ), D f (u(x, t )) = {3[1 + (2π )α ] sin(2π x ) − [1 + (6π x )α ] sin(6π x )} 4 3 + αt α −1 sin (2π x ) − Du
with initial condition u(x, 0 ) = 100 exp(−100x2 ). The accuracy of these schemes are measured by reporting in Table 1, their error norm L∞ defined as:
The ETD2RK scheme (3.29) is asymptotically stable if y
The exact solution with the corresponding source term in x ∈ (0, 1) are given by
(4.40)
(4.41)
where ui and Ui are the respective ith exact and numerical solutions. It was observed that all the schemes performed well at α = 1.1 with computational time (CPUT) of less than 1-s, but the overall credit is given to the ETD2RK scheme. Hence, we experiment in high dimensions with the ETD2RK scheme. Again, we numerically simulate (4.40) by setting the source term
f (u(x, t )) = β u2 ,
(4.42)
so that the equation becomes fractional nonlinear KiSS model, where β is the plankton growth rate and D is the turbulent diffusivity, and the environmental conditions outside the patch are assumed unsuitable for the growth of plankton, so that u(0, t ) = u(L, t ) = 0 for all t and initial condition u(x, t) ≡ 0 for x, 0 and x > L. Here 0 < x < L, for L defines the patch radius [29,30]. Eq. (4.40) with (4.42) in classical form was first considered by Kierstead and Slobodkin (1953) [42] as problem of critical aggregation in application to the dynamics of the spread of plankton patches in marine ecosystems. It has been observed in various field of studies that small and large patches behave differently, while large patches tend to grow, the small patches tend to disappear. This type of behaviour occurs as a result of the interplay interaction between marine turbulence and plankton multiplication. We must ensure
176
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
Fig. 2. Numerical solution of the fractional KiSS model at some instances of fractional power α shown the oscillatory behaviour of the species at t = 400, D = 10 and β = 0.5.
Fig. 3. Surface plot showing the 2D spatiotemporal distribution of species at some instances of α = 1.1(0.1 )1.9 for the fractional Enzyme kinetic Eq. (4.43), with N = 300, β = 0.05, k = h = 0.25 and D = 0.5 via ETD2RK at final time t = 150.
that the patchy size L > 0 is chosen large enough to give room for the waves to propagate. The spatiotemporal oscillations displayed in Fig. 2 is obtained at parameter values L = 150, β = 0.5 and simulation time t = 150.
4.2. Enzyme kinetic model of Michaelis–Menten type Our primary interest is not really in one dimensional simulations, these are relatively easily undertaken using either method of lines coupled with spatial adaptive schemes. It is in higher dimensions that the ideas presented in this paper really become of serious value. We consider to illustrate the numerical algorithms proposed in Section 3 using a couple of non- trivial examples which we formulated in the form of space fractional reaction-diffusion equations. Here, we consider the two-dimensional fractional-in-space analogue of three-dimensional enzyme kinetics of Michaelis–Menten type of classical reaction-diffusion equation considered in [12] is
given as:
α ∂u u ∂ u ∂αu =D + −β ∂t ∂ xα ∂ yα 1+u
(4.43)
subject to the homogeneous Dirichlet boundary condition and the random the initial data (rand(N)∗ ρ ), where ρ is given as a small noise of orders in the interval [0.1, 0.5], and N defines the number points on the interval of integration. The simulation results in Fig. 3 which we obtained at different instances of fractional power α on interval 1 ≤ α ≤ 2, show different chaotic (unusual behaviour) and spatiotemporal structures arising from spurious oscillations in the solution. We extend the applicability of the schemes to solve threedimensional fractional enzyme kinetic problem (4.43). For the simulations, we utilize the initial condition:
u(x, y, z, 0 ) = exp(−20[(x − h ¯) ¯ /3 )2 + ( y − h ¯ /3 )2 + ( z − h ¯ /3 )2 ]/ h − exp(−20[(x − h ¯ /2 )2 + ( y − h ¯ /2 )2 + ( z − h ¯ /2 )2 ]/ h ¯) + exp(−20[(x− h ¯ )2 + (y− h ¯ )2 + (z− h ¯ )2 ]/ h ¯)
(4.44)
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
177
Fig. 4. Three dimensional simulation of fractional enzyme kinetic equation on a fixed computational domain (x, y, z) ∈ [0, L]3 , at various α values with D = 0.5, β = 0.5, h ¯ = 20 and final time t = 20.
Fig. 5. Numerical results for two dimensional fractional Brusselator system obtained at various α . Other parameters are: D1 = 0.05, D2 = 2, A = 1, B = 2 and final time t = 10 0 0.
Fig. 6. Numerical results for two dimensional fractional Brusselator system obtained at various α , with initial condition (4.47). Other parameters are: D1 = 0.05, D2 = 2, A = 1, B = 2 and final time t = 10 0 0. Take note of variation in the amplitudes.
subject to the boundary conditions set at the extreme ends of the computational domain of fixed size (x, y, z ) ∈ [0, L]3 , L = 20. In Fig. 4, different kinds of complex spatiotemporal structures that combine both oscillation and quasi-chaotic features are observed for the enzyme kinetic model. By increasing the fractional derivative order α in the equation results to even more complex spatiotemporal patterns. It should be noted that, apart from the patterns produced in Fig. 4, other phenomenon are obtainable depending on the choices of initial function and other parameter values. 4.3. Fractional Brusselator model Brusselator system in classical sense is considered as one of the most important reaction-diffusion equations used to describe the mechanisms of chemical reaction-diffusion with nonlinear
oscillations [12,23]. The Brusselator system consists of a pair of two chemical species, whose concentrations are governed by the given conditions. The Brusselator system occurs generally in the formation of ozone by atomic oxygen via a triple collision. It also used in the modelling of certain chemical reaction diffusion processes, like laser physics, plasma, and other enzymatic reactions. In this section, we begin our numerical simulations with the fractional two-dimensional generalization of the existing onedimensional classical Brusselator system [22]:
α ∂u ∂ u ∂αu = D1 + + u2 v − (A + 1 )u + B, ∂t ∂ xα ∂ yα α ∂v ∂ v ∂αv = D2 + − u2 v + Au ∂t ∂ xα ∂ yα
(4.45)
178
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
Fig. 7. Numerical results for three dimensional fractional Brusselator system obtained at various α at t = 5, with initial condition (4.49). Other parameters are given in Fig. 6. Take note of variation in the amplitudes.
subject to the initial conditions
u(x, y, 0 ) = 1 − 0.5 exp(−0.05(x2 + y2 )),
v(x, y, 0 ) = 0.25 exp(−0.05(x2 + y2 )).
(4.46)
The boundary conditions are set at the extremes of domain size ±L, L = 20. The 2D temporal surface distributions for chemical species u and v are displayed in Fig. 5. Obviously, the chemical species display a similar behaviour but in opposite direction. These type of behaviours may have good interpretation in biology and physics. It should be mentioned that depending on the choices of the initial conditions, the two chemical species can be made to distribute in phase in a similar fashion. For instance, if we considered the initial functions
u(x, y, 0 ) = 1 − 0.5 exp(−0.05(x2 + y2 )),
v(x, y, 0 ) = 0.25 exp(−0.05(x2 + y2 )).
(4.47)
The results at varying α are displayed in Fig. 6. It is obvious that both species distribute in phase in a similar fashion, readers are to take note of the variation in the amplitudes. Here, we also extend numerical simulation of the fractional Brusselator model to three dimensions by considering the generalization of existing one-dimensional reaction-diffusion system
α ∂u ∂ u ∂αu ∂αu = D1 + + + u2 v − (A + 1 )u + B, ∂t ∂ xα ∂ yα ∂ z α α ∂v ∂ v ∂αv ∂αu = D2 + + − u2 v + Au ∂t ∂ xα ∂ yα ∂ z α
this assertion is justified in Figs. 4–7. It should also be mentioned that, other spatiotemporal and chaotic patterns other than the ones reported in this paper are obtainable depending on the choices of the initial and parameter values. 5. Conclusion In this paper, we consider numerical solution of a range of fractional-in-space reaction diffusion equations. The spatial discretization of the fractional derivatives is formulated via the Caputo–Fabrizio derivative in Riemann–Liouville sense. For the Temporal discretization, we formulate two family of the exponential time differencing schemes in the fractional form. The stability regions of these schemes are verified and displayed. Accuracy and suitability of these schemes are examined by reporting the maximum norm error of fractional diffusion equation, in comparison with other two existing schemes of the same order. A range of spatiotemporal patterns emerged from the computer simulations in 1D, 2D and 3D which show the evolution of species in space. The results presented in this paper further suggest that pattern formation in the standard reaction-diffusion processes are practically the same as in the superdiffusive scenarios. Our assertion is justified in the results presented. The mathematical approach used here can be extended to solve multispecies problems of integer and noninteger orders in high dimensions. References
(4.48)
subject to the initial conditions
u(x, y, 0 ) = 1 − exp(−10(x − ν /2 )2 + (y − ν /2 )2 + (z − ν /2 )2 ),
v(x, y, 0 ) = exp(−10(x − ν /2 )2 + 2(y − ν /2 )2 + (z − ν /2 )2 ). (4.49) and the no flux boundary conditions fixed at ±L3 , L = 5. Spatial evolution of species u and v are shown in the first and second rows of Fig. 7 respectively, where a range of patterns such as star-like, diamond-like and oval-like shapes are obtained. Numerical simulation results have shown that complex dynamical patterns can be obtained in the super-diffusive (1 < α ≤ 2) fractional reaction-diffusion system other than the standard reaction-diffusion case when α = 2 . The work done in this paper have given insightful evidence through computer simulations that pattern formation in the classical reaction-diffusion models is practicably the same as in fractional reaction-diffusion scenarios,
[1] Algahtani OJJ. Comparing the Atangana-Baleanu and Caputo-Fabrizio derivative with fractional order: Allen cahn model. Chaos Solit Fract 2016;89:552–9. [2] Alkahtani BST. Chua’s circuit model with Atangana-Baleanu derivative with fractional order. Chaos Solit Fract 2016;89:547–51. [3] Alkahtani BST, Atangana A. Controlling the wave movement on the surface of shallow water with the Caputo-Fabrizio derivative with fractional order. Chaos Solit Fract 2016;89:539–46. [4] Asante-Asamani EO, Khaliq AQM, Wade BA. A real distinct poles exponential time differencing scheme for reaction-diffusion systems. J Comput Appl Math 2016;299:24–34. [5] Atangana A, Alkahtani BST. New model of groundwater owing within a confine aquifer: application of Caputo-Fabrizio derivative. Arabian J Geosci 2016;9:1–6. [6] Atangana A, Alqahtani RT. Numerical approximation of the space-time Caputo– Fabrizio fractional derivative and application to groundwater pollution equation. Adv Difference Equ 2016;2016(1):1–13. [7] Atangana A, Baleanu D. New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model. Thermal Sci 2016;20:763–9. [8] Atangana A, Koca I. Chaos in a simple nonlinear system with Atangana-Baleanu derivatives with fractional order. Chaos Solit Fract 2016;89:447–54. [9] Bakkyaraj T, Sahadevan R. Invariant analysis of nonlinear fractional ordinary differential equations with Riemann-Liouville fractional derivative. Nonlinear Dyn 2015;80:447–55.
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179 [10] Baleanu D, Caponetto R, Machado JT. Challenges in fractional dynamics and control theory. J Vib Control 2016;22:2151–2. [11] Baleanu D, Diethelm K, Scalas E. Fractional calculus: models and numerical methods, series on complexity, nonlinearity and chaos. World Scientific; 2012. [12] Bhatt HP, Khaliq AQM. The locally extrapolated exponential time differencing LOD scheme for multidimensional reaction-diffusion systems. J Comput Appl Math 2015;285:256–78. [13] Bueno-Orovio A, Kay D, Burrage K. Fourier spectral methods for fractional-in-space reaction-diffusion equations. BIT Numer Math 2014;54:937–54. [14] Caputo M, Fabrizio M. Applications of new time and spatial fractional derivatives with exponential kernels. Prog Fract Differ Appl 2016;2:1–11. [15] Coronel-Escamilla A, Gómez-Aguilar JF, López-López MG, Alvarado– Martínez VM, Guerrero-Ramírez GV. Triple pendulum model involving fractional derivatives with different kernels. Chaos Solit Fract 2016;91:248–61. [16] Cox SM, Matthews PC. Exponential time differencing for stiff systems. J Comput Phys 2002;176:430–55. [17] Gómez-Aguilar JF, Córdova-Fraga T, Escalante-Martínez JE, Calderón-Ramón C, Escobar-Jiménez RF. Electrical circuits described by a fractional derivative with regular kernel. Rev Mex Fis 2016;62:144–54. [18] López-López MG, Alvarado-Martínez VM, Reyes-Reyes J, Adam-Medina M, Gómez-Aguilar JF. Modeling diffusive transport with a fractional derivative without singular kernel. Physica A 2016;447:467–81. [19] Gómez-Aguilar JF, Torres L, Yépez-Martínez H, Baleanu D, Reyes JM, Sosa IO. Fractional liénard type model of a pipeline within the fractional derivative without singular kernel. Adv Difference Equ 2016;2016(1):1–13. [20] Hu Z, Wise SM, Wang C, Lowengrub JS. Stable and efficient finite-difference nonlinear-multigrid schemes for the phase field crystal equation. J Comput Phys 2009;228:5323–39. [21] Kilbas AA, Srivastava HM, Trujillo JJ. Theory and applications of fractional differential equations. Netherlands: Elsevier; 2006. [22] Kleefeld B, Khaliq AQM, Wade BA. An ETD cranknicolson method for reaction diffusion systems. Numer Methods Partial Differ Equ 2012;28:1309–35. [23] Lefever R, Nicolis G. Chemical instabilities and sustained oscillations. J Theor Biol 1971;30:267–84. [24] Liu F, Zhuang P, Anh V, Turner I, Burrage K. Stability and convergence of the finite difference method for the space-time fractional advection-diffusion equation. Appl Math Comput 2007;91:12–20. [25] Meerschaert MM, Tadjeran C. Finite difference approximations for fractional advection-dispersion flow equations. J Comput Appl Math 2004;172:65–77. [26] Meerschaert MM, Tadjeran C. Finite difference approximations for two-sided space-fractional partial differential equations. Appl Numer Math 2006;56:80–90.
179
[27] Owolabi KM, Patidar KC. Higher-order time-stepping methods for time-dependent reaction-diffusion equations arising in biology. Appl Math Comput 2014;240:30–50. [28] Owolabi KM, Patidar KC. Numerical solution of singular patterns in one-dimensional gray-scott-like models. Int J Nonlinear Sci Numer Simul 2014;15:437–62. [29] Owolabi KM, Patidar KC. Existence and permanence in a diffusive kiSS model with robust numerical simulations. Int J Differ Equ 2015;2015(485860):8. doi:10.1155/2015/485860. [30] Owolabi KM, Patidar KC. Effect of spatial configuration of an extended nonlinear Kierstead-Slobodkin reaction-transport model with adaptive numerical scheme. Springer Plus 2016;5:303. doi:10.1186/s40064- 016- 1941- y. [31] Owolabi KM, Atangana A. Numerical solution of fractional-in-space nonlinear Schrödinger equation with the riesz fractional derivative. Eur Phys J Plus 2016;131:335. doi:10.1140/epjp/i2016-16335-8. [32] Owolabi KM. Robust and adaptive techniques for numerical simulation of nonlinear partial differential equations of fractional order. Commun Nonlinear Sci Numer Simul 2017;44:304–17. [33] Owolabi KM. Mathematical study of two-variable systems with adaptive numerical methods. Numer Anal Appl 2016;9:218–30. [34] Petropoulos PG. Analysis of exponential time-differencing for FDTD in lossy dielectrics. IEEE Trans Antennas Propag 1997;45:1054. [35] Pindza E, Owolabi KM. Fourier spectral method for higher order space fractional reaction-diffusion equations. Commun Nonlinear Sci Numer Simul 2016;40:112–28. [36] Podlubny I. Fractional differential equations. San Diego: Academic Press; 1999. [37] Sousan E, Li C. A weighted finite difference method for the fractional diffusion equation based on the Riemann-Liouville derivative. Appl Numer Math 2015;90:22–37. [38] Taflove A. Computational electrodynamics: the finite-difference time-domain method. Artech House, London: Artech House Antenna Library; 1995. [39] Wu G, Baleanu D, Zeng S. Discrete chaos in fractional sine and standard maps. Phys Lett A 2014;378:484–7. [40] Zheng S. Nonlinear evolution equations. Chapman & Hall/CRC monographs and surveys in pure and applied mathematics; 2004. [41] Zheng M, Liu F, Turner I, Anh V. A novel high order space-time spectral method for the time fractional fokker-planck equation. SIAM J Sci Comput 2015;37:A701–24. [42] Kierstead H, Slobodkin LB. The size of water masses containing plankton looms. J Marine Res 1953;12:141–7.
Contents lists available at ScienceDirect
Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Numerical approximation of nonlinear fractional parabolic differential equations with Caputo–Fabrizio derivative in Riemann–Liouville sense Kolade M. Owolabi∗, Abdon Atangana Institute for Groundwater Studies, Faculty of Natural and Agricultural Sciences, University of the Free State, Bloemfontein 9300, South Africa
a r t i c l e
i n f o
Article history: Received 2 November 2016 Revised 4 April 2017 Accepted 4 April 2017 Available online 11 April 2017 MSC: 34A34 35A05 35K57 65L05 65M06 93C10
a b s t r a c t This paper considers the Caputo–Fabrizio derivative in Riemann–Liouville sense for the spatial discretization fractional derivative. We formulate two notable exponential time differencing schemes based on the Adams–Bashforth and the Runge–Kutta methods to advance the fractional derivatives in time. Our approach is tested on a number of fractional parabolic differential equations that are of current and recurring interest, and which cover pitfalls and address points and queries that may naturally arise. The effectiveness and suitability of the proposed techniques are justified via numerical experiments in one and higher dimensions. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Caputo–Fabrizio derivative Exponential decay-law Exponential time differencing Finite difference method Fractional nonlinear PDEs Numerical simulations Riemann–Liouville definition
1. Introduction
In this paper, we consider the general fractional-in-space reaction-diffusion equation of the form
The theory of fractional calculus or fractional derivative has a long standing history with applications in various scientific disciplines like chemistry, physics, signal and image processing, aerodynamics, economics, polymer rheology, electrodynamics, economics, control theory, biophysics, and blood flow phenomena [11,21,36,39]. Even though the discussion on derivatives of non-integer order could be dated back as far as the development of the classical theory of integer order differential calculus, Many researchers have investigated and still developing fractional calculus concepts in various ways as well as the numerical solutions of fractional-order differential equations using different numerical techniques [12,14,35].
∂ u(x, t ) = K ()α u(x, t ) + f (u(x, t )), t > 0, ∂t
∗
Corresponding author. E-mail addresses: [email protected] (K.M. Owolabi), [email protected] (A. Atangana). http://dx.doi.org/10.1016/j.chaos.2017.04.008 0960-0779/© 2017 Elsevier Ltd. All rights reserved.
(1.1)
subject to the initial condition u(x, 0 ) = u0 (x ), x ∈ [a, b] and the homogeneous Dirichlet or Neumann boundary conditions. Here, κ > 0 is the diffusion tensor or conductivity, α is called the laplacian operator with order α in superdiffusive domain 1 < α ≤ 2. The function f(u(x,t)) is nonlinear term that accounts for the local kinetics. Several numerical methods have been used to discretize the Laplacian operator of Eq. (1.1) [8,9,13,35,37,41]. Knowing well that the fractional differential operator is nonlocal in nature, which has posed a lot numerical and computational challenges other than the issue of stiffness encountered in the standard reaction-diffusion equations. In this paper, we utilize the Caputo–Fabrizio fractional derivative in Riemann–Liouville sense to discretize the space fractional derivative. We formulate some exponential time difference schemes based on Adams–Bashforth and Runge–Kutta schemes to advance in time.
172
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
The paper is organized as follows. In Section 2, we briefly review the main definitions concerning fractional differential equations and fractional derivatives. Section 3 presents the discretization scheme for the spatial fractional derivative of Caputo–Fabrizio in Riemnn-Liouville sense with the exponential decay law in the superdiffusive context. We also formulate some exponential time differencing schemes here for the temporal discretization. Numerical experiments in one, two and three dimensions are performed in Section 4. We conclude the paper with Section 5. 2. Some basic definitions of fractional differentiation There are many mathematical definitions about fractional derivative [7,21,36]. In relation with the Laplace transform, we show in this section some definitions of the Riemann–Liouville (R–L) fractional differentiation and their Fourier transform properties. Definition 2.1. The left and right R–L fractional derivatives of order α , (n − 1 < α < n ) of the function u(x) on a closed interval [a, b] are given as [36]
for 1 ≤ α ≤ 2 [26,36]. Definition 2.2. Assume that f is a differentiable function, Let α be a real number bounded in the interval 0 ≤ α ≤ 2, then the derivative with order α with power law is defined as [8,36] RL α 0 Dx
{ f ( x )} =
1 d (1 − α ) dx
α u (x ) =
1 dn (n − α ) dxn
x a
u (ξ ) dξ , (x − ξ )α−n+1
x
0
f (x )(x − τ )−α dτ .
(2.10)
Definition 2.3. Assume that f is a differentiable function, Let α be a real number bounded in the interval 0 ≤ α ≤ 2, then the derivative with order α with exponential decay law is defined as [8,36] AD α 0 Dx
{ f ( x )} =
B (α ) d 1 − α dx
x
0
f (x ) exp −α
t −x dx. 1−α
The Laplace transform of (2.11) is given as
L{
AD α 0 Dt
B (α ) d ( f (t ))}( p) = L 1 − α dt
b
t
(2.11)
t −x f (x ) exp −α dx ( p), 1−α
B(α ) pL{ f (t )}( p) = . 1 − α p + 1−αα
1. left R–L fractional derivative: a Dx
(2.12)
(2.2) 3. Numerical method of approximations
2. right R–L fractional derivative: α
b Db u
(x ) =
(−1 )n dn (n − α ) dxn
If α = n, then a Dα x u (x ) =
b
x
u (ξ ) dξ . (ξ − x )α−n+1
dn u (x ) dxn
(2.3)
d and x Dα u(x ) = (−1 )n dx n u ( x ). b n
If α > 0, u ∈ C0∞ (), ∈ R then, the Fourier transforms of the left and right R–L fractional derivatives are found satisfying
F (−∞ Dαx u(x )) = (iω )α uˆ (ω ),
and F (x Dα∞ u(x )) = (−iω )α uˆ (ω ), (2.4)
where the Fourier transform of u is denoted by uˆ (ω ),
uˆ (ω ) =
R
e−ωx u(x )dx.
The shifted Grünwald operator [25]
A h¯ ,ρ α u(x ) =
∞ 1 α gk u ( x − ( k − ρ ) h ¯ ), h ¯α
(2.5)
k=0
approximates the R–L fractional derivative uniformly with firstorder accuracy, that is,
Aαh¯ ,ρ u(x ) = −∞ Dαx u(x ) + O (h ),
(2.6)
In this section, we first consider methods of approximation of the spatial fractional derivative with exponential decay law using the Caputo–Fabrizio derivative in Riemann–Liouville sense in the superdiffusive interval 1 < α ≤ 2. We then formulate some second order exponential time differencing schemes based on Adams– Bashforth and Runge–Kutta methods to advance in time. 3.1. Approximation of space fractional derivative with exponential decay law Owing to the fact that most physical problems cannot be modelled by the power law, Caputo and Fabrizio [1,3,14] suggested a more vibrant approach of differentiation using the exponential decay law as an alternative to the power law kernel. This new development has caught many scholars attention, see [2,5,6,10,15,17–19]. Here, we briefly present the discretization scheme for the spatial fractional derivative of Caputo–Fabrizio in Riemann–Liouville sense with the exponential decay law in the superdiffusive region 1 < α ≤ 2. Definition 3.1. Let f be a function (not necessarily differentiable), then the Caputo–Fabrizio fractional derivative in Riemann–Liouville sense of order bounded in 1 < α ≤ 2 is defined as
where ρ is a positive integer and g(kα ) = (−1 )k (α ). Obviously, the k coefficients g(kα ) in Eq. (2.5) are the same as the coefficients of the power series function (1 − z )α , where
( 1 − z )α =
∞
∞ (−1 ) ( )z = g(kα ) zk , k α k
k
k=0
k=0
α + 1 k
g(kα−1) , k = 1, 2, . . . , .
g(1α ) = −α < 0,
(2.8)
⎫ ⎪ ⎬
1 ≥ g(2α ) ≥ g(3α ) ≥ · · · ≥ 0, ⎪ ∞ ( α ) ∞ ( α ) ⎭ = 0, < 0, m ≥ 1, k=0 gk k=0 gk
α d2 √ (1 − α ) π dx2 1 < α ≤ 2.
x 0
f (y ) exp −
α2 (x −y )2 dy, (1 − α ) (3.13)
In what follows, we suggest numerical approximation of this derivation, to achieve this, we let
The coefficients in (2.5) are also found satisfying the properties
g(0α ) = 1,
{ f (t )} =
(2.7)
for all |z| ≤ 1, which can be recursively evaluated using
g(0α ) = 1, g(kα ) = 1 −
CF R α 0 Dx
(2.9)
α d2 U (x ) = √ (1 − α ) π dx2
x 0
α2 2 f (y ) exp − (x − y ) dy, (1 − α ) (3.14)
where CF R α 0 Dx
{ f (t )} =
d 2 U (x ) . dx2
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
By replacing the spatial second derivative with a central finite difference approximation, we have
d 2U (x ) U (xi+1 ) − 2U (xi ) + U (xi−1 ) CF R α = . 0 Dt { f (x )} = dx2 2 ( x ) 2
(3.15)
With the above, Eq. (3.15) becomes
CF R α 0 Dt
1 { f ( x )} = 2(x )2
α U (x j+1 ) = √ (1 − α ) π
x j+1 0
f (τ ) exp −
2
α
(1 − α )
+
k=0
(x j+1 − τ ) dτ ,
2
j x k+1 α f (xk+1 ) √ (1 − α ) π k=0 xk α 2 × exp − (x j+1 − τ )2 dτ + G1 , (1 − α ) xk+1 α = √ f (xk+1 ) (1 − α ) π xk α 2 2 × exp − (x j+1 − τ ) dτ + G1 , (1 − α )
=
j f (xk+1 ) α ,1 H j,k + G1 , 2
Hαj,k,2 Hαj,k,3
Hαj,k,1 = erfc −α
x j+1 − xk+1 1−α
G =
x j+1 − xk , 1−α
and
G1 =
j x k+1 α ( f (τ ) − f (xk+1 )) √ (1 − α ) π k=0 xk α 2 2 × exp − (x j+1 − τ ) dτ . (1 − α )
j−1 f (xk+1 ) α ,2 H j,k + G2 2
(3.17)
k=0
where
Hαj,k,2 = erfc −α
x j − xk+1 1−α
− erfc −α
x j − xk , 1−α
and
G2 =
j−1 x k+1 α ( f (τ ) − f (x j )) √ (1 − α ) π k=0 xk α 2 2 × exp − (x j − τ ) d τ . (1 − α )
j−1 f (xk+1 ) α ,3 H j,k + G3 , 2
(3.18)
k=0
with
Hαj,k,3 = erfc −α
x j−1 − xk+1 1−α
− erfc −α
G3 =
α √ (1 − α ) π
× exp −
j−1 k=1
α
tk+1
tk
(1 − α )
( f (τ ) − f (x j−1 ))
2
x j−1 − xk , 1−α
and
(x j−1 − τ )2 dτ .
j x k+1 α ( f (τ ) − f (xk+1 )) √ (1 − α ) π k=0 xk j−1 x k+1 α 2 2 × exp − (x j+1 − τ ) dτ −2 ( f (τ ) − f (x j )) (1 − α ) xk k=0 j−1 t k+1 α 2 2 × exp − (x j − τ ) d τ + ( f (τ ) − f (x j−1 )) (1 − α ) t k=1 k α 2 × exp − (x j−1 − τ )2 dτ , (3.20) (1 − α )
j x 2 k+1 α − ( (1−αα ) ) (x j+1−τ )2 |G1 | = ( f (τ ) − f (xk+1 ))e dτ √ (1 − α ) π k=0 xk j x k+1 ( f (τ ) − f (x α k+1 ))(τ − xk+1 ) = √ (τ − xk+1 ) (1 − α ) π k=0 xk 2 − ( (1−αα ) ) (x j+1 −τ )2 × e dτ j x 2 α 2 k+1 − ( x − τ ) ( (1−α ) ) j+1 = f (y )(τ − xk+1 )e dτ , k=0 xk τ < y < x j+1
j x k+1 α x max f (x ) √ (1 − α ) π 0≤x≤x j+1 xk k=0 2 α × exp − (x j−1 − τ )2 dτ (1 − α ) j x −x x −x x j+1 j+1 k+1 k ≤ max f (x ) erfc −α −erfc −α , 2 0≤x≤x j+1 1−α 1−α k=0
x ( j + 1 )x ⇒ |G1 | ≤ max f (x ) efrc −α . (3.21) 2 0≤x≤x j+1 1−α
≤
Also,
U (x j−1 ) =
the erfc{·} denotes the error function. Conveniently, we now evaluate the first component of the reminder G as
Also, by adopting a similar technique, we get
U (x j ) =
+ G,
(3.19)
(3.16)
− erfc −α
k=0
k=0
where
k=0
x j+1 − xk+1 x j+1 − xk − erfc −α , 1−α 1−α x −x x −x j j k+1 k = erfc −α − erfc −α , 1−α 1−α x −x x −x j−1 j−1 k+1 k = erfc −α − erfc −α , 1−α 1−α
Hαj,k,1 = erfc −α
=
j j−1 f (xk+1 ) α ,1 f (xk+1 ) α ,2 H j,k − 2 H j,k 2 2
f (xk+1 ) Hαj,k,3 2 j
So also,
173
Without loss of generality, similar expressions can be obtained for the second and third terms as:
|G2 | ≤ and
|G3 | ≤
f (x )efrc −α ( j )x , 2 0≤x≤x j+1 1−α
x
max
f (x )efrc −α ( j − 1 )x . 2 0≤x≤x j+1 1−α
x
max
174
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
3.2. Exponential time differencing schemes (ETD) Let B = D (Lα ) ⊂ X, for 0 < α ≤ 1 with L represents a closed linear operator. The space B is a Banach space, Xα with norm υ B = Lα υ . We assume the following property on F holds: Let Y be open subset of R × Xα and F: Y → X holds: For every (t, u) ∈ Y there exists B ⊂ Y and constants c ≤ 0, 0 < θ ≤ 1 such that F (u1 , t1 ) − F (u2 , t2 ) c(|u1 − u2 |α + |t1 − t2 |θ ) for all (ui , ti ) ∈ B. From theory of standard semigroup [4,40] says that for every initial data (u0 , t0 ) ∈ Y, the initial value problem possess a unique local solution of the form
u ∈ C ([t0 , t ) : X ) ∩ C 1 ((t0 , t1 ) : X )
(3.22)
where t1 = t1 (u0 , t0 ) > t0 . It is assumed that the solution satisfies the integral formulation of our initial value problem which we reformulated into the recurrence form
u(tn+1 ) = e−Lh u(tn ) +
tn+1
tn
e−L(tn+1 −s ) F (u(s ), s )ds.
(3.23)
By putting s = tn + τ h with tn = nh(0 ≤ h ≤ h0 ; 0 ≤ n ≤ N ) and τ ∈ [0, 1] the recurrence relation is reduced to the form
u(tn+1 ) = e−Lh u(tn ) + h
1 0
e−L(1−τ F (u(tn + τ h ), tn + τ h )dτ , (3.24)
which form the basis for deriving the exponential time differencing methods. Eq. (3.24) is exact, method and various order ETD schemes come be derived from it depending on how one approximates the integral [12,16,27,28]. We focus primarily in this work on the second order ETD schemes, which adapt a linear approximation of the nonlinear function
F (u(tn + s ), tn + s ) ≈ F (u(tn ), tn )
+s
F (u(tn+1 , tn+1 )) − F (u(tn ), tn ) , h
un+1 = un eLh + Fn ((1 + hL )eLh − 1 − 2hL )/hLα + Fn−1 (−eLh + 1 − hL )/hLα .
(3.28)
Fractional second-order ETD method of Runge–Kutta type, denoted for brevity as ETD2RK:
un+1 = an + (F (an , tn + h ) − Fn )(eLh − 1 − hL )/hLα
(3.29)
where
an = un eLh + Fn (eLh − 1 )/L. To justify the accuracy and efficiency of these schemes, we consider a family of second order standard integrating factor (IF) methods whose formulations are also based on the Adams– Bashforth and Runge–Kutta schemes: Integrating factor method based on second-order Adams– Bashforth scheme (IF2AB):
un+1 = un eLh +
3h Lh h Fn e − Fn−1 e2Lh . 2 2
(3.30)
Integrating factor method based on second-order Runge–Kutta method (IF2RK)
un+1 = un eLh +
h (Fn eLh + F ((un + hFn )eLh , tn + h )). 2
(3.31)
Details of derivation of these schemes can be found in [16]. For any of these schemes to be asymptotically stable, then the ratio un+1 /un ≤ 1. For instance, to show that scheme (3.28) is asymptotically stable for 1 < α ≤ 2, we consider the test equation of the form
du(t ) = Lu(t ) + λu(t ). dt
(3.32)
Next, one writes Eq. (3.28) in the form
s = τ h ∈ [0, h],
Un+1 {((Lh + 1 )eLh − 1 − 2Lh )Fn + (−eLh + Lh + 1 )Fn−1 } = eLh + un Lα hun (3.33)
to yield the semidiscrete formula
The fractional second-order Adams–Bashforth, denoted as the ETD2ADAMS:
un+1 = e−Lh un + L−1 I − e−Lh F (un )
L−2 + hL − I + eLh [F (un+1 ) − F (un )]. h
On substituting Fn = λun into the equation above yields
(3.25)
By employing a locally second order approximation of un+1 :
u∗ = e−Lh un + L−1 (I − e−Lh )F (un ),
(3.26)
Un+1 {((Lh + 1 )eLh −1 −2Lh )λun + (−eLh + Lh + 1 )Fn−1 } = eLh + , un Lα hun (3.34)
the semidiscrete formula then results to
taking α x = λh, y = Lh and r = un+1 /un in the last expression we have
un+1 = e−Lh un + L−1 I − e−Lh F (un )
y2 r 2 − (y2 ey + [(y + 1 )ey −2y −1]x )r + (ey −y −1 )x = 0
L−2 + h L − I + e L h [ F ( u ∗ ) − F ( u n )] , h u∗ = e−Lh un + L−1 (I − e−Lh )F (un ).
therefore,
(3.27)
We assume here that the spatial operator L has been discretized to form a matrix through Caputo–Fabrizio fraction derivative in Riemann–Liouville sense or any of the classical techniques such as spectral, finite elements and finite difference methods [9,13,20,24,26,31,32,35], to mention a few. For the temporal discretization, we focus on two important second order ETD schemes whose formulations are based on the Adams–Bashforth and Runge–Kutta methods. This type of schemes have been applied in computational electrodynamics [34,38]. The following schemes can be used to efficiently integrate the fractional-in-space reaction-diffusion equation of superdiffusive order α for 1 < α ≤ 2 in time:
r = ey +
ey − 1 x+ y
( ey − y − 1 ) y2
Eq. (3.28) is asymptotically stable if
r = ey +
ey − 1 x+ y
( ey − y − 1 ) y2
(3.35)
x2 .
(3.36)
x2 ≤ 1.
Similarly, Eq. (3.29) can be written in the form
un+1 = eLh un +
{((Lh −2 )eLh + Lh + 2 )Fn + 2(eLh −Lh −1 )F (an , tn + h ) + h/2 )} un Lα h
.
(3.37)
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
175
Table 1 Numerical results of problem (4.40) for various schemes obtained at some different instances of fractional index α with final time t = 1 β = 0.75, D = 10 and h = 0.125.
α
ETD2AB
CPUT
ETD2RK
CPUT
IF2AB
CPUT
IF2RK
CPUT
1.1 1.3 1.5 1.7 1.9
2.4341e−05 3.0823e−05 3.2590e−05 4.0334e−05 4.5059e−05
0.174 0.174 0.172 0.183 0.173
5.0120e−08 1.2757e−07 1.7443e−07 1.9539e−07 1.9716e−07
0.184 0.175 0.173 0.180 0.185
3.5670e−05 4.2082e−05 4.8526e−05 5.4999e−05 6.1492e−05
0.181 0.184 0.173 0.180 0.176
4.2278e−05 5.0045e−05 5.7776e−05 6.5477e−05 7.3156e−05
0.191 0.169 0.172 0.173 0.170
Fig. 1. Stability regions of the ETD2 and ETD2RK schemes in the x − y plane.
Again, we let for simplicity α = 2, r = un+1 /un , x = λh and y = Lh to have
r = ey +
+
2(ey − y − 1 )ey/2 + (y − 2 )ey + y + 2 x y2
2(ey − y − 1 )(ey/2 − 1 ) 2 x . y3
(3.38)
r=e +
+
2(ey − y − 1 )ey/2 + (y − 2 )ey + y + 2 x y2
L∞ = max |ui − Ui |,
2(ey − y − 1 )(ey/2 − 1 ) 2 x ≤ 1. y3
i
(3.39)
Fig. 1 shows the stability regions of the schemes (3.28) and (3.29) in the x − y plane. The stability regions of schemes (3.30) and (3.31) follow similar process, see [16,33] for details. 4. Numerical experiments and discussion Applicability and reliability of numerical methods discussed in above section are examined here by considering some illustrative examples which are still of current and recurring interest in the fields of applied mathematics and physics to cover pitfalls and naturally arise questions. 4.1. Fractional diffusion-like models Under this section, we first justify the performance of the novel schemes by considering the numerical experiments of the space fractional one-dimensional problem with the homogeneous Dirichlet boundary conditions [13,35] given by
∂ u(x, t ) ∂ α u(x, t ) =D + f (u(x, t )), 1 < α ≤ 2, t > 0. ∂t ∂ xα
3 u(x, t ) = t α sin (2π x ), D f (u(x, t )) = {3[1 + (2π )α ] sin(2π x ) − [1 + (6π x )α ] sin(6π x )} 4 3 + αt α −1 sin (2π x ) − Du
with initial condition u(x, 0 ) = 100 exp(−100x2 ). The accuracy of these schemes are measured by reporting in Table 1, their error norm L∞ defined as:
The ETD2RK scheme (3.29) is asymptotically stable if y
The exact solution with the corresponding source term in x ∈ (0, 1) are given by
(4.40)
(4.41)
where ui and Ui are the respective ith exact and numerical solutions. It was observed that all the schemes performed well at α = 1.1 with computational time (CPUT) of less than 1-s, but the overall credit is given to the ETD2RK scheme. Hence, we experiment in high dimensions with the ETD2RK scheme. Again, we numerically simulate (4.40) by setting the source term
f (u(x, t )) = β u2 ,
(4.42)
so that the equation becomes fractional nonlinear KiSS model, where β is the plankton growth rate and D is the turbulent diffusivity, and the environmental conditions outside the patch are assumed unsuitable for the growth of plankton, so that u(0, t ) = u(L, t ) = 0 for all t and initial condition u(x, t) ≡ 0 for x, 0 and x > L. Here 0 < x < L, for L defines the patch radius [29,30]. Eq. (4.40) with (4.42) in classical form was first considered by Kierstead and Slobodkin (1953) [42] as problem of critical aggregation in application to the dynamics of the spread of plankton patches in marine ecosystems. It has been observed in various field of studies that small and large patches behave differently, while large patches tend to grow, the small patches tend to disappear. This type of behaviour occurs as a result of the interplay interaction between marine turbulence and plankton multiplication. We must ensure
176
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
Fig. 2. Numerical solution of the fractional KiSS model at some instances of fractional power α shown the oscillatory behaviour of the species at t = 400, D = 10 and β = 0.5.
Fig. 3. Surface plot showing the 2D spatiotemporal distribution of species at some instances of α = 1.1(0.1 )1.9 for the fractional Enzyme kinetic Eq. (4.43), with N = 300, β = 0.05, k = h = 0.25 and D = 0.5 via ETD2RK at final time t = 150.
that the patchy size L > 0 is chosen large enough to give room for the waves to propagate. The spatiotemporal oscillations displayed in Fig. 2 is obtained at parameter values L = 150, β = 0.5 and simulation time t = 150.
4.2. Enzyme kinetic model of Michaelis–Menten type Our primary interest is not really in one dimensional simulations, these are relatively easily undertaken using either method of lines coupled with spatial adaptive schemes. It is in higher dimensions that the ideas presented in this paper really become of serious value. We consider to illustrate the numerical algorithms proposed in Section 3 using a couple of non- trivial examples which we formulated in the form of space fractional reaction-diffusion equations. Here, we consider the two-dimensional fractional-in-space analogue of three-dimensional enzyme kinetics of Michaelis–Menten type of classical reaction-diffusion equation considered in [12] is
given as:
α ∂u u ∂ u ∂αu =D + −β ∂t ∂ xα ∂ yα 1+u
(4.43)
subject to the homogeneous Dirichlet boundary condition and the random the initial data (rand(N)∗ ρ ), where ρ is given as a small noise of orders in the interval [0.1, 0.5], and N defines the number points on the interval of integration. The simulation results in Fig. 3 which we obtained at different instances of fractional power α on interval 1 ≤ α ≤ 2, show different chaotic (unusual behaviour) and spatiotemporal structures arising from spurious oscillations in the solution. We extend the applicability of the schemes to solve threedimensional fractional enzyme kinetic problem (4.43). For the simulations, we utilize the initial condition:
u(x, y, z, 0 ) = exp(−20[(x − h ¯) ¯ /3 )2 + ( y − h ¯ /3 )2 + ( z − h ¯ /3 )2 ]/ h − exp(−20[(x − h ¯ /2 )2 + ( y − h ¯ /2 )2 + ( z − h ¯ /2 )2 ]/ h ¯) + exp(−20[(x− h ¯ )2 + (y− h ¯ )2 + (z− h ¯ )2 ]/ h ¯)
(4.44)
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
177
Fig. 4. Three dimensional simulation of fractional enzyme kinetic equation on a fixed computational domain (x, y, z) ∈ [0, L]3 , at various α values with D = 0.5, β = 0.5, h ¯ = 20 and final time t = 20.
Fig. 5. Numerical results for two dimensional fractional Brusselator system obtained at various α . Other parameters are: D1 = 0.05, D2 = 2, A = 1, B = 2 and final time t = 10 0 0.
Fig. 6. Numerical results for two dimensional fractional Brusselator system obtained at various α , with initial condition (4.47). Other parameters are: D1 = 0.05, D2 = 2, A = 1, B = 2 and final time t = 10 0 0. Take note of variation in the amplitudes.
subject to the boundary conditions set at the extreme ends of the computational domain of fixed size (x, y, z ) ∈ [0, L]3 , L = 20. In Fig. 4, different kinds of complex spatiotemporal structures that combine both oscillation and quasi-chaotic features are observed for the enzyme kinetic model. By increasing the fractional derivative order α in the equation results to even more complex spatiotemporal patterns. It should be noted that, apart from the patterns produced in Fig. 4, other phenomenon are obtainable depending on the choices of initial function and other parameter values. 4.3. Fractional Brusselator model Brusselator system in classical sense is considered as one of the most important reaction-diffusion equations used to describe the mechanisms of chemical reaction-diffusion with nonlinear
oscillations [12,23]. The Brusselator system consists of a pair of two chemical species, whose concentrations are governed by the given conditions. The Brusselator system occurs generally in the formation of ozone by atomic oxygen via a triple collision. It also used in the modelling of certain chemical reaction diffusion processes, like laser physics, plasma, and other enzymatic reactions. In this section, we begin our numerical simulations with the fractional two-dimensional generalization of the existing onedimensional classical Brusselator system [22]:
α ∂u ∂ u ∂αu = D1 + + u2 v − (A + 1 )u + B, ∂t ∂ xα ∂ yα α ∂v ∂ v ∂αv = D2 + − u2 v + Au ∂t ∂ xα ∂ yα
(4.45)
178
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179
Fig. 7. Numerical results for three dimensional fractional Brusselator system obtained at various α at t = 5, with initial condition (4.49). Other parameters are given in Fig. 6. Take note of variation in the amplitudes.
subject to the initial conditions
u(x, y, 0 ) = 1 − 0.5 exp(−0.05(x2 + y2 )),
v(x, y, 0 ) = 0.25 exp(−0.05(x2 + y2 )).
(4.46)
The boundary conditions are set at the extremes of domain size ±L, L = 20. The 2D temporal surface distributions for chemical species u and v are displayed in Fig. 5. Obviously, the chemical species display a similar behaviour but in opposite direction. These type of behaviours may have good interpretation in biology and physics. It should be mentioned that depending on the choices of the initial conditions, the two chemical species can be made to distribute in phase in a similar fashion. For instance, if we considered the initial functions
u(x, y, 0 ) = 1 − 0.5 exp(−0.05(x2 + y2 )),
v(x, y, 0 ) = 0.25 exp(−0.05(x2 + y2 )).
(4.47)
The results at varying α are displayed in Fig. 6. It is obvious that both species distribute in phase in a similar fashion, readers are to take note of the variation in the amplitudes. Here, we also extend numerical simulation of the fractional Brusselator model to three dimensions by considering the generalization of existing one-dimensional reaction-diffusion system
α ∂u ∂ u ∂αu ∂αu = D1 + + + u2 v − (A + 1 )u + B, ∂t ∂ xα ∂ yα ∂ z α α ∂v ∂ v ∂αv ∂αu = D2 + + − u2 v + Au ∂t ∂ xα ∂ yα ∂ z α
this assertion is justified in Figs. 4–7. It should also be mentioned that, other spatiotemporal and chaotic patterns other than the ones reported in this paper are obtainable depending on the choices of the initial and parameter values. 5. Conclusion In this paper, we consider numerical solution of a range of fractional-in-space reaction diffusion equations. The spatial discretization of the fractional derivatives is formulated via the Caputo–Fabrizio derivative in Riemann–Liouville sense. For the Temporal discretization, we formulate two family of the exponential time differencing schemes in the fractional form. The stability regions of these schemes are verified and displayed. Accuracy and suitability of these schemes are examined by reporting the maximum norm error of fractional diffusion equation, in comparison with other two existing schemes of the same order. A range of spatiotemporal patterns emerged from the computer simulations in 1D, 2D and 3D which show the evolution of species in space. The results presented in this paper further suggest that pattern formation in the standard reaction-diffusion processes are practically the same as in the superdiffusive scenarios. Our assertion is justified in the results presented. The mathematical approach used here can be extended to solve multispecies problems of integer and noninteger orders in high dimensions. References
(4.48)
subject to the initial conditions
u(x, y, 0 ) = 1 − exp(−10(x − ν /2 )2 + (y − ν /2 )2 + (z − ν /2 )2 ),
v(x, y, 0 ) = exp(−10(x − ν /2 )2 + 2(y − ν /2 )2 + (z − ν /2 )2 ). (4.49) and the no flux boundary conditions fixed at ±L3 , L = 5. Spatial evolution of species u and v are shown in the first and second rows of Fig. 7 respectively, where a range of patterns such as star-like, diamond-like and oval-like shapes are obtained. Numerical simulation results have shown that complex dynamical patterns can be obtained in the super-diffusive (1 < α ≤ 2) fractional reaction-diffusion system other than the standard reaction-diffusion case when α = 2 . The work done in this paper have given insightful evidence through computer simulations that pattern formation in the classical reaction-diffusion models is practicably the same as in fractional reaction-diffusion scenarios,
[1] Algahtani OJJ. Comparing the Atangana-Baleanu and Caputo-Fabrizio derivative with fractional order: Allen cahn model. Chaos Solit Fract 2016;89:552–9. [2] Alkahtani BST. Chua’s circuit model with Atangana-Baleanu derivative with fractional order. Chaos Solit Fract 2016;89:547–51. [3] Alkahtani BST, Atangana A. Controlling the wave movement on the surface of shallow water with the Caputo-Fabrizio derivative with fractional order. Chaos Solit Fract 2016;89:539–46. [4] Asante-Asamani EO, Khaliq AQM, Wade BA. A real distinct poles exponential time differencing scheme for reaction-diffusion systems. J Comput Appl Math 2016;299:24–34. [5] Atangana A, Alkahtani BST. New model of groundwater owing within a confine aquifer: application of Caputo-Fabrizio derivative. Arabian J Geosci 2016;9:1–6. [6] Atangana A, Alqahtani RT. Numerical approximation of the space-time Caputo– Fabrizio fractional derivative and application to groundwater pollution equation. Adv Difference Equ 2016;2016(1):1–13. [7] Atangana A, Baleanu D. New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model. Thermal Sci 2016;20:763–9. [8] Atangana A, Koca I. Chaos in a simple nonlinear system with Atangana-Baleanu derivatives with fractional order. Chaos Solit Fract 2016;89:447–54. [9] Bakkyaraj T, Sahadevan R. Invariant analysis of nonlinear fractional ordinary differential equations with Riemann-Liouville fractional derivative. Nonlinear Dyn 2015;80:447–55.
K.M. Owolabi, A. Atangana / Chaos, Solitons and Fractals 99 (2017) 171–179 [10] Baleanu D, Caponetto R, Machado JT. Challenges in fractional dynamics and control theory. J Vib Control 2016;22:2151–2. [11] Baleanu D, Diethelm K, Scalas E. Fractional calculus: models and numerical methods, series on complexity, nonlinearity and chaos. World Scientific; 2012. [12] Bhatt HP, Khaliq AQM. The locally extrapolated exponential time differencing LOD scheme for multidimensional reaction-diffusion systems. J Comput Appl Math 2015;285:256–78. [13] Bueno-Orovio A, Kay D, Burrage K. Fourier spectral methods for fractional-in-space reaction-diffusion equations. BIT Numer Math 2014;54:937–54. [14] Caputo M, Fabrizio M. Applications of new time and spatial fractional derivatives with exponential kernels. Prog Fract Differ Appl 2016;2:1–11. [15] Coronel-Escamilla A, Gómez-Aguilar JF, López-López MG, Alvarado– Martínez VM, Guerrero-Ramírez GV. Triple pendulum model involving fractional derivatives with different kernels. Chaos Solit Fract 2016;91:248–61. [16] Cox SM, Matthews PC. Exponential time differencing for stiff systems. J Comput Phys 2002;176:430–55. [17] Gómez-Aguilar JF, Córdova-Fraga T, Escalante-Martínez JE, Calderón-Ramón C, Escobar-Jiménez RF. Electrical circuits described by a fractional derivative with regular kernel. Rev Mex Fis 2016;62:144–54. [18] López-López MG, Alvarado-Martínez VM, Reyes-Reyes J, Adam-Medina M, Gómez-Aguilar JF. Modeling diffusive transport with a fractional derivative without singular kernel. Physica A 2016;447:467–81. [19] Gómez-Aguilar JF, Torres L, Yépez-Martínez H, Baleanu D, Reyes JM, Sosa IO. Fractional liénard type model of a pipeline within the fractional derivative without singular kernel. Adv Difference Equ 2016;2016(1):1–13. [20] Hu Z, Wise SM, Wang C, Lowengrub JS. Stable and efficient finite-difference nonlinear-multigrid schemes for the phase field crystal equation. J Comput Phys 2009;228:5323–39. [21] Kilbas AA, Srivastava HM, Trujillo JJ. Theory and applications of fractional differential equations. Netherlands: Elsevier; 2006. [22] Kleefeld B, Khaliq AQM, Wade BA. An ETD cranknicolson method for reaction diffusion systems. Numer Methods Partial Differ Equ 2012;28:1309–35. [23] Lefever R, Nicolis G. Chemical instabilities and sustained oscillations. J Theor Biol 1971;30:267–84. [24] Liu F, Zhuang P, Anh V, Turner I, Burrage K. Stability and convergence of the finite difference method for the space-time fractional advection-diffusion equation. Appl Math Comput 2007;91:12–20. [25] Meerschaert MM, Tadjeran C. Finite difference approximations for fractional advection-dispersion flow equations. J Comput Appl Math 2004;172:65–77. [26] Meerschaert MM, Tadjeran C. Finite difference approximations for two-sided space-fractional partial differential equations. Appl Numer Math 2006;56:80–90.
179
[27] Owolabi KM, Patidar KC. Higher-order time-stepping methods for time-dependent reaction-diffusion equations arising in biology. Appl Math Comput 2014;240:30–50. [28] Owolabi KM, Patidar KC. Numerical solution of singular patterns in one-dimensional gray-scott-like models. Int J Nonlinear Sci Numer Simul 2014;15:437–62. [29] Owolabi KM, Patidar KC. Existence and permanence in a diffusive kiSS model with robust numerical simulations. Int J Differ Equ 2015;2015(485860):8. doi:10.1155/2015/485860. [30] Owolabi KM, Patidar KC. Effect of spatial configuration of an extended nonlinear Kierstead-Slobodkin reaction-transport model with adaptive numerical scheme. Springer Plus 2016;5:303. doi:10.1186/s40064- 016- 1941- y. [31] Owolabi KM, Atangana A. Numerical solution of fractional-in-space nonlinear Schrödinger equation with the riesz fractional derivative. Eur Phys J Plus 2016;131:335. doi:10.1140/epjp/i2016-16335-8. [32] Owolabi KM. Robust and adaptive techniques for numerical simulation of nonlinear partial differential equations of fractional order. Commun Nonlinear Sci Numer Simul 2017;44:304–17. [33] Owolabi KM. Mathematical study of two-variable systems with adaptive numerical methods. Numer Anal Appl 2016;9:218–30. [34] Petropoulos PG. Analysis of exponential time-differencing for FDTD in lossy dielectrics. IEEE Trans Antennas Propag 1997;45:1054. [35] Pindza E, Owolabi KM. Fourier spectral method for higher order space fractional reaction-diffusion equations. Commun Nonlinear Sci Numer Simul 2016;40:112–28. [36] Podlubny I. Fractional differential equations. San Diego: Academic Press; 1999. [37] Sousan E, Li C. A weighted finite difference method for the fractional diffusion equation based on the Riemann-Liouville derivative. Appl Numer Math 2015;90:22–37. [38] Taflove A. Computational electrodynamics: the finite-difference time-domain method. Artech House, London: Artech House Antenna Library; 1995. [39] Wu G, Baleanu D, Zeng S. Discrete chaos in fractional sine and standard maps. Phys Lett A 2014;378:484–7. [40] Zheng S. Nonlinear evolution equations. Chapman & Hall/CRC monographs and surveys in pure and applied mathematics; 2004. [41] Zheng M, Liu F, Turner I, Anh V. A novel high order space-time spectral method for the time fractional fokker-planck equation. SIAM J Sci Comput 2015;37:A701–24. [42] Kierstead H, Slobodkin LB. The size of water masses containing plankton looms. J Marine Res 1953;12:141–7.