Numerical solutions of space-fractional advection-diffusion equations with nonlinear source term

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.1 (1-10)

Applied Numerical Mathematics ••• (••••) •••–•••

Contents lists available at ScienceDirect

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

Numerical solutions of space-fractional advection-diffusion equations with nonlinear source term Alessandra Jannelli a,∗ , Marianna Ruggieri b , Maria Paola Speciale a a

Department of Mathematical and Computer Sciences, Physical Sciences and Earth Sciences, University of Messina, Viale F. Stagno d’Alcontres 31, 98166 Messina, Italy b Faculty of Engineering and Architecture, Kore University of Enna, Via delle Olimpiadi, Cittadella Universitaria, 94100, Enna, Italy

a r t i c l e

i n f o

Article history: Received 30 January 2019 Received in revised form 9 January 2020 Accepted 18 January 2020 Available online xxxx Keywords: Fractional derivatives Advection-diffusion-reaction equations Lie symmetries Implicit trapezoidal method

a b s t r a c t In this paper, numerical solutions of space-fractional advection-diffusion equations, involving the Riemann-Liouville derivative with a nonlinear source term, are presented. We propose a procedure that combines the fractional Lie symmetries analysis, to reduce the original fractional partial differential equations into fractional ordinary differential equations, with a numerical method. By adopting the Caputo definition of derivative, the reduced fractional ordinary equations are solved by applying the implicit trapezoidal method. The numerical results confirm the applicability and the efficiency of the proposed approach. © 2020 IMACS. Published by Elsevier B.V. All rights reserved.

1. Introduction In recent years, the fractional partial differential equations (FPDEs) have played an important role in the modeling of physical and natural phenomena arising in many areas of the applied sciences, for the description of nonlinear phenomena such as diffusion processes and wave propagation problems. Several methodologies have been proposed to find analytical and numerical solutions. Efficient analytical and semianalytical methods, as the Adomian decomposition methods [4,3], the classical Laplace transform for linear fractional equations [38,29], the homotopy perturbation method [12,37] and the variational iteration method [30] have been developed. Many numerical methods have also been implemented, including the finite difference methods [35,7], the finite volume methods [13], the finite element methods [40,41] and the spectral methods [6]. As for the integer order differential equations, finite difference methods are one of the most important classes of numerical methods for solving FDEs. In [33], several numerical methods are presented and some numerical approximations to solutions to different fractional differential equations are given, experimentally verified on various examples. In [43], an implicit difference approximation to solve time-fractional diffusion equations is obtained. In [26,25], numerical solutions by finite difference/spectral approximations for fractional diffusion equations are proposed. Liu et al. [27] developed an explicit difference method and an implicit difference method for solving a space–time fractional advection dispersion equation on a finite domain. In [24,2,23], high order numerical difference schemes were constructed in order to solve the fractional advection–diffusion equations. Zhang et al. [42] obtained a finite difference method for FPDEs on a non-uniform mesh.

*

Corresponding author. E-mail addresses: [email protected] (A. Jannelli), [email protected] (M. Ruggieri), [email protected] (M.P. Speciale).

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

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.2 (1-10)

A. Jannelli et al. / Applied Numerical Mathematics ••• (••••) •••–•••

2

As well known, the Lie symmetries of a differential equation are a powerful tool for the determination of similarity solutions, conservation laws and so on [15,16,36]. So that, some extensions of the Lie symmetries to FPDEs have been widely explored and some approaches have been proposed in [9–11,22]. In this paper, we study the following space-fractional advection-diffusion (SF–ADR) equation with a nonlinear source term β+1

∂t u (t , x) − k1 ∂x β+1

where ∂x

β+1

∂x

u (t , x) + k2 ∂x u (t , x) + f (t , x, u ) = 0,

0 < β < 1,

(1)

is the Riemann-Liouville (R-L) fractional derivative operator

1

∂ u (t , x) = ∂x ∂x u (t , x) = (1 − β) ∂ x β

x 0

∂r u (t , r ) dr (x − r )β

(2)

u (t , x) being the field variable with t and x independent variables, k1 > 0 the diffusion coefficient, k2 = 0 the advection coefficient, and f = f (t , x, u ) the nonlinear source term. In order to get analytical and numerical solutions of the FPDE (1), we apply a procedure recently proposed in [17–19]. It combines the Lie symmetry analysis (developed by Gazizov et al. [9–11,22]) that reduces the original FPDE into a fractional ordinary differential equation (FODE), with the numerical methods that solve the obtained FODE. We solve the reduced β + 1 order FODE that allows to get solution of SF–ADR with nonlinear source term, defining, in this way, a wider class of problem, unlike the results obtained in [17–19], where the authors solved time fractional, space fractional and time and space fractional advection-diffusion equations by solving reduced β order FODEs. When the reduced equation does not admit exact solutions, a numerical method is considered. In this case, we introduce the Caputo’s derivative, that permits to rewrite the FODE in terms of this derivative so that we are able to determine numerical solutions, assigning initial conditions having a physical meaning. The numerical solutions are obtained by applying the implicit unconditionally stable trapezoidal method. We remark that the proposed approach allows to determine solutions of FPDEs that may be difficult to find by their direct integration. The strong point of the approach is the low computational effort: it is evident that to solve a fractional ordinary differential equation leads to a lower cost respect to the integration of a FPDE. The plan of the paper is as follows. In Section 2, we present the Lie symmetry theory and the Lie transformations admitted by the model that reduce the SF–ADR equation (1) into β + 1 order FODEs. In Section 3, we describe the numerical approach and, finally, in the last Section some numerical examples are reported in order to confirm the applicability and robustness of the proposed procedure. 2. Lie symmetries and reduction into FODEs In this Section, we recall the main results concerning the Lie theory applied to FPDEs and then we introduce the Lie transformations to reduce the original equation FPDE (1) into a FODE. 2.1. Lie symmetries for FPDEs This Section is devoted to some basic definitions concerning the Lie symmetry for the DEs and its extension to FDEs. Assigned a differential equation

(t , x, u , . . . , u (k) ) = 0,

(3)

where t and x are the independent variables, u is the dependent variable, and u (k) are all partial derivatives of the u up to order k, we recall that an invertible transformation of the variables t , x, u, depending on a continuous parameter a,

T = T (t , x, u , a),

X = X (t , x, u , a),

U = U (t , x, u , a),

(4)

is said to be an one-parameter (a) Lie point symmetry transformation of the equation (3) if the equation (3) has the same form in the new variables T , X and U . According to the Lie theory, by expanding the finite transformation (4) in Taylor’s series around a = 0, we get the infinitesimal transformation

T = t + a ξ1 (t , x, u ) + o(a2 ), where the infinitesimals ξ1 , ξ2 and

X = x + a ξ2 (t , x, u ) + o(a2 ),

U = u + a η(t , x, u ) + o(a2 )

η characterize the infinitesimal operator [16,32,31],

 = ξ1 (t , x, u )∂t + ξ2 (t , x, u )∂x + η(t , x, u )∂u . The infinitesimals ξ1 , ξ2 and

k  = 0|=0 .

(5)

(6)

η are found by means of the Lie’s algorithm, requiring that (7)

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.3 (1-10)

A. Jannelli et al. / Applied Numerical Mathematics ••• (••••) •••–•••

3

The k is the k-order prolongation of the operator (6) acting on (3), that involves the prolonged infinitesimals ηut , ηu x , upon to ηu (k) , linked to the transformations of all the partial derivatives ut , u x , . . . , u (k) , given by prolongation formulae

η u t = D t η − u t D t ξ1 − u x D t ξ2 , η u x = D x η − u t D x ξ1 − u x D x ξ2 , ηu xx = D x ηu x − utx D x ξ1 − u xx D x ξ2 . The prolongation formulae for the remaining extensions, defined recursively, may be found in [1]. The invariance condition (7) leads to obtain an overdetermined set of linear differential equations (determining equations) for the infinitesimals, whose integration gives the generators of Lie point symmetries admitted by the equation (3), see [32,31]. The extension of the Lie approach for DEs to FDEs was proposed in the papers [9–11,22], where a new infinitesimal ηuαt of fractional derivative, when the considered model involves a time fractional derivative, is defined as

ηuαt = D tα (η) + ξ2 D tα (u x ) − D tα (ξ2 u x )

(8)

+ D tα ( D t (ξ1 ) u ) − D tα +1 (ξ1 u ) + ξ1 D tα +1 (u )

denoting with D tα (·) the operator of Riemann-Liouville fractional differentiation with respect to t and D t (·) the classical total derivative. So that, in the context of fractional differential equations, the invariance condition (7) involves also the previous infinitesimal ηuαt , and the coefficients of the determining equations depend on all the derivatives of variable u and also on the D tα u. For finding the Lie symmetries admitted by fractional differential equations we use an algorithm implemented in MAPLE, symbolic computing platform, developed according to the Gazizov method and requiring that the FPDEs are written in terms of the R-L derivatives. The package FracSym [39,21] uses some routines of the MAPLE symmetry packages as DESOLVII and ASP, that automate the search of symmetries for FDEs. The infinitesimals characterizing the Lie symmetries of the SF–ADR equation (1) are

ξ1 = a 1 ,

ξ2 = 0 ,

where the functions

χ = χ (t , x) and ψ(t ) satisfy the constraint

β+1

∂t χ − k1 ∂x

η = χ (t , x) + ψ(t )u

χ + k2 ∂x χ + ψ  u + a1 ∂t f + (χ + u ψ)∂u f − ψ f = 0,

(9)

that links the arbitrary functions χ = χ (t , x) and ψ = ψ(t ) to f = f (t , x, u ), with a1 being an arbitrary parameter. The Lie transformation (assumed a1 = 1), involving arbitrary functions, is

X = x,



U = u (t , x)e − ψ(τ )dτ −





e − ψ(τ )dτ χ (τ , x)dτ ;

in order to integrate the equation (9) suitable assumptions on the functions ψ(t ) and

(10)

χ (t , x) are needed.

2.2. Reductions of FPDEs to FODEs In the following, two different examples of reductions of FPDE are presented, for suitable assumptions on the functions

ψ(t ) and χ (t , x) in (10). Case 1. When we assume ψ(t ) = a, a = const, through the transformation (10), we get the exact solution of the SF–ADR problem



u (t , x) = eat

 U (x) +

e −aτ χ (τ , x)dτ



,

and by integration of (9), we obtain



f (t , x, u ) = e

at

 φ(x, U ) −

e

−aτ



β+1

∂τ χ (τ , x) − k1 ∂x

  χ (τ , x) + k2 ∂x χ (τ , x) dτ

where φ(x, U ) is an arbitrary function of its arguments. When the transformation (10) and the previous form of f (t , x, u ) are inserted into the equation (1), it is reduced into the following fractional nonlinear ordinary differential equation β+1

−k1 D X

U ( X ) + k2 D X U ( X ) + φ( X , U ( X )) + aU ( X ) = 0.

(11)

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.4 (1-10)

A. Jannelli et al. / Applied Numerical Mathematics ••• (••••) •••–•••

4

The choice of the arbitrary function φ( X , U ) characterizes the solution of the equation (11) that with a suitable choice of

χ = χ (t , x) specializes the source term f (t , x, u ) and the classes of solutions. By setting φ( X , U ) = −k1 D X φ ∗ ( X ) − aU ( X ), we get the following FODE β

D X U (X) −

k2 k1

U ( X ) + φ ∗ ( X ) = 0,

β+1

β

(12)

β

being D X U ( X ) = D X ( D X U ( X )) = D X ( D X U ( X )) when U (0) = 0 (see Eq. (2.123) in [33]). In [20], we numerically solved the Eq. (12) and for some choices of function φ ∗ ( X ) we reported the comparison with the exact solutions, written in terms of the Mittag-Leffler function. The efficiency and the accuracy of the proposed procedure were shown by the comparison of the numerical solutions with the exact ones. Case 2. When we assume

χ (t , x) = 0, the transformation (10) leads to get the exact solution of the SF–ADR equation (1)



u (t , x) = U (x)e ψ(τ )dτ . In this case, we are able to integrate (9) and get 

f (t , x, u ) = −ψ(t )u + e ψ(τ )dτ φ(x, U ) where φ(x, U ) is an arbitrary function of its arguments. When we insert the previous form of f (t , x, u ) into the equation (1), it reduces to the following fractional nonlinear ordinary differential equation β+1

−k1 D X

U ( X ) + k2 D X U ( X ) + φ( X , U ( X )) = 0.

(13)

The choice of the arbitrary functions φ( X , U ) characterizes the solution of the equation (13) and specializes the source term f (t , x, u ). The suitable choice of function φ( X , U ) = −k1 D X φ ∗ ( X ) allows to get β

D X U (X) −

k2 k1

U ( X ) + φ ∗ ( X ) = 0,

β+1

β

(14)

β

being D X U ( X ) = D X ( D X U ( X )) = D X ( D X U ( X )) with U (0) = 0 (see Eq. (2.123) in [33]). In [19], the numerical solutions of (14) were compared to the exact solutions, obtained with some choices of function φ ∗ ( X ) and written in terms of the Mittag-Leffler function. The analysis of the errors and the convergence order proved the accuracy, the efficiency and the reliability of the proposed procedure. We remark that the FODEs (11) and (13), obtained by the Lie transformations of the original FPDE, are equations of β + 1 order, that for suitable assumptions on the function φ( X , U ( X )) are reduced to β order equations. In [19,20] β order FODEs have been investigated and the obtained solutions allowed to get solutions of the original SF–ADR equation with linear source terms. In this paper, we investigate the β + 1 order FODEs (11) and (13) that allow to get solutions of SF–ADR equation with nonlinear source terms defining, in this way, a wider classes of problems that we can solve by the proposed procedure. 3. Numerical approach In this section, we present the numerical approach used to solve the β + 1 order FODE (13). According to the procedure followed in [17,18], we introduce the Caputo derivative and rewrite the reduced FODE in terms of this derivative. Then, we use a numerical method in order to find approximate solutions. The obtained numerical solutions are used for evaluating the solutions of the original advection–diffusion model (1). We consider the Caputo’s approach because its main advantage is that the initial conditions take the same form as the one for integer-order differential equations, i.e., they contain the limit values of integer-order derivatives of unknown functions at the lower point. So that the initial conditions assume a physical meaning in agreement with many natural phenomena. We recall the definition of the Caputo fractional derivative of a function V ( X ) ∗

1

β

D X V (X) =

(1 − β)

X 0

1

d

( X − S )β dS

V ( S )dS

and its link with the Riemann-Liouville fractional derivative ∗

β

β

D X V ( X ) = D X ( V ( X ) − V (0)) .

(15)

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.5 (1-10)

A. Jannelli et al. / Applied Numerical Mathematics ••• (••••) •••–•••

5

In this context, we have ∗

β+1

DX

U ( X ) = ∗ D X ( D X U ( X )) = D X ( D X U ( X ) − D X U (0)) β

β

β+1

= DX

U (X) −

D X U (0)

(1 − β) X β

,

0 < β < 1,

D X U (0) being the initial condition imposed on the first derivative of the solution. Then, the equation (13) can be written as follows

−k∗1 D X

β+1

U ( X ) + k2 D X U ( X ) + φ( X , U ( X )) = 0

(16)

where we set D X U (0) = 0. Eq. (16) is a multi-term FODE, that, in a general form, is described by the following equation

α I∗ D βXI U ( X ) + · · · + α2∗ D βX2 U ( X ) + α1∗ D βX1 U ( X ) = g ( X , U ( X ))

(17)

where α1 , α2 , · · · , α I are real coefficients and β1 , β2 , · · · , β I , with β I > · · · > β2 > β1 , are the orders of the fractional derivatives. The function g ( X , U ( X )) is assumed to be nonlinear. The initial conditions related to equation are m I , with mi = βi  for i = 1, · · · , I (1 )

U ( X0) = U 0,

D X U ( X0) = U 0 ,

··· ,

m −1

DXI

(m I −1)

U ( X0) = U 0

(18)

.

In order to numerically solve the fractional initial value problem (FIVP) (17)-(18), let us introduce a mesh of points X j = X 0 + jh, for j = 0, . . . , J , with uniform spacing h, where naturally X 0 = 0 and X J = X max . We denote the numerical approximation of exact solution U ( X j ) of the problem (16) at all points by U j . A lot of methods have been developed and can be found in the specialized literature. Our main aim is to propose a numerical procedure simple to implement, with high accuracy and good stability properties. For this reason, in this paper, we focus only on the following fractional trapezoidal method

U j = T (X) −

I −1  αi i =1

αI

⎛

hβ I −βi ⎝a(β I −βi ) U 0 +

j  l =1

⎛

⎞

(β −β ) a j −Il i U l ⎠ +

⎞

 (β ) 1 a j −Il g ( X j , U j )⎠ + hβ I ⎝a(β I ) g ( X 0 , U 0 ) + j

αI

(19)

l =1

with

a (γ ) = (γ )

a0

=

( j − 1)γ +1 − j γ ( j − γ − 1) , (γ + 2) 1

(γ + 2)

,

(γ )

aj

=

( j − 1)(γ +1) − 2 j (γ +1) + ( j + 1)(γ +1) (γ + 2)

j = 1, 2, · · ·

and

T ( X ) = T m I −1 [U , X 0 ]( X ) +

m i −1 I −1  αi  ( X − X 0 )k+β I −βi (k) U ( X0) , αI (k + β I − βi + 1) i =1

k =0

where T m I −1 [U , X 0 ]( X ) is the Taylor polynomial of degree m − 1 for the function U ( X ) centered at X 0

T m I −1 [U , X 0 ]( X ) =

m −1  k =0

( X − X 0 )k (k) U ( X0) . k!

The method (19) is a generalization to the multi-term FODEs of the classical trapezoidal method, a widely used numerical scheme for solving the linear and nonlinear ordinary differential equations. It is an unconditionally stable method with convergence order O ((h)min(β+1,2) ), for h → 0, on a uniform mesh. Generally, the convergence order of the fractional trapezoidal method usually is β + 1 when 0 < β < 1, and only when the solution is sufficiently smooth or when β > 1, the expected order two is reached. The fractional trapezoidal method assures a high accuracy and has good stability property but it is implicit: at each step of the integration of the fractional ordinary differential equations by the numerical method, a nonlinear equation must be solved and therefore an algorithm for the solution of nonlinear equations is required. In this paper, we use the classical Newton’s method, one of the most popular iterative methods for solving nonlinear equations; at each iteration step it requires the evaluation of the Jacobian matrix that, often, is too expensive to compute or not easily

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.6 (1-10)

A. Jannelli et al. / Applied Numerical Mathematics ••• (••••) •••–•••

6

Fig. 1. Numerical solution U j of the IVP (21)-(22). Left frame: β = 0.1. Right frame: β = 0.8.

obtainable analytically and, then, it has to be evaluated numerically. In this paper, we do not use a built-in procedure in order to implement the Newton method and, moreover, to avoid a loss of accuracy by introducing a further source of error, the Jacobian matrix is analytically computed. We refer to the papers [8], [5] and [28] for the technical details on the proposed method. In this paper, the FIVP (17)-(18) reduces to

−k∗1 D X

β+1

U ( X ) + k2 D X U ( X ) = −φ( X , U ( X )) (1 )

U ( X0) = U 0,

D X U (0) = U 0 = 0

(20)

β with I = 2, β2 = β + 1 and β1 = 1 so that ∗ D X1 = D X . Moreover, we set −φ( X , U ( X )).

α2 = −k1 and α1 = k2 and, finally g ( X , U ( X )) =

4. Applications In this Section, by using the fractional trapezoidal method, we numerically solve the FIVP (20) and, then, we use the computed solution in order to derive an approximate solution of the original problem by means of the transformation (10). We find numerical solutions of the models under study for arbitrary values of the fractional order β , for which the exact solutions are unknown. In the examples, the function φ( X , U ( X )) is chosen in such a way that the exact solutions of the model (1) with β = 1 have a physical meaning. Then, we compare them with the numerical solutions with β approaching 1 with the aim of showing the applicability and efficiency of the proposed numerical procedure described in previous Section. All numerical simulations are performed on Intel Core 7 by using Matlab software. Example 1. By assuming

 

φ( X , U ( X )) = ak2 1 − U ( X ) + b ,

with a and b arbitrary constants, and replacing in the Eq. (16) we obtain the following FODE

−k∗1 D X

β+1



U ( X ) + k2 D X U ( X ) + ak2 1 −



U ( X ) + b = 0,

(21)

subject to the initial conditions

U (0) = U 0

D X U (0) = 0.

The exact solution of the initial value problem (IVP) (21)-(22) when β = 1 and a = k2

U ( X ) = (e 2k1

X

− 1)2 − b.

(22) k2 2k1

is

(23)

In Fig. 1, we report the numerical solution of the IVP with U 0 = −b obtained for β = 0.1 and β = 0.8 and the exact solution of the IVP with β = 1. The results are performed on the interval [0, 5] with J = 100 grid points by setting k1 = 0.5, k2 = −2 and b = −0.1. For both choices of β , the average number of Newton’s iterations is equal to 3 at the first computation steps. After, the average number of Newton’s iterations reduces to 2 and 1 for β = 0.1 and β = 0.8, respectively.

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.7 (1-10)

A. Jannelli et al. / Applied Numerical Mathematics ••• (••••) •••–•••

7

Fig. 2. Numerical solution unj of the model (1). Left frame: β = 0.1. Right frame: β = 0.8.

Now, by using the transformation (10) we get 

u (t , x) = U (x)e ψ(t )dt

(24)

that is the solution of the original SF–ADR model whose source term, obtained by integration of the (9), is

  k22  ψ(t )dt f (t , x, u ) = −ψ(t )u + e (1 − ue − ψ(t )dt + b). 2k1

In the following, we report two possible choices of the arbitrary function ψ(t , x) that allow to characterize the solution u (t , x) and specialize the source term f (t , x, u ).

• If we assume ψ(t ) = c, with c arbitrary constant, the transformation reads u (t , x) = U (x)ect and the source term assumes the form

f (t , x, u ) = −cu +

k22 2k1

ect (1 −

ue −ct + b).

The above nonlinear function involves a linear term and a nonlinear mth-order reaction term where m = 1/2, that describes, being 0 < m < 1, absorption processes like chemical reaction ones. The numerical solution unj , approximation of the exact solution u (t n , x j ) of SF–ADR equation, is reported in Fig. 2. For the computations, we set c = −2 and consider a domain [0, 5] × [0, 5], discretized in time with uniformly spaced mesh of N + 1 points t n , for n = 0, . . . , N with N = 100. 1 • If we assume ψ(t ) = β− , the transformation reads t +t c

u (t , x) = U (x)(t + t c )β−1

(25)

and the source term of the original SF–ADR model is given by

f (t , x, u ) = −

 

k2 β −1 u + 2 (t + t c )β−1 1 − (t + t c )1−β u + b . t + tc 2k1

In Fig. 3, the numerical solution with t c = 2 is shown. When β = 1, the exact solution u (t , x) = U (x) is given by (23) and the source term reads

f (t , x, u ) =

k22 2k1

(1 −

u + b).

In this case the source term characterizes a model belonging to the classes of models studied by Herrera, Minzoni and Ondarza [14] describing annihilation processes. Example 2. By assuming

  φ( X , U ) = aU ( X ) e (β−1) X − U ( X )

with a arbitrary constant, we obtain the following nonlinear FODE

−k∗1 D X

β+1





U ( X ) + k2 D X U ( X ) + aU ( X ) e (β−1) X − U ( X ) = 0

(26)

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.8 (1-10)

A. Jannelli et al. / Applied Numerical Mathematics ••• (••••) •••–•••

8

Fig. 3. Numerical solution unj of the model (1). Left frame: β = 0.1. Right frame: β = 0.8.

Fig. 4. Numerical solution U j of the IVP (26)-(27). Top: Left frame: β = 0.1. Right frame: β = 0.8. Bottom: Left frame: β = 0.98. Right frame: β = 0.999.

subject to the initial conditions

U (0) = U 0 ,

D X U (0) = 0.

(27)

In Fig. 4, we report the numerical solutions obtained for β = 0.1, 0.8, 0.98 and 0.999. The results are performed on the interval [0, 40] with J = 100 grid points by setting k1 = 1, k2 = −1.2 and a = 5. For all choices of fractional order β , the average number of Newton’s iterations is equal to 3 at the first computation steps; at the end of the process, the average number of Newton’s iterations reduces to 2.

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.9 (1-10)

A. Jannelli et al. / Applied Numerical Mathematics ••• (••••) •••–•••

9

Fig. 5. Numerical solution unj of the model (1). Top: Left frame: β = 0.1. Right frame: β = 0.8. Bottom: Left frame: β = 0.98. Right frame: β = 0.999.

By means of the transformation (10) we get the following solution of the original SF–ADR model (1) 

u (t , x) = U (x)e ψ(τ )dτ whose source term, obtained by integration of the (9), is given by







f (t , x, u ) = −u ψ(t ) − a(e (β−1)x − ue − ψ(τ )dτ ) . β−1

By assuming ψ(t ) = t +t , the transformation reads c

u (t , x) = U (x)(t + t c )β−1

(28)

and the source term is given by



 β −1 (β−1)x 1−β f (t , x, u ) = −u − a(e − (t + t c ) u) . t + tc

(29)

In Fig. 5, the numerical solution with t c = .2 is shown for different values of β . We point out that the behavior of the numerical solution of the FPDE, computed for β → 1, (see Fig. 5), is in agreement with the solution of the integer order ADR equation with a logistic-type nonlinear source term

∂t u (t , x) − k1 ∂xx u (t , x) + k2 ∂x u (t , x) + au (t , x)(1 − u (t , x)) = 0, obtained by the original FPDE (1) with source term given by (29) when β = 1 (see Sect. 1.1.5 and 1.1.1. of [34]) 5. Conclusion In this paper, we investigate the SF–ADR models, in terms of the Riemann-Liouville derivative, with nonlinear source terms. We apply an approach that combines the fractional Lie symmetry analysis with the numerical method and allows to solve the original model by means of the resolution of a β + 1 order FODE. The reported numerical results, obtained by solving the FODE by the implicit trapezoidal method, confirm the efficiency and the applicability of the proposed procedure to a wide class of FPDEs.

JID:APNUM AID:3756 /FLA

[m3G; v1.261; Prn:27/01/2020; 12:52] P.10 (1-10)

A. Jannelli et al. / Applied Numerical Mathematics ••• (••••) •••–•••

10

Acknowledgements A. J. is a member of the INdAM Research group GNCS. M.R. and M.P.S. are members of the INdAM Research group GNFM. References [1] G.W. Bluman, S. Kumei, Symmetries and Differential Equations, Springer-Verlag Inc., New York, 1989. [2] J. Cao, C. Li, Y. Chen, High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations (II), Fract. Calc. Appl. Anal. 18 (2015) 735–761. [3] J.F. Cheng, Y.M. Chu, Solution to the linear fractional differential equation using Adomian decomposition method, Math. Probl. Eng. (2011) 587068, https://doi.org/10.1155/2011/587068. [4] V. Daftardar-Geji, H. Jafari, Adomian decomposition: a tool for solving a system of fractional differential equations, J. Math. Anal. Appl. 301 (2005) 508–518. [5] K. Diethelm, The Analysis of Fractional Differential Equations, Springer-Verlag, Berlin, 2004. [6] E.H. Doha, A.H. Bhrawy, D. Baleanu, S.S. Ezz-Eldien, Numerical schemes with high spatial accuracy for a variable-order anomalous subdiffusion equation, Appl. Math. Comput. 219 (2013) 8042–8056. [7] R. Fazio, A. Jannelli, S. Agreste, A finite difference method on non-uniform meshes for time-fractional advection-diffusion equations with a source term, Appl. Sci. 8 (960) (2018). [8] R. Garrappa, Numerical solution of fractional differential equations: a survey and a software tutorial, Mathematics 6 (16) (2018). [9] R.K. Gazizov, A.A. Kasatkin, S.Y. Lukashchuk, Continuous transformation groups of fractional differential equations, Vestn. USATU 9 (2007) 125–135. [10] R.K. Gazizov, A.A. Kasatkin, S.Y. Lukashchuk, Symmetry properties of fractional diffusion equations, Phys. Scr. T 136 (2009) 014016. [11] R.K. Gazizov, A.A. Kasatkin, S.Y. Lukashchuk, Group-invariant solutions of fractional differential equations, in: Nonlinear Science and Complexity, 2011, pp. 51–59. [12] J.H. He, A coupling method of a homotopy technique and a perturbation technique for nonlinear problems, Int. J. Non-Linear Mech. 35 (2000) 37–43. [13] H. Hejazi, T. Moroney, F. Liu, Stability and convergence of a finite volume method for the space fractional advection-dispersion equation, J. Comput. Appl. Math. 255 (2014) 684–697. [14] J.J.E. Herrera, A. Minzoni, R. Ondarza, Physica D 57 (1992) 249–266. [15] N.H. Ibragimov, Transformation Groups Applied to Mathematical Physics, Reidel, Dordrecht, 1985. [16] N.H. Ibragimov, CRC Handbook of Lie Group Analysis of Differential Equations, CRC Press, Boca Raton, 1994–1995–1996. [17] A. Jannelli, M. Ruggieri, M.P. Speciale, Analytical and numerical solutions of fractional type advection-diffusion equation, AIP Conf. Proc. 1863 (1) (2017) 530005. [18] A. Jannelli, M. Ruggieri, M.P. Speciale, Exact and numerical solutions of time-fractional advection–diffusion equation with a nonlinear source term by means of the Lie symmetries, Nonlinear Dyn. 92 (2018) 543–555. [19] A. Jannelli, M. Ruggieri, M.P. Speciale, Analytical and numerical solutions of time and space fractional advection–diffusion–reaction equation, Commun. Nonlinear Sci. Numer. Simul. 70 (2019) 89–101. [20] A. Jannelli, M. Ruggieri, M.P. Speciale, Numerical solutions of space fractional advection–diffusion equation with a source term, AIP Conf. Proc. 2116 (2019) 280007. [21] G.F. Jefferson, J. Carminati, ASP: automated symbolic computation of approximate symmetries of differential equations, Comput. Phys. Commun. 184 (2013) 1045–1063. [22] R.A. Leo, G. Sicuro, P. Tempesta, A theorem on the existence of symmetries of fractional PDEs, C. R. Math. 352 (3) (2014) 219–222. [23] C. Li, J. Cao, H. Li, High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations (III), J. Comput. Appl. Math. 299 (2016) 159–175. [24] C. Li, R. Wu, H. Ding, High-order approximation to Caputo derivatives and Caputo-type advection-diffusion equations (I), Commun. Appl. Ind. Math. 6 (2014) e-536. [25] Y. Lin, X. Li, C. Xu, Finite difference/spectral approximations for the fractional cable equation, Math. Comput. 80 (2011) 1369–1396. [26] Y. Lin, C. Xu, Finite difference/spectral approximations for the time-fractional diffusion equation, J. Comput. Phys. 225 (2007) 1533–1552. [27] F. Liu, F.P. Zhuang, V. Anh, I. Turner, K. Burrage, Stability and convergence of the difference methods for the space-time fractional advection-diffusion equation, Appl. Math. Comput. 191 (2007) 12–20. [28] C. Lubich, Discretized fractional calculus, SIAM J. Numer. Anal. 17 (3) (1986) 704–719. [29] K.S. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley and Sons, 1993. [30] S. Momani, Z. Odibat, Analytical approach to linear fractional partial differential equations arising in fluid mechanics, Phys. Lett. A 355 (4–5) (2006) 271–279. [31] P.J. Olver, Applications of Lie Groups to Differential Equations, Springer, New York, 1986. [32] L.V. Ovsiannikov, Group Analysis of Differential Equations, Academic Press, New York, 1982. [33] I. Podlubny, Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, Some Methods of Their Solution and Some of Their Applications, Academic Press, San Diego, 1999. [34] A.D. Polyanin, V.F. Zaitsev, Handbook of Nonlinear Partial Differential Equations, Chapman and Hall/CRC, 2004. [35] J. Ren, Z. Sun, X. Zhao, Compact difference scheme for the fractional sub-diffusion equation with Neumann boundary conditions, J. Comput. Phys. 232 (2013) 456–467. [36] M. Ruggieri, M.P. Speciale, On the construction of conservation laws: a mixed approach, J. Math. Phys. 58 (2017) 023510. [37] J. Saberi-Nadjafi, A. Ghorbani, He’s homotopy perturbation method: an effective tool for solving nonlinear integral and integro-differential equations, Comput. Math. Appl. 58 (2009) 1345–1351. [38] S. Samko, A.A. Kilbas, O. Marichev, Fractional Integrals and Derivatives, Taylor and Francis, 1993. [39] K.T. Vu, G.F. Jefferson, J. Carminati, Finding generalized symmetries of differential equations using the MAPLE package DESOLVII, Comput. Phys. Commun. 183 (2012) 1044–1054. [40] N. Whang, W. Deng, Y. Wu, Finite difference/element method for a two-dimensional modified fractional diffusion equation, Adv. Appl. Math. Mech. 4 (2012) 496–518. [41] F. Zeng, C. Li, F. Liu, I. Turner, The use of finite difference/element approaches for solving the time-fractional subdiffusion equation, SIAM J. Sci. Comput. 35 (2013) 2976–3000. [42] Y. Zhang, Z. Sun, H. Liao, Finite difference methods for the time fractional advection diffusion equation on non-uniform meshes, Comput. Phys. 265 (2014) 195–210. [43] P. Zhuang, F. Liu, Implicit difference approximation for the time fractional diffusion equation, J. Appl. Math. Comput. 22 (2006) 87–89.