- Email: [email protected]
A Green’s function approach for the numerical solution of a class of fractional reaction–diffusion equations
Accepted Manuscript A Green’s function approach for the numerical solution of a class of fractional reaction-diffusion equations Eliseo Hernandez-Martinez, Francisco Valdes-Parada, Jose Alvarez-Ramirez, Hector Puebla, Epifanio Morales-Zarate PII: DOI: Reference:
S0378-4754(15)00192-5 http://dx.doi.org/10.1016/j.matcom.2015.09.004 MATCOM 4236
To appear in:
Mathematics and Computers in Simulation
Received date: 30 July 2014 Revised date: 13 September 2015 Accepted date: 14 September 2015 Please cite this article as: E. Hernandez-Martinez, F. Valdes-Parada, J. Alvarez-Ramirez, H. Puebla, E. Morales-Zarate, A Green’s function approach for the numerical solution of a class of fractional reaction-diffusion equations, Math. Comput. Simulation (2015), http://dx.doi.org/10.1016/j.matcom.2015.09.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
(Revised Manuscript No. MATCOM-D-14-00272R1) A Green’s function approach for the numerical solution of a class of fractional reaction-diffusion equations Eliseo Hernandez-Martinez1∗, Francisco Valdes-Parada2 , Jose Alvarez-Ramirez2 Hector Puebla3 and Epifanio Morales-Zarate1 1 Facultad
de Ciencias Qu´ımicas
Universidad Veracruzana, Xalapa Veracruz, M´exico 2 Departamento
de Ingenier´ıa de Procesos e Hidraulica
Universidad Aut´onoma Metropolitana-Iztapalapa, M´exico D.F., M´exico 3 Departamento
de Energ´ıa,
Universidad Aut´onoma Metropolitana-Azcapotzalco, M´exico D.F., M´exico ∗ Email:
[email protected]
September 13, 2015
Abstract Reaction-diffusion equations with spatial fractional derivatives are increasingly used in different science and engineering fields to describe spatial patterns arising from the interaction of chemical or biochemical reactions and anomalous diffusive transport mechanisms. Most numerical schemes to solve fractional reaction-diffusion equations use finite difference schemes based on the Gr¨ unwaldLetnikov formula. This work introduces a new systematic approach based on Green’s function formulations to obtain numerical schemes for fractional reaction-diffusion equations. The idea is to pose an integral formulation of the equation in terms of the underlying Green’s function of the fractional operator to subsequently use numerical quadrature to obtain a set of ordinary differential equations. To illustrate the numerical accuracy of the method, dynamic and steady-state situations are considered and compared with analytical and numerical solutions via Gr¨ unwald finite differences schemes. Numerical simulations show that the scheme proposed improves the performance ∗
Corresponding author
1
and convergence of traditional finite differences schemes based on Gr¨ unwald formula. Keywords. Anomalous diffusion; reaction-diffusion equations; Green’s function; numerical solution.
1
Introduction
Reaction-diffusion systems have been used in different science and engineering fields. However, recent research indicates that the classical diffusion equation is inadequate to model many real situations. For example, in processes where the particles spread at a faster rate than that predicted by Fickian diffusion and can present asymmetry. This phenomenon is commonly known as anomalous diffusion [23, 27]. Fractional reaction-diffusion equations can describe a large number of physical processes controlled by anomalous diffusion. Some example include, solute transport [31], transport in plasma turbulence [7], in the pattern formation [13, 12], transport at the Earth surface [30], diffusion through porous materials [3, 35]. Lenzi et al. [18] showed that the anomalous diffusion is a suitable framework for modeling transport in catalytic particles. Recently, Zheng et al. [38] considered fractal theory for modeling diffusion in porous media, finding good agreement between experimental data and theoretical predictions. On the other hand, Das [6] computed approximate analytical solutions for fractional diffusion equations that can be used in the modeling of reaction-diffusion in catalytic systems. Liu et al. [21] used a fractional differential equation to describe anomalous diffusion in a mobile/immobile transport model. Friess et al. [11] showed that non-Fickian diffusion phenomenon should be considered for describing transport of some alcohol vapors in teflon material. Different numerical methods have been used to solve both linear and nonlinear fractional reactiondiffusion equations such as, standard discretization of the fractional-in-space operator [16], Fourier methods [4], homotopy methods [17], finite element [37] and Gr¨ unwald fractional difference (GrFD). Due to its simplicity of computational implementation, the latter is the most commonly used [22, 8, 25, 5]. However, it is well-known that the discretization of fractional derivatives by means of GrFD is of first-order accurate O(h) and it may lead to instabilities [33]. For this reason, several works are focused in the derivation of stable and high-order schemes [9, 10, 36]. In [19], Li and Zeng presented a review about GrFD schemes where they discussed their application to different types of fractional differential equations. Due to the extensive applications of the fractional models, it is necessary to dispose of stable numerical methodologies with high approximation order. In this sense, recent studies have demonstrated that numerical schemes based on Green’s functions are an alternative for the solving of reaction2
diffusion equations [34]. In fact, systematic derivation of nonstandard FD schemes can be obtained by means of Green’s function formulations [1, 15]. Whereby, in this work a new systematic approach for the discretization of fractional reaction-diffusion equations based on Green’s function formulations of the fractional operator is proposed. The numerical solution of proposed scheme is based on method of lines (MOL), where a Runge-Kutta method is used for solving the ordinary differential equations (ODE) system obtained. The paper is organized as follows, first we present the integral formulation of the fractional differential equation [2, 32]. Subsequently, a quadrature rule is used to obtain a set of ordinary differential equations. In addition, depending on the discretization of the integral, a numerical scheme with equivalent structure to finite differences based on Gr¨ unwald formula can be obtained. To illustrate the numerical accuracy of the method, the results are compared with the analytical and numerical solution via GrFD.
2
Preliminaries
For convenience, we start our derivations by presenting some definitions and properties of the fractional integral and differential operators [2, 32]. Later, we describe the derivation of integral formulation for solving fractional reaction-diffusion systems. For the study of fractional differential equations, different definitions of the derivatives and integrals have been proposed [29, 28, 26]. The integral operator of fractional of order α > 0 of a function u(x) is defined as [2], Iaα u(x)
1 = Γ(α)
Z
x a
(x − z)α−1 u(x)dz
(1)
where Iaα is the so-called Riemann-Liouville fractional integral operator of order α. On the other hand, the fractional differential operator for a function u(x) is defined as [32], Daα u(x)
1 = Γ(n − α)
Z
x a
(x − z)n−α−1
dn u(z) dz dz n
(2)
where Daα is the Caputo derivative of order α and n is an integer such that n − 1 < α ≤ n. Note that
when α = 2 one recovers the classical diffusion operator. It is assumed that u(x) is a continuously
differentiable function and piecewise smooth, in order for the following considerations to be true [2, 32], • Lemma 1. For α > 1, if we assume that u(x) ∈ C(a, b) ∩ L(a, b), then the Caputo homogeneous differential equation (Daα u(x) = 0) has the following solution
u(x) = c0 + c1 x + c2 x2 + . . . + cn−1 xn−1 3
• Remark 1. For α > 1, if we assume u(x) ∈ C(a, b), then Daα Iaα u(x) = u(x) almost everywhere. • Lemma 2. If α > 1 and assuming that u(x) ∈ C(a, b) ∩ L(a, b), then Iaα Daα u(x) = u(x) + c1 + c2 x + . . . + cn−1 xn−1 for some ci ∈ R, i = 1, . . . , n, where n is the smallest integer greater than or equal to α. Anomalous diffusion coupled with chemical reactions may be described by fractional diffusion operators such as, ∂u(x, t) = Daα u(x, t) − f (u(x, t)), ∂t
x ∈ [a, b]
(3)
with initial condition as u(x, 0) = φ(x)
(4)
and the boundary conditions ∂u(a, t) + βa u(a, t) + γa = 0 ∂x ∂u(b, t) + βb u(b, t) + γb = 0 αb ∂x
αa
(5)
where f (u(x, t)) is a continuously differentiable function of its arguments and α ∈ (1, 2] is the order
of the derivative. Eq. (3) arises from a mass balance and a fractional Fick’s law where the mass flux is proportional to the fractional gradient −∇α−1 u(x, t) [35]. Note that α = 2 corresponds to the classical diffusion equation. The case α ∈ (1, 2) models a superdiffusive flow in which a cloud of diffusing particles spreads at a faster rate than the classical diffusion model. For developing the integral formulation, let us rewrite Eq. (3) as Daα u(x, t) = Ψ(x, t) where Ψ =
∂u(x,t) ∂t
x ∈ [a, b]
(6)
+ f (u(x, t)). For Lemma 2 with n = 2, if one multiplies Eq. (6) by the Riemann-
Liouville integral operator, one obtains Iaα Daα u(x, t) = u(x, t) + c0 + c1 x = Iaα Ψ(u, t) or, equivalently
1 u(x, t) = Γ(α)
Z
a
x
(x − z)α−1 Ψ(x, t)dz − c0 − c1 x
Evaluating the boundary condition given by Eqs. (5) into Eq. (8) it results that γa αa + βa βa Z b Z b αb (b − z)α−1 (b − z)α−2 αb αa 1 Ψ(z, t)dz + Ψ(z, t)dz + − b − a + αβbb − αβaa βb a Γ(α − 1) Γ(α) βb βa a
c0 = −c1 a − c1 c1 =
4
(7)
(8)
leading to u(x, t) =
a−x+
αa βa
αb βb
Z
b
(b − z)α−2 Ψ(z, t)dz + Γ(α − 1)
b − a + αβbb − αβaa a Z x 1 γa + (x − z)α−1 Ψ(z, t)dz − Γ(α) a βa
Z
b a
(b − z)α−1 αb αa Ψ(z, t)dz + − Γ(α) βb βa
(9)
After some algebraic manipulations, Eq. (9) takes the form Z b γb γa GC (z, x)Ψ(z, t)dz − σ(x) + u(x, t) = (σ(x) + 1) βa βb a where σ(x) =
x−a−
b−a+
αb βb
Γ(α)
(10)
αa βa
−
αa βa
and GC (z, x) is the Caputo fractional Green’s function, − αb (b−z)α−2 σ(x) − (b−z)α−1 σ(x) + βb Γ(α−1) Γ(α) GC (z, x) = − αb (b−z)α−2 σ(x) − (b−z)α−1 σ(x) βb Γ(α−1)
(x−z)α−1 Γ(α)
if z < x if z ≥ x
(11)
Eq. (11) is the general Green’s function, which, as a whole with Eq. (10), denote the Caputo
Green’s function formulation (Caputo-GFF) for fractional reaction-diffusion equation given by Eq. (1). In order to observe the effect of anomalous transport on the Green’s function, we consider Dirichlet-type boundary conditions (αa = αb = 0, βa = βb = 1), for x = 0.5 and z ∈ [0, 1], Figure
1a shows the fractional Green’s function distribution in the domain [0, 1], where we can observe that
for 1 < α < 2 the Green’s function is not symmetrical (anomalous diffusion) whereas for α = 2.0 the symmetrical profile is recovered (Fickian diffusion). Green’s functions tend to zero only at z = 1.0, whereas at z = 0.0 the Green’s functions are inversely proportional to α, this behavior is compensated with the contribution of the boundary conditions. For mixed-type boundary conditions (αb = βa = 0, αa = βb = 1), Figure 1b shows that the Green’s function distribution is dissipated at the Neumanntype boundary, which increase as α decreases. On the other hand, unlike traditional methods for the numerical solution of boundary value problems (e.g., finite differences), Caputo-GFF can provide an interpretation of how to carry out the transport and reaction processes. Some comments regarding the structure of Green’s function approach derived above are in order: i) Caputo Green’s function GC (x, z) can be reduced to the standard Green’s functions as α → 2 and according to the parameters αi and βi , it is possible to obtain Green’s function for different
5
types of boundary conditions (e.g., Dirichlet, Neumann, Robin). For instance, if we consider α = 2 the Green’s function is G(z, x) =
1 b − a + αβbb −
αa βa
z − a − αa x−b− βa x − a − αa z−b− βa
αb βb αb βb
if z < x if z ≥ x
(12)
obtaining the integral formulation for the classical reaction-diffusion system with general-type boundary conditions. This result was reported in [14], indicating that Eq. (11) is the general Green’s function for reaction-diffusion systems with general-type boundary conditions in a slab. ii) In the original boundary value problem the diffusive transport is described by the differential operator, while in the integral formulation the transport phenomena is described through the Green’s function. iii) The proposed scheme identifies that the profile u(x, t) consists of the interaction of the phenomena involved (accumulation, reaction and boundary conditions) and how they are distributed by Green’s function. iv) The first two terms on the right side of equation (10) describe the influence of inhomogeneous boundary conditions. v) The integral contains two contributions: the memory of the system that indicates how the concentration profiles change over the time and how the fractional Green’s function distributes the source in the domain (i.e., Green’s function describes the anomalous diffusion phenomena in time and space). When
∂u(x,t) ∂t
= 0, the profile of u(x, t) does not depend on time, then Eq.
(10) describes the profiles of Eq. (1) in steady-state. vi) The proposed scheme is a systematic methodology that can be applied to different definitions of fractional derivatives. For example, the Riemann-Liouville derivative is given by Daα u(x) = Rx 1 dn α+1−n u(z)dz and its properties described in Lemmas 1 and 2 are u(x) = Γ(n−α) dxn a (x − z) c1 xα−1 +c2 xα−2 +. . .+cn−1 xα−n−1 and Iaα Daα u(x) = u(x)+c1 xα−1 +c2 xα−2 +. . .+cN −1 xα−N −1 ,
respectively [2]. Following the procedure described in Eqs. (6)-(11), one can obtain that the Riemann-Liouville Green’s function formulation (RL-GFF) is given by Z b xα−1 γa γb η4 GRL (z, x)Ψ(z, t)dz σ(x) + − σ(x) + u(x, t) = η2 η2 βa βb a η1 xα−1 −η2 xα−2 , η1 = αβaa (α − 2)aα−3 + aα−2 , η2 η1 η4 −η2 η3 αb α−3 +bα−2 , η = αb (α−1)bα−2 +bα−1 and G 4 RL (z, x) is βb (α−2)b βb
where σ(x) =
6
=
αa βa (α
(13)
− 1)aα−2 + aα−1 , η3 =
the Riemann-Liouville fractional
Green’s function,
− αb (b−z)α−2 σ(x) − βb Γ(α−1) GRL (z, x) = − αb (b−z)α−2 σ(x) − βb Γ(α−1)
3
(b−z)α−1 Γ(α) σ(x) (b−z)α−1 Γ(α) σ(x)
+
(x−z)α−1 Γ(α)
if z < x if z ≥ x
(14)
Numerical scheme
In this section, the numerical methodologies, Gr¨ unwald finite differences scheme and Green’s function formulations, for solving fractional reaction-diffusion models given by Eqs. (3)-(5) are presented. For this, we consider ∆x = h > 0 as the spatial grid size, ∆x = (b − a)/ (N + 1), with xi = a + i∆x and
ui (t) = u(xi , t), for i = 1, 2, ..., N .
3.1
Gr¨ unwald finite differences scheme
Finite difference approximations are commonly used to numerically compute the solution of the fractional boundary value problem given by Eq. (3). So, in order to compare the proposed methodology, we used a spatial discretization of α-order fractional derivative, by mean of right-shifted Gr¨ unwald formula [24], such as
Daα ui (t) ≈ δx,a ui (t) = where the weights gα,k are given by gα,k =
i+1 1 X gα,k ui−k+1 (t) hα
(15)
k=0
Γ(k−1) Γ(−α)Γ(k+1) .
So, if the shifted Gr¨ unwald estimates are
substituted in Eq. (3), it is obtained the following i+1 1 X dui (t) gα,k ui−k+1 (t) − f (ui (t)), = α dt h
i = 1, 2, . . . , N
(16)
k=0
Note that for the classical diffusion equation, (α = 2), the weights g2,k are g2,0 = g2,2 = 1, g2,1 = −2
and for k > 2 g2,k = 0. Eq. (16) is a set of N differential equations that can be solved numerically by means of numerical integration methods to obtain the approximate solution of u(x, t).
3.2
Green’s function formulation
For this case, we rewrite the Caputo-GFF of u(x, t) given by (10) as Z b GC (z, x)Ψ(z, t)dz u(x, t) = B(x, t) +
(17)
a
where the boundary contribution B(x, t) is given by
B(x, t) = (σ(x) + 1) 7
γb γa − σ(x) βa βb
(18)
and the Green’s function G(z, x) is given by either Eq. (11). The proposed numerical scheme is obtained by forcing to zero the residual of Eq. (17) at the grid points x = xi . That is, Z b GC (z, xi )Ψ(z, t)dz, i = 1, ..., N ui (t) = B(xi , t) +
(19)
a
In the second step, the integral in Eq. (19) is approximated by a numerical quadrature. For instance, the composite trapezoidal rule gives N X ui (t) = ri (t) + h GC (zj , xi )Ψj (t), i = 1, ..., N
(20)
j=1
where ri (t) = B(xi , t) +
h [GC (a, xi )Ψa (t) + GC (b, xi )Ψb (t)])] 2
(21)
Let u = (u1 , u2 , ..., uN )T , r = (r1 , r2 , ..., rN )T and Ψ = (Ψ1 , Ψ2 , ..., ΨN )T . Then, Eq. (21) can be rewritten as u(t) = r(t) + hGΨ(t)
(22)
where the components of the N × N -dimensional matrix G are gij = GC (xj , xi ), i, j = 1, ..., N .
Substituting Ψ(t) =
du(t) dt
+ f (u(t)), we obtain the following set of N initial-value differential equations du(t) (u(t) − r(t)) = G−1 − f (x, u(t)) dt h
(23)
with initial conditions u(t0 ) = (φ(x1 ), φ(x2 ), ...., φ(xN ))T . Here, one can use a standard package for integration of initial value problems to obtain the approximate solution u(t). The MOL technique is based on the discretization of the spatial operator (using schemes discussed above), leading to a set of ODE. So, the system reduced is integrated in time using an ODE code such as, Runge-Kutta methods. The advantage of the MOL technique is that temporal accuracy can be specified by the users, therefore the error checking, robustness, order selection and time-step adaptive features available in sophisticated ODE codes can be applied to the time integration of the partial differential equations. For example, Liu et al. [20], proposed the MOL for solving of fractional Fokker-Planck equation, where backward differentiation formulas were used. Their results showed a good agreement between the numerical solution and analytic solution. Other alternative is to perform a full discretization that includes the temporal derivative. In this way, let δ = ∆t be a time step, tk = t0 + kδ and uki = ui (tk ). For example, assuming an implicit Euler discretization of the time derivative in Eq. (22); namely, du(tk ) dt
≈
uk −uk−1 . δ
This leads to the following expression: h i h k k u =r + G uk − uk−1 + δf (x, uk ) , k = 1, 2, .... δ 8
(24)
which, for each time tk , can be solved for the vector uk by means of, e.g., fixed-point or NewtonRaphson iterations. On the other hand, considering the integral formulation given by Eq. (10) and following the ideas of the Gr¨ unwald formula, one can propose a backward scheme with an equidistant spatial grid size, obtaining ui (t) =
xi − xi+1 −
xi+1 − a +
with GC (z, x) =
αb βb
αb βb
−
xi − a − αβaa γa − βa xi+1 − a + αβbb −
αa βa
a xi −a− α β
a α xi+1 −a+ βb b αa xi −a− β a α xi+1 −a+ βb b
h h
α−2
−z) − αβbb (xi+1 Γ(α−1)
−
α−2
− αβbb (x−i+1−z) Γ(α−1)
αa βa
γb + βb
(xi+1 −z)α−1 Γ(α)
−
Z
i
xi+1
+ i α−1
(xi+1 −z) Γ(α)
GC (z, x)Ψ(z, t)dz
(25)
a
(xi −z)α−1 Γ(α)
si z < x si z ≥ x
(26)
Now, one can perform the discretization of the integral in (25) by means of numerical approximations. Using the trapezoidal rule to approximate the integral and substituting Ψ(zj , t) =
duj (t) dt
+
f (uj (t)), one obtains the following nonlinear differential system N X duj (t) + f (uj (t)) , i = 1, ..., N ui (t) = ri (t) + h GC (zj , xi ) dt
(27)
j=1
where
αb βb
xi − a − αβaa γa γb ri (t) = − xi+1 − a + αβbb − αβaa βa xi+1 − a + αβbb − αβaa βb h dub (t) h dua (t) + f (ua (t)) GC (a, xi ) + + f (ub (t)) GC (b, xi ) + 2 dt 2 dt xi − xi+1 −
Let u=(u1 , . . . , uN )T , r=(r1 , . . . , rN )T and f (u(t)) = (f (u1 (t)), f (u2 (t)), ..., f (uN (t)))T . Then, Eq. (27) can be rewritten as
du(t) u(t) = r(t) + hG + f(u(t)) dt
(28)
where G is a lower triangular matrix of N xN with components di,j = GC (xj , xi ), i, j = 1, . . . , N . Eq. (28) can be presented as Eq. (23) and can be solved using numerical methods for solving of ODEs. Notice that the scheme (28) is equivalent to the scheme (16), where, in order to calculate ui , both schemes consider the values of [ua , u1 , . . . , ui , ub ] and the weights, gα,k /hα in the GrFD scheme are replaced for G−1 /h in the Green’s function.
4
Numerical examples
In this section we explore the ability of the proposed scheme to produce accurate solutions for fractional reaction-diffusion equations. To this end, we consider the integral formulation based on Caputo and 9
Riemann-Liouville definitions, given by Eqs. (10)-(11) and (13)-(14), respectively. The ODE systems was integrated by using Runge-Kutta-Fehlberg subroutine from Matlab package and the numerical results are compared with a classical GrFD schemes and the corresponding analytical solution. As an index for measuring approximation errors, the absolute maximum error Em was considered, which is defined as Em = max |ue (xi ) − uap (xi )|
(29)
where uap (xi ) is the numerical approximation and ue (xi ) is the analytical solution.
4.1
Steady state simulations
Consider the steady state fractional equation: dα u(x) = f (u(x)) dxα
x ∈ [0, 1]
(30)
To illustrate the effect of the fractional order (α) in the numerical accuracy of the proposed schemes, we study three cases: α = 1.8, α = 1.5 and α = 1.2. 4.1.1
Case 1
Consider Eq. (30) with α = 1.8, D(x) =
Γ(2.2) 2.8 6 x ,
f (x) = u(x) + x3 (x − 1) D(x)−1 and Dirichlet-
type boundary conditions, u(0) = 0 and u(1) = 1. The exact solution for this problem is u(x) = x3 . For 100 nodes, Figure 2a shows the profile u(x) obtained from Green’s function formulation using Caputo and Riemann-Liouville definitions and traditional Gr¨ unwald finite differences. The numerical solutions compare well (Em < 10−2 ) with the exact analytic solution to the fractional differential equation. Figure 2b shows the absolute maximum error (Em ) as function of the number of nodes (N ), the approximation errors are smaller for the Green’s function schemes than for Gr¨ unwald finite differences. Notice that, for GrFD the approximation order is O (h), whereas the for Caputo-GFF is O (h1.71 ) and Riemann-Liouville-GFF scheme is O (h1.76 ), indicating better numerical accuracy that GrFD. 4.1.2
Case 2
Let us fix now α = 1.5, f (x) =
128 √ x3.5 7 π
−
64 √ x2.5 5 π
with non-homogeneous Neumann-type boundary
conditions (du(0)/du = 0 and u(1) = 2) and exact solution as u(x) = x4 (x + 1). The numerical and
10
analytical profiles u(x) are shown in Figure 3a. The approximation errors are shown in Figure 3b; one can observe that the approximation order calculated with GrFD converges to O (h), whereas for Green’s function formulations it converges to ∼ O (h1.5 ). As in the previous case the minor Em is for
Caputo-GFF scheme. However, the approximation order for RL-GFF and Caputo-GFF schemes are the same. 4.1.3
Case 3
Consider Eq. (30) with α = 1.2, f (x) = x3 + x2 −
6 1.8 Γ(2.8) x
+
2 0.8 Γ(1.8) x
− u(x) with homogeneous
boundary conditions. The exact solution is u(x) = x2 (1 − x). Similarly to the previous examples the GFF schemes shows major performance than GrFD (Figure 4a,b). However, it is appreciated that
the approximation order calculated with numerical schemes based on Green’s functions decreases with the fractional order of the differential equation; in other words, the approximation order for Green’s function finite differences schemes is ∼ O (hα ) (Table 1).
4.2
Dynamical simulations
Consider Eq. (3) with α = 1.5, f (x, t) = e−t
128 √ 7 π
with boundary conditions u(0, t) = 0, u(1, t) = 2et
and initial condition u(x, 0) = (x−1)x4 . For this problem the analytical solution is u(x, t) = et x4 (x−1). We used the scheme given in (21) to obtain a set of ODEs. A Runge-Kutta-Fehlberg method with ∆t = 0.01 was used to integrate the resulting set of equations. For N = 100, Figure 5a shows the profile u(x, t) considering five different times. It is appreciated that both GFF schemes show better performance than Gr¨ unwald finite differences. For instance, the error (Em ) in t = 1.0 is 0.019, 0.0044 and 0.0060 for GrFD, Caputo-GFF and RL-GFF schemes, respectively. h 24 x2.2 − Next, let us consider the case in which α = 1.8, f (x) = e−t Γ(3.2)
6 1.2 Γ(2.2) x
i + x4 − x3 with
homogeneous boundary conditions. The exact solution for this problem is u(x) = e−t x3 (x − 1). We
consider the same conditions from the previous example, Figure 5b shows the analytical and numerical solutions. It is noted that the GFF formulations yield smaller approximation errors than the GrFD scheme. In fact, for t = 1.0 the maximum error is smaller than 0.01% for the two GFF schemes whereas for the GrFD the error is four times smaller than that obtained from GFF schemes. In general, the numerical schemes based on Green’s functions formulations exhibit convergence ∼ O (hα ), whereas that GrFD exhibit a O (h). In addition, the approximation order of the proposed
schemes can be increased if we incorporate sequence acceleration methods (e.g. Richardson extrapolation), these methods have been recently applied to proposed schemes to increase the approximation 11
order [31, 33].
5
Conclusions
In this paper we proposed a new approach for the numerical solution of fractional reaction-diffusion equations. The numerical scheme is based on Green’s function formulations of the fractional operator, yielding an integral form of the equation solution. In a second step, the integral equation is approximated by means of standard quadratures (e.g., composite trapezoidal rule), yielding a set of ordinary differential equations that is numerically solved by efficient algorithms (e.g., Runge-Kutta-Fehlberg). Numerical tests with fractional equations having unique analytical solution were used to illustrate the performance of the numerical scheme, showing that the Green’s function approaches are superior in approximation error and convergence rate to finite difference schemes based on Gr¨ unwald formula.
References [1] J. Alvarez-Ramirez, F.J. Valdes-Parada, J. Alvarez, J.A. Ochoa-Tapia, A Green’s function formulation for finite-differences schemes, Chem. Eng. Sci. 62 (2009) 3083-3091. [2] Z. Bai, H. Lu, Positive solutions for boundary value problems of nonlinear fractional differential equation, J. Math. Anal. Appl. 311 (2005) 495-505. [3] B. Berkowitz, H. Scher, Anomalous transport in random fracture networks, Phys. Rev. Lett. 79 (1997) 4038-4041. [4] C.M. Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing subdiffusion, J. Comp. Physics. 227 (2007) 886-897. [5] M. Cui, Compact finite difference method for the fractional diffusion equation, J. Comp. Physics, 228 (2009) 7792-7804. [6] S. Das, Approximate solution of fractional diffusion equation-revisited, Int. Rev. Chem. Eng, 4 (2012) 501-504. [7] D. del-Castillo-Negrete, B.A. Carreras, V.E. Lynch, Non-diffusive transport in plasma turbulence: A fractional diffusion approach, Phys. Rev. Lett. 94 (2005) 065003.
12
[8] Z.Q. Deng, V.P. Singh, L. Bengtsson, Numerical solution of fractional advection-dispersion equation, J. Hydraulic Engineering. 130 (2004) 422-431. [9] H. Ding, C. Li, Y. Chen, High-order algorithms for Riesz derivative and their applications (I). Abstract Appl. Anal. 653797 (2014) 1-17. [10] H. Ding, C. Li, Y. Chen, High-order algorithms for Riesz derivative and their applications (II). J. Comput. Phys. In press, doi:10.1016/j.jcp.2014.06.007. [11] K. Friess, J.C. Jansen, J. Pozivil, V. Hanta, V. Hynek, O. Vopicka, E. Drioli, Anomalous phenomena occurring during permeation and sorption of C1 -C6 alcohol vapors in teflon AF 2400, Ind. Eng. Chem. Res. 52 (2013) 10406-10417. [12] V.V. Gafiychuk, B.Y. Datsko, Pattern formation in a fractional reaction-diffusion system, Physica A, 365 (2006) 300-306. [13] B.I. Henry, S.L. Wearne, Existence of Turing instabilities in a two-species fractional reactiondiffusion system, SIAM J. Appl. Math. 62 (2002) 870-887. [14] E. Hernandez-Martinez, J. Alvarez-Ramirez, F.J. Valdes-Parada, Formulaciones integrales para ecuaciones reacci´ on-difusi´ on generalizadas, Rev. Mex. Ing. Qui., 10 (2011) 363-373. [15] E. Hernandez-Martinez, F.J. Valdes-Parada, J. Alvarez-Ramirez, A Green’s function formulation of nonlocal finite-difference schemes for reaction-diffusion equations, J. Comput. Appl. Math., 235 (2011) 3096-3103. [16] M. Ilic, I.W. Turner, F. Liu, V. Anh, Analytical and numerical solutions of a one-dimensional fractional-in-space diffusion equation in a composite medium, Appl. Math. Comput. 216 (2010) 2248-2262. [17] H. Jafari, S. Seifi, Homotopy analysis method for solving linear and nonlinear fractional diffusionwave equation, Commun. Nonlinear Sci. 14 (2009) 2006-2012. [18] E.K. Lenzi, H.V. Ribeiro, J. Martins, M.K. Lenzi, G.G. Lenzi, S. Specchia, Non-Markovian diffusion equation and diffusion in a porous catalyst, Chem. Eng. J., 172 (2011) 1083-1087. [19] C. Li, F. Zeng, Finite differences methods for fractional differential equations. Int. J. Bifurcat. Chaos. 22 (2012) 1230014.
13
[20] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker-Planck equation, J. Comput. Appl. Math. 166 (2004) 209-219. [21] Q. Liu, F. Liu, I. Turner, V. Anh, Y.T. Gu, A RBF meshless approach for modeling a fractal mobile/immobile transport model, Appl. Math. Comput. 226 (2014) 336-347. [22] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided space-fractional partial differential equations, Appl. Numer. Math. 56 (2006) 80-90. [23] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: A fractional dynamics approach, Phys. Rep. 339 (2002) 1-77. [24] K. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential, Wiley, New York, 1993. [25] S. Momani, Z. Odibat, A novel method for nonlinear fractional partial differential equations: Combination of DTM and generalized Taylor’s formula, J. Comput. Appl. Math. 220 (2008) 85-95. [26] K.B. Oldham, J. Spanier, The Fractional Calculus, Dover Publications, New York, 2006. [27] P. Paradisi, R. Cesari, F. Mainardi, F. Tampieri, The fractional Fick’s law for non-local transport processes, Physica A 293 (2001) 130-142. [28] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [29] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives: Theory and Dpplications, Gordon and Breach Science, New York, 1993. [30] R. Schumer, M.M. Meerschaert, B. Baeumer, Fractional advection-dispersion equations for modeling transport at the Earth surface, J. Geophys. Res. 114 (2009) 1-15. [31] C. Shen, M.S. Phanikumar, An efficient space-fractional dispersion approximation for stream solute transport modeling, Adv. Water Resour. 32 (2009) 1482-1494. [32] X. Su, S. Zhang, Solution to boundary-value problems for nonlinear differential equations of fractional order, Electron. J. Diff. Equations, 26 (2009) 1-15. [33] C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comp. Physics. 213 (2006) 205-213. 14
[34] F.J. Valdes-Parada, M. Sales-Cruz, J.A. Ochoa-Tapia, J. Alvarez-Ramirez, On Green’s function methods to solve nonlinear reaction-diffusion systems, Comp. Chem. Eng. 32 (2008) 503-511. [35] F.J. Valdes-Parada, J.A. Ochoa-Tapia, J. Alvarez-Ramirez, Effective medium equations for fractional Fick’s law in porous media, Physica A 373 (2009) 339-353. [36] Y. Zhang, Z. Sunb, H. Liao, Finite difference methods for the time fractional diffusion equation on non-uniform meshes. J. Comptu. Phys. 265 (2014) 195-210. [37] Y. Zheng, C. Li, Z. Zhao, A note on the finite element method for the space-fractional advection diffusion equation, Comput. Math. Appl. 59 (2010) 1718-1726. [38] Q. Zheng, B. Yu, S. Wang, L. Luo, A diffusivity model for gas diffusion through fractal porous media, Chem. Eng. Sci., 68 (2012) 650-655.
15
6
Caption to Figures
Figure 1 Distribution of fractional Green’s function, for x = 0.5 and domain [0,1]. Green’s function obtain of a) Dirichlet-type and b) Mixed-type boundary conditions. Figure 2 Exact and numerical solutions for Eq. (30), a) profile u(x) and b) maximum error vs. number of nodes. Figure 3 Solution for Eq. (30) with homogeneous Dirichlet-type boundary conditions; a) Exact and approximate solutions and b) maximum error for integral and finite difference schemes. 6 x1.8 + Figure 4 a) Exact and numerical solutions for Eq. (30) with α = 1.2 and f (x) = x3 +x2 − Γ(2.8) 2 0.8 Γ(1.8) x
− u(x), b) maximum error vs. number of nodes .
Figure 5 Dynamic solution for Eq. (3), a) nonhomogeneous and b) homogeneous boundary conditions.
7
Caption to Tables
Table 1 Approximation order for different values of α in steady state.
Case
α
Finite differences
Caputo-GFF
Riemann-Liouville-GFF
1
1.8
0.9876
1.7154
1.7650
2
1.5
0.9950
1.4750
1.4707
3
1.2
0.9817
1.1957
1.1858
Table 1
16
Fractional Green's function, G(z x)
0.4
0.2
0.0
-0.2
-0.4
a
0.0
0.2
0.4
0.6
0.8
1.0
0.6
0.8
1.0
Fractional Green's function, G(z x)
z
0.0
-0.2
-0.4
-0.6
-0.8
b -1.0 0.0
0.2
0.4
z
Figure 1
17
1.0 Exact Solution
0.8
GrFD Caputo-GFF
u(x t)
Riemann-Liouville-GFF
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
Length,
0.8
1.0
x
1E-3
1E-4
GrFD
M
aximum error, E
m
0.01
Caputo-GFF
1E-5
Riemann-Liouville-GFF
20
40
N
umber of nodes,
N
60
80
100
Figure 2
18
0.00
u(x t)
-0.02
-0.04
-0.06
Exact Solution GrFD Caputo-GFF
-0.08
0.0
Riemann-Liouville-GFF
0.2
0.4
0.6
M
aximum error, E
m
Length,
0.8
1.0
x
0.01
1E-3
GrFD Caputo-GFF Riemann-Liouville-GFF
20
40
N
umber of nodes,
N
60
80
100
Figure 3
19
0.15
u(x t)
0.10
0.05 Exact Solution GrFD Caputo-GFF Riemann-Liouville-GFF
0.00 0.0
0.2
0.4
0.6
1.0
x
0.01
GrFD
M
aximum error, E
m
Length,
0.8
Caputo-GFF Riemann-Liouville-GFF
1E-3 20
40
N
umber of nodes,
N
60
80
100
Figure 4
20
2.0 Exact solution GrFD
u(x t)
1.5
Caputo-GFF
t=0.0
Riemann-Liouville-GFF
1.0
0.5
t=1.0
0.0 0.0
0.2
0.4
0.6
Length,
0.8
1.0
x
0.14 Exact solution
0.12
GrFD
t=0.0
Caputo-GFF
u(x t)
0.10
Riemann-Liouville-GFF
0.08
0.06
0.04 t=1.0
0.02
0.00 0.0
0.2
0.4
0.6
Length,
0.8
1.0
x
Figure 5
21
Highlights (for review)
Highlights
• Numerical solutions of fractional RD systems are given using an Green’s function formulations. • The scheme proposed exhibit global approximation orders of O(hα ). • Proposed scheme exhibit better numerical approximation than traditional schemes.
1
S0378-4754(15)00192-5 http://dx.doi.org/10.1016/j.matcom.2015.09.004 MATCOM 4236
To appear in:
Mathematics and Computers in Simulation
Received date: 30 July 2014 Revised date: 13 September 2015 Accepted date: 14 September 2015 Please cite this article as: E. Hernandez-Martinez, F. Valdes-Parada, J. Alvarez-Ramirez, H. Puebla, E. Morales-Zarate, A Green’s function approach for the numerical solution of a class of fractional reaction-diffusion equations, Math. Comput. Simulation (2015), http://dx.doi.org/10.1016/j.matcom.2015.09.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Manuscript Click here to view linked References
(Revised Manuscript No. MATCOM-D-14-00272R1) A Green’s function approach for the numerical solution of a class of fractional reaction-diffusion equations Eliseo Hernandez-Martinez1∗, Francisco Valdes-Parada2 , Jose Alvarez-Ramirez2 Hector Puebla3 and Epifanio Morales-Zarate1 1 Facultad
de Ciencias Qu´ımicas
Universidad Veracruzana, Xalapa Veracruz, M´exico 2 Departamento
de Ingenier´ıa de Procesos e Hidraulica
Universidad Aut´onoma Metropolitana-Iztapalapa, M´exico D.F., M´exico 3 Departamento
de Energ´ıa,
Universidad Aut´onoma Metropolitana-Azcapotzalco, M´exico D.F., M´exico ∗ Email:
[email protected]
September 13, 2015
Abstract Reaction-diffusion equations with spatial fractional derivatives are increasingly used in different science and engineering fields to describe spatial patterns arising from the interaction of chemical or biochemical reactions and anomalous diffusive transport mechanisms. Most numerical schemes to solve fractional reaction-diffusion equations use finite difference schemes based on the Gr¨ unwaldLetnikov formula. This work introduces a new systematic approach based on Green’s function formulations to obtain numerical schemes for fractional reaction-diffusion equations. The idea is to pose an integral formulation of the equation in terms of the underlying Green’s function of the fractional operator to subsequently use numerical quadrature to obtain a set of ordinary differential equations. To illustrate the numerical accuracy of the method, dynamic and steady-state situations are considered and compared with analytical and numerical solutions via Gr¨ unwald finite differences schemes. Numerical simulations show that the scheme proposed improves the performance ∗
Corresponding author
1
and convergence of traditional finite differences schemes based on Gr¨ unwald formula. Keywords. Anomalous diffusion; reaction-diffusion equations; Green’s function; numerical solution.
1
Introduction
Reaction-diffusion systems have been used in different science and engineering fields. However, recent research indicates that the classical diffusion equation is inadequate to model many real situations. For example, in processes where the particles spread at a faster rate than that predicted by Fickian diffusion and can present asymmetry. This phenomenon is commonly known as anomalous diffusion [23, 27]. Fractional reaction-diffusion equations can describe a large number of physical processes controlled by anomalous diffusion. Some example include, solute transport [31], transport in plasma turbulence [7], in the pattern formation [13, 12], transport at the Earth surface [30], diffusion through porous materials [3, 35]. Lenzi et al. [18] showed that the anomalous diffusion is a suitable framework for modeling transport in catalytic particles. Recently, Zheng et al. [38] considered fractal theory for modeling diffusion in porous media, finding good agreement between experimental data and theoretical predictions. On the other hand, Das [6] computed approximate analytical solutions for fractional diffusion equations that can be used in the modeling of reaction-diffusion in catalytic systems. Liu et al. [21] used a fractional differential equation to describe anomalous diffusion in a mobile/immobile transport model. Friess et al. [11] showed that non-Fickian diffusion phenomenon should be considered for describing transport of some alcohol vapors in teflon material. Different numerical methods have been used to solve both linear and nonlinear fractional reactiondiffusion equations such as, standard discretization of the fractional-in-space operator [16], Fourier methods [4], homotopy methods [17], finite element [37] and Gr¨ unwald fractional difference (GrFD). Due to its simplicity of computational implementation, the latter is the most commonly used [22, 8, 25, 5]. However, it is well-known that the discretization of fractional derivatives by means of GrFD is of first-order accurate O(h) and it may lead to instabilities [33]. For this reason, several works are focused in the derivation of stable and high-order schemes [9, 10, 36]. In [19], Li and Zeng presented a review about GrFD schemes where they discussed their application to different types of fractional differential equations. Due to the extensive applications of the fractional models, it is necessary to dispose of stable numerical methodologies with high approximation order. In this sense, recent studies have demonstrated that numerical schemes based on Green’s functions are an alternative for the solving of reaction2
diffusion equations [34]. In fact, systematic derivation of nonstandard FD schemes can be obtained by means of Green’s function formulations [1, 15]. Whereby, in this work a new systematic approach for the discretization of fractional reaction-diffusion equations based on Green’s function formulations of the fractional operator is proposed. The numerical solution of proposed scheme is based on method of lines (MOL), where a Runge-Kutta method is used for solving the ordinary differential equations (ODE) system obtained. The paper is organized as follows, first we present the integral formulation of the fractional differential equation [2, 32]. Subsequently, a quadrature rule is used to obtain a set of ordinary differential equations. In addition, depending on the discretization of the integral, a numerical scheme with equivalent structure to finite differences based on Gr¨ unwald formula can be obtained. To illustrate the numerical accuracy of the method, the results are compared with the analytical and numerical solution via GrFD.
2
Preliminaries
For convenience, we start our derivations by presenting some definitions and properties of the fractional integral and differential operators [2, 32]. Later, we describe the derivation of integral formulation for solving fractional reaction-diffusion systems. For the study of fractional differential equations, different definitions of the derivatives and integrals have been proposed [29, 28, 26]. The integral operator of fractional of order α > 0 of a function u(x) is defined as [2], Iaα u(x)
1 = Γ(α)
Z
x a
(x − z)α−1 u(x)dz
(1)
where Iaα is the so-called Riemann-Liouville fractional integral operator of order α. On the other hand, the fractional differential operator for a function u(x) is defined as [32], Daα u(x)
1 = Γ(n − α)
Z
x a
(x − z)n−α−1
dn u(z) dz dz n
(2)
where Daα is the Caputo derivative of order α and n is an integer such that n − 1 < α ≤ n. Note that
when α = 2 one recovers the classical diffusion operator. It is assumed that u(x) is a continuously
differentiable function and piecewise smooth, in order for the following considerations to be true [2, 32], • Lemma 1. For α > 1, if we assume that u(x) ∈ C(a, b) ∩ L(a, b), then the Caputo homogeneous differential equation (Daα u(x) = 0) has the following solution
u(x) = c0 + c1 x + c2 x2 + . . . + cn−1 xn−1 3
• Remark 1. For α > 1, if we assume u(x) ∈ C(a, b), then Daα Iaα u(x) = u(x) almost everywhere. • Lemma 2. If α > 1 and assuming that u(x) ∈ C(a, b) ∩ L(a, b), then Iaα Daα u(x) = u(x) + c1 + c2 x + . . . + cn−1 xn−1 for some ci ∈ R, i = 1, . . . , n, where n is the smallest integer greater than or equal to α. Anomalous diffusion coupled with chemical reactions may be described by fractional diffusion operators such as, ∂u(x, t) = Daα u(x, t) − f (u(x, t)), ∂t
x ∈ [a, b]
(3)
with initial condition as u(x, 0) = φ(x)
(4)
and the boundary conditions ∂u(a, t) + βa u(a, t) + γa = 0 ∂x ∂u(b, t) + βb u(b, t) + γb = 0 αb ∂x
αa
(5)
where f (u(x, t)) is a continuously differentiable function of its arguments and α ∈ (1, 2] is the order
of the derivative. Eq. (3) arises from a mass balance and a fractional Fick’s law where the mass flux is proportional to the fractional gradient −∇α−1 u(x, t) [35]. Note that α = 2 corresponds to the classical diffusion equation. The case α ∈ (1, 2) models a superdiffusive flow in which a cloud of diffusing particles spreads at a faster rate than the classical diffusion model. For developing the integral formulation, let us rewrite Eq. (3) as Daα u(x, t) = Ψ(x, t) where Ψ =
∂u(x,t) ∂t
x ∈ [a, b]
(6)
+ f (u(x, t)). For Lemma 2 with n = 2, if one multiplies Eq. (6) by the Riemann-
Liouville integral operator, one obtains Iaα Daα u(x, t) = u(x, t) + c0 + c1 x = Iaα Ψ(u, t) or, equivalently
1 u(x, t) = Γ(α)
Z
a
x
(x − z)α−1 Ψ(x, t)dz − c0 − c1 x
Evaluating the boundary condition given by Eqs. (5) into Eq. (8) it results that γa αa + βa βa Z b Z b αb (b − z)α−1 (b − z)α−2 αb αa 1 Ψ(z, t)dz + Ψ(z, t)dz + − b − a + αβbb − αβaa βb a Γ(α − 1) Γ(α) βb βa a
c0 = −c1 a − c1 c1 =
4
(7)
(8)
leading to u(x, t) =
a−x+
αa βa
αb βb
Z
b
(b − z)α−2 Ψ(z, t)dz + Γ(α − 1)
b − a + αβbb − αβaa a Z x 1 γa + (x − z)α−1 Ψ(z, t)dz − Γ(α) a βa
Z
b a
(b − z)α−1 αb αa Ψ(z, t)dz + − Γ(α) βb βa
(9)
After some algebraic manipulations, Eq. (9) takes the form Z b γb γa GC (z, x)Ψ(z, t)dz − σ(x) + u(x, t) = (σ(x) + 1) βa βb a where σ(x) =
x−a−
b−a+
αb βb
Γ(α)
(10)
αa βa
−
αa βa
and GC (z, x) is the Caputo fractional Green’s function, − αb (b−z)α−2 σ(x) − (b−z)α−1 σ(x) + βb Γ(α−1) Γ(α) GC (z, x) = − αb (b−z)α−2 σ(x) − (b−z)α−1 σ(x) βb Γ(α−1)
(x−z)α−1 Γ(α)
if z < x if z ≥ x
(11)
Eq. (11) is the general Green’s function, which, as a whole with Eq. (10), denote the Caputo
Green’s function formulation (Caputo-GFF) for fractional reaction-diffusion equation given by Eq. (1). In order to observe the effect of anomalous transport on the Green’s function, we consider Dirichlet-type boundary conditions (αa = αb = 0, βa = βb = 1), for x = 0.5 and z ∈ [0, 1], Figure
1a shows the fractional Green’s function distribution in the domain [0, 1], where we can observe that
for 1 < α < 2 the Green’s function is not symmetrical (anomalous diffusion) whereas for α = 2.0 the symmetrical profile is recovered (Fickian diffusion). Green’s functions tend to zero only at z = 1.0, whereas at z = 0.0 the Green’s functions are inversely proportional to α, this behavior is compensated with the contribution of the boundary conditions. For mixed-type boundary conditions (αb = βa = 0, αa = βb = 1), Figure 1b shows that the Green’s function distribution is dissipated at the Neumanntype boundary, which increase as α decreases. On the other hand, unlike traditional methods for the numerical solution of boundary value problems (e.g., finite differences), Caputo-GFF can provide an interpretation of how to carry out the transport and reaction processes. Some comments regarding the structure of Green’s function approach derived above are in order: i) Caputo Green’s function GC (x, z) can be reduced to the standard Green’s functions as α → 2 and according to the parameters αi and βi , it is possible to obtain Green’s function for different
5
types of boundary conditions (e.g., Dirichlet, Neumann, Robin). For instance, if we consider α = 2 the Green’s function is G(z, x) =
1 b − a + αβbb −
αa βa
z − a − αa x−b− βa x − a − αa z−b− βa
αb βb αb βb
if z < x if z ≥ x
(12)
obtaining the integral formulation for the classical reaction-diffusion system with general-type boundary conditions. This result was reported in [14], indicating that Eq. (11) is the general Green’s function for reaction-diffusion systems with general-type boundary conditions in a slab. ii) In the original boundary value problem the diffusive transport is described by the differential operator, while in the integral formulation the transport phenomena is described through the Green’s function. iii) The proposed scheme identifies that the profile u(x, t) consists of the interaction of the phenomena involved (accumulation, reaction and boundary conditions) and how they are distributed by Green’s function. iv) The first two terms on the right side of equation (10) describe the influence of inhomogeneous boundary conditions. v) The integral contains two contributions: the memory of the system that indicates how the concentration profiles change over the time and how the fractional Green’s function distributes the source in the domain (i.e., Green’s function describes the anomalous diffusion phenomena in time and space). When
∂u(x,t) ∂t
= 0, the profile of u(x, t) does not depend on time, then Eq.
(10) describes the profiles of Eq. (1) in steady-state. vi) The proposed scheme is a systematic methodology that can be applied to different definitions of fractional derivatives. For example, the Riemann-Liouville derivative is given by Daα u(x) = Rx 1 dn α+1−n u(z)dz and its properties described in Lemmas 1 and 2 are u(x) = Γ(n−α) dxn a (x − z) c1 xα−1 +c2 xα−2 +. . .+cn−1 xα−n−1 and Iaα Daα u(x) = u(x)+c1 xα−1 +c2 xα−2 +. . .+cN −1 xα−N −1 ,
respectively [2]. Following the procedure described in Eqs. (6)-(11), one can obtain that the Riemann-Liouville Green’s function formulation (RL-GFF) is given by Z b xα−1 γa γb η4 GRL (z, x)Ψ(z, t)dz σ(x) + − σ(x) + u(x, t) = η2 η2 βa βb a η1 xα−1 −η2 xα−2 , η1 = αβaa (α − 2)aα−3 + aα−2 , η2 η1 η4 −η2 η3 αb α−3 +bα−2 , η = αb (α−1)bα−2 +bα−1 and G 4 RL (z, x) is βb (α−2)b βb
where σ(x) =
6
=
αa βa (α
(13)
− 1)aα−2 + aα−1 , η3 =
the Riemann-Liouville fractional
Green’s function,
− αb (b−z)α−2 σ(x) − βb Γ(α−1) GRL (z, x) = − αb (b−z)α−2 σ(x) − βb Γ(α−1)
3
(b−z)α−1 Γ(α) σ(x) (b−z)α−1 Γ(α) σ(x)
+
(x−z)α−1 Γ(α)
if z < x if z ≥ x
(14)
Numerical scheme
In this section, the numerical methodologies, Gr¨ unwald finite differences scheme and Green’s function formulations, for solving fractional reaction-diffusion models given by Eqs. (3)-(5) are presented. For this, we consider ∆x = h > 0 as the spatial grid size, ∆x = (b − a)/ (N + 1), with xi = a + i∆x and
ui (t) = u(xi , t), for i = 1, 2, ..., N .
3.1
Gr¨ unwald finite differences scheme
Finite difference approximations are commonly used to numerically compute the solution of the fractional boundary value problem given by Eq. (3). So, in order to compare the proposed methodology, we used a spatial discretization of α-order fractional derivative, by mean of right-shifted Gr¨ unwald formula [24], such as
Daα ui (t) ≈ δx,a ui (t) = where the weights gα,k are given by gα,k =
i+1 1 X gα,k ui−k+1 (t) hα
(15)
k=0
Γ(k−1) Γ(−α)Γ(k+1) .
So, if the shifted Gr¨ unwald estimates are
substituted in Eq. (3), it is obtained the following i+1 1 X dui (t) gα,k ui−k+1 (t) − f (ui (t)), = α dt h
i = 1, 2, . . . , N
(16)
k=0
Note that for the classical diffusion equation, (α = 2), the weights g2,k are g2,0 = g2,2 = 1, g2,1 = −2
and for k > 2 g2,k = 0. Eq. (16) is a set of N differential equations that can be solved numerically by means of numerical integration methods to obtain the approximate solution of u(x, t).
3.2
Green’s function formulation
For this case, we rewrite the Caputo-GFF of u(x, t) given by (10) as Z b GC (z, x)Ψ(z, t)dz u(x, t) = B(x, t) +
(17)
a
where the boundary contribution B(x, t) is given by
B(x, t) = (σ(x) + 1) 7
γb γa − σ(x) βa βb
(18)
and the Green’s function G(z, x) is given by either Eq. (11). The proposed numerical scheme is obtained by forcing to zero the residual of Eq. (17) at the grid points x = xi . That is, Z b GC (z, xi )Ψ(z, t)dz, i = 1, ..., N ui (t) = B(xi , t) +
(19)
a
In the second step, the integral in Eq. (19) is approximated by a numerical quadrature. For instance, the composite trapezoidal rule gives N X ui (t) = ri (t) + h GC (zj , xi )Ψj (t), i = 1, ..., N
(20)
j=1
where ri (t) = B(xi , t) +
h [GC (a, xi )Ψa (t) + GC (b, xi )Ψb (t)])] 2
(21)
Let u = (u1 , u2 , ..., uN )T , r = (r1 , r2 , ..., rN )T and Ψ = (Ψ1 , Ψ2 , ..., ΨN )T . Then, Eq. (21) can be rewritten as u(t) = r(t) + hGΨ(t)
(22)
where the components of the N × N -dimensional matrix G are gij = GC (xj , xi ), i, j = 1, ..., N .
Substituting Ψ(t) =
du(t) dt
+ f (u(t)), we obtain the following set of N initial-value differential equations du(t) (u(t) − r(t)) = G−1 − f (x, u(t)) dt h
(23)
with initial conditions u(t0 ) = (φ(x1 ), φ(x2 ), ...., φ(xN ))T . Here, one can use a standard package for integration of initial value problems to obtain the approximate solution u(t). The MOL technique is based on the discretization of the spatial operator (using schemes discussed above), leading to a set of ODE. So, the system reduced is integrated in time using an ODE code such as, Runge-Kutta methods. The advantage of the MOL technique is that temporal accuracy can be specified by the users, therefore the error checking, robustness, order selection and time-step adaptive features available in sophisticated ODE codes can be applied to the time integration of the partial differential equations. For example, Liu et al. [20], proposed the MOL for solving of fractional Fokker-Planck equation, where backward differentiation formulas were used. Their results showed a good agreement between the numerical solution and analytic solution. Other alternative is to perform a full discretization that includes the temporal derivative. In this way, let δ = ∆t be a time step, tk = t0 + kδ and uki = ui (tk ). For example, assuming an implicit Euler discretization of the time derivative in Eq. (22); namely, du(tk ) dt
≈
uk −uk−1 . δ
This leads to the following expression: h i h k k u =r + G uk − uk−1 + δf (x, uk ) , k = 1, 2, .... δ 8
(24)
which, for each time tk , can be solved for the vector uk by means of, e.g., fixed-point or NewtonRaphson iterations. On the other hand, considering the integral formulation given by Eq. (10) and following the ideas of the Gr¨ unwald formula, one can propose a backward scheme with an equidistant spatial grid size, obtaining ui (t) =
xi − xi+1 −
xi+1 − a +
with GC (z, x) =
αb βb
αb βb
−
xi − a − αβaa γa − βa xi+1 − a + αβbb −
αa βa
a xi −a− α β
a α xi+1 −a+ βb b αa xi −a− β a α xi+1 −a+ βb b
h h
α−2
−z) − αβbb (xi+1 Γ(α−1)
−
α−2
− αβbb (x−i+1−z) Γ(α−1)
αa βa
γb + βb
(xi+1 −z)α−1 Γ(α)
−
Z
i
xi+1
+ i α−1
(xi+1 −z) Γ(α)
GC (z, x)Ψ(z, t)dz
(25)
a
(xi −z)α−1 Γ(α)
si z < x si z ≥ x
(26)
Now, one can perform the discretization of the integral in (25) by means of numerical approximations. Using the trapezoidal rule to approximate the integral and substituting Ψ(zj , t) =
duj (t) dt
+
f (uj (t)), one obtains the following nonlinear differential system N X duj (t) + f (uj (t)) , i = 1, ..., N ui (t) = ri (t) + h GC (zj , xi ) dt
(27)
j=1
where
αb βb
xi − a − αβaa γa γb ri (t) = − xi+1 − a + αβbb − αβaa βa xi+1 − a + αβbb − αβaa βb h dub (t) h dua (t) + f (ua (t)) GC (a, xi ) + + f (ub (t)) GC (b, xi ) + 2 dt 2 dt xi − xi+1 −
Let u=(u1 , . . . , uN )T , r=(r1 , . . . , rN )T and f (u(t)) = (f (u1 (t)), f (u2 (t)), ..., f (uN (t)))T . Then, Eq. (27) can be rewritten as
du(t) u(t) = r(t) + hG + f(u(t)) dt
(28)
where G is a lower triangular matrix of N xN with components di,j = GC (xj , xi ), i, j = 1, . . . , N . Eq. (28) can be presented as Eq. (23) and can be solved using numerical methods for solving of ODEs. Notice that the scheme (28) is equivalent to the scheme (16), where, in order to calculate ui , both schemes consider the values of [ua , u1 , . . . , ui , ub ] and the weights, gα,k /hα in the GrFD scheme are replaced for G−1 /h in the Green’s function.
4
Numerical examples
In this section we explore the ability of the proposed scheme to produce accurate solutions for fractional reaction-diffusion equations. To this end, we consider the integral formulation based on Caputo and 9
Riemann-Liouville definitions, given by Eqs. (10)-(11) and (13)-(14), respectively. The ODE systems was integrated by using Runge-Kutta-Fehlberg subroutine from Matlab package and the numerical results are compared with a classical GrFD schemes and the corresponding analytical solution. As an index for measuring approximation errors, the absolute maximum error Em was considered, which is defined as Em = max |ue (xi ) − uap (xi )|
(29)
where uap (xi ) is the numerical approximation and ue (xi ) is the analytical solution.
4.1
Steady state simulations
Consider the steady state fractional equation: dα u(x) = f (u(x)) dxα
x ∈ [0, 1]
(30)
To illustrate the effect of the fractional order (α) in the numerical accuracy of the proposed schemes, we study three cases: α = 1.8, α = 1.5 and α = 1.2. 4.1.1
Case 1
Consider Eq. (30) with α = 1.8, D(x) =
Γ(2.2) 2.8 6 x ,
f (x) = u(x) + x3 (x − 1) D(x)−1 and Dirichlet-
type boundary conditions, u(0) = 0 and u(1) = 1. The exact solution for this problem is u(x) = x3 . For 100 nodes, Figure 2a shows the profile u(x) obtained from Green’s function formulation using Caputo and Riemann-Liouville definitions and traditional Gr¨ unwald finite differences. The numerical solutions compare well (Em < 10−2 ) with the exact analytic solution to the fractional differential equation. Figure 2b shows the absolute maximum error (Em ) as function of the number of nodes (N ), the approximation errors are smaller for the Green’s function schemes than for Gr¨ unwald finite differences. Notice that, for GrFD the approximation order is O (h), whereas the for Caputo-GFF is O (h1.71 ) and Riemann-Liouville-GFF scheme is O (h1.76 ), indicating better numerical accuracy that GrFD. 4.1.2
Case 2
Let us fix now α = 1.5, f (x) =
128 √ x3.5 7 π
−
64 √ x2.5 5 π
with non-homogeneous Neumann-type boundary
conditions (du(0)/du = 0 and u(1) = 2) and exact solution as u(x) = x4 (x + 1). The numerical and
10
analytical profiles u(x) are shown in Figure 3a. The approximation errors are shown in Figure 3b; one can observe that the approximation order calculated with GrFD converges to O (h), whereas for Green’s function formulations it converges to ∼ O (h1.5 ). As in the previous case the minor Em is for
Caputo-GFF scheme. However, the approximation order for RL-GFF and Caputo-GFF schemes are the same. 4.1.3
Case 3
Consider Eq. (30) with α = 1.2, f (x) = x3 + x2 −
6 1.8 Γ(2.8) x
+
2 0.8 Γ(1.8) x
− u(x) with homogeneous
boundary conditions. The exact solution is u(x) = x2 (1 − x). Similarly to the previous examples the GFF schemes shows major performance than GrFD (Figure 4a,b). However, it is appreciated that
the approximation order calculated with numerical schemes based on Green’s functions decreases with the fractional order of the differential equation; in other words, the approximation order for Green’s function finite differences schemes is ∼ O (hα ) (Table 1).
4.2
Dynamical simulations
Consider Eq. (3) with α = 1.5, f (x, t) = e−t
128 √ 7 π
with boundary conditions u(0, t) = 0, u(1, t) = 2et
and initial condition u(x, 0) = (x−1)x4 . For this problem the analytical solution is u(x, t) = et x4 (x−1). We used the scheme given in (21) to obtain a set of ODEs. A Runge-Kutta-Fehlberg method with ∆t = 0.01 was used to integrate the resulting set of equations. For N = 100, Figure 5a shows the profile u(x, t) considering five different times. It is appreciated that both GFF schemes show better performance than Gr¨ unwald finite differences. For instance, the error (Em ) in t = 1.0 is 0.019, 0.0044 and 0.0060 for GrFD, Caputo-GFF and RL-GFF schemes, respectively. h 24 x2.2 − Next, let us consider the case in which α = 1.8, f (x) = e−t Γ(3.2)
6 1.2 Γ(2.2) x
i + x4 − x3 with
homogeneous boundary conditions. The exact solution for this problem is u(x) = e−t x3 (x − 1). We
consider the same conditions from the previous example, Figure 5b shows the analytical and numerical solutions. It is noted that the GFF formulations yield smaller approximation errors than the GrFD scheme. In fact, for t = 1.0 the maximum error is smaller than 0.01% for the two GFF schemes whereas for the GrFD the error is four times smaller than that obtained from GFF schemes. In general, the numerical schemes based on Green’s functions formulations exhibit convergence ∼ O (hα ), whereas that GrFD exhibit a O (h). In addition, the approximation order of the proposed
schemes can be increased if we incorporate sequence acceleration methods (e.g. Richardson extrapolation), these methods have been recently applied to proposed schemes to increase the approximation 11
order [31, 33].
5
Conclusions
In this paper we proposed a new approach for the numerical solution of fractional reaction-diffusion equations. The numerical scheme is based on Green’s function formulations of the fractional operator, yielding an integral form of the equation solution. In a second step, the integral equation is approximated by means of standard quadratures (e.g., composite trapezoidal rule), yielding a set of ordinary differential equations that is numerically solved by efficient algorithms (e.g., Runge-Kutta-Fehlberg). Numerical tests with fractional equations having unique analytical solution were used to illustrate the performance of the numerical scheme, showing that the Green’s function approaches are superior in approximation error and convergence rate to finite difference schemes based on Gr¨ unwald formula.
References [1] J. Alvarez-Ramirez, F.J. Valdes-Parada, J. Alvarez, J.A. Ochoa-Tapia, A Green’s function formulation for finite-differences schemes, Chem. Eng. Sci. 62 (2009) 3083-3091. [2] Z. Bai, H. Lu, Positive solutions for boundary value problems of nonlinear fractional differential equation, J. Math. Anal. Appl. 311 (2005) 495-505. [3] B. Berkowitz, H. Scher, Anomalous transport in random fracture networks, Phys. Rev. Lett. 79 (1997) 4038-4041. [4] C.M. Chen, F. Liu, I. Turner, V. Anh, A Fourier method for the fractional diffusion equation describing subdiffusion, J. Comp. Physics. 227 (2007) 886-897. [5] M. Cui, Compact finite difference method for the fractional diffusion equation, J. Comp. Physics, 228 (2009) 7792-7804. [6] S. Das, Approximate solution of fractional diffusion equation-revisited, Int. Rev. Chem. Eng, 4 (2012) 501-504. [7] D. del-Castillo-Negrete, B.A. Carreras, V.E. Lynch, Non-diffusive transport in plasma turbulence: A fractional diffusion approach, Phys. Rev. Lett. 94 (2005) 065003.
12
[8] Z.Q. Deng, V.P. Singh, L. Bengtsson, Numerical solution of fractional advection-dispersion equation, J. Hydraulic Engineering. 130 (2004) 422-431. [9] H. Ding, C. Li, Y. Chen, High-order algorithms for Riesz derivative and their applications (I). Abstract Appl. Anal. 653797 (2014) 1-17. [10] H. Ding, C. Li, Y. Chen, High-order algorithms for Riesz derivative and their applications (II). J. Comput. Phys. In press, doi:10.1016/j.jcp.2014.06.007. [11] K. Friess, J.C. Jansen, J. Pozivil, V. Hanta, V. Hynek, O. Vopicka, E. Drioli, Anomalous phenomena occurring during permeation and sorption of C1 -C6 alcohol vapors in teflon AF 2400, Ind. Eng. Chem. Res. 52 (2013) 10406-10417. [12] V.V. Gafiychuk, B.Y. Datsko, Pattern formation in a fractional reaction-diffusion system, Physica A, 365 (2006) 300-306. [13] B.I. Henry, S.L. Wearne, Existence of Turing instabilities in a two-species fractional reactiondiffusion system, SIAM J. Appl. Math. 62 (2002) 870-887. [14] E. Hernandez-Martinez, J. Alvarez-Ramirez, F.J. Valdes-Parada, Formulaciones integrales para ecuaciones reacci´ on-difusi´ on generalizadas, Rev. Mex. Ing. Qui., 10 (2011) 363-373. [15] E. Hernandez-Martinez, F.J. Valdes-Parada, J. Alvarez-Ramirez, A Green’s function formulation of nonlocal finite-difference schemes for reaction-diffusion equations, J. Comput. Appl. Math., 235 (2011) 3096-3103. [16] M. Ilic, I.W. Turner, F. Liu, V. Anh, Analytical and numerical solutions of a one-dimensional fractional-in-space diffusion equation in a composite medium, Appl. Math. Comput. 216 (2010) 2248-2262. [17] H. Jafari, S. Seifi, Homotopy analysis method for solving linear and nonlinear fractional diffusionwave equation, Commun. Nonlinear Sci. 14 (2009) 2006-2012. [18] E.K. Lenzi, H.V. Ribeiro, J. Martins, M.K. Lenzi, G.G. Lenzi, S. Specchia, Non-Markovian diffusion equation and diffusion in a porous catalyst, Chem. Eng. J., 172 (2011) 1083-1087. [19] C. Li, F. Zeng, Finite differences methods for fractional differential equations. Int. J. Bifurcat. Chaos. 22 (2012) 1230014.
13
[20] F. Liu, V. Anh, I. Turner, Numerical solution of the space fractional Fokker-Planck equation, J. Comput. Appl. Math. 166 (2004) 209-219. [21] Q. Liu, F. Liu, I. Turner, V. Anh, Y.T. Gu, A RBF meshless approach for modeling a fractal mobile/immobile transport model, Appl. Math. Comput. 226 (2014) 336-347. [22] M.M. Meerschaert, C. Tadjeran, Finite difference approximations for two-sided space-fractional partial differential equations, Appl. Numer. Math. 56 (2006) 80-90. [23] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: A fractional dynamics approach, Phys. Rep. 339 (2002) 1-77. [24] K. Miller, B. Ross, An Introduction to the Fractional Calculus and Fractional Differential, Wiley, New York, 1993. [25] S. Momani, Z. Odibat, A novel method for nonlinear fractional partial differential equations: Combination of DTM and generalized Taylor’s formula, J. Comput. Appl. Math. 220 (2008) 85-95. [26] K.B. Oldham, J. Spanier, The Fractional Calculus, Dover Publications, New York, 2006. [27] P. Paradisi, R. Cesari, F. Mainardi, F. Tampieri, The fractional Fick’s law for non-local transport processes, Physica A 293 (2001) 130-142. [28] I. Podlubny, Fractional Differential Equations, Academic Press, San Diego, 1999. [29] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives: Theory and Dpplications, Gordon and Breach Science, New York, 1993. [30] R. Schumer, M.M. Meerschaert, B. Baeumer, Fractional advection-dispersion equations for modeling transport at the Earth surface, J. Geophys. Res. 114 (2009) 1-15. [31] C. Shen, M.S. Phanikumar, An efficient space-fractional dispersion approximation for stream solute transport modeling, Adv. Water Resour. 32 (2009) 1482-1494. [32] X. Su, S. Zhang, Solution to boundary-value problems for nonlinear differential equations of fractional order, Electron. J. Diff. Equations, 26 (2009) 1-15. [33] C. Tadjeran, M.M. Meerschaert, H.P. Scheffler, A second-order accurate numerical approximation for the fractional diffusion equation, J. Comp. Physics. 213 (2006) 205-213. 14
[34] F.J. Valdes-Parada, M. Sales-Cruz, J.A. Ochoa-Tapia, J. Alvarez-Ramirez, On Green’s function methods to solve nonlinear reaction-diffusion systems, Comp. Chem. Eng. 32 (2008) 503-511. [35] F.J. Valdes-Parada, J.A. Ochoa-Tapia, J. Alvarez-Ramirez, Effective medium equations for fractional Fick’s law in porous media, Physica A 373 (2009) 339-353. [36] Y. Zhang, Z. Sunb, H. Liao, Finite difference methods for the time fractional diffusion equation on non-uniform meshes. J. Comptu. Phys. 265 (2014) 195-210. [37] Y. Zheng, C. Li, Z. Zhao, A note on the finite element method for the space-fractional advection diffusion equation, Comput. Math. Appl. 59 (2010) 1718-1726. [38] Q. Zheng, B. Yu, S. Wang, L. Luo, A diffusivity model for gas diffusion through fractal porous media, Chem. Eng. Sci., 68 (2012) 650-655.
15
6
Caption to Figures
Figure 1 Distribution of fractional Green’s function, for x = 0.5 and domain [0,1]. Green’s function obtain of a) Dirichlet-type and b) Mixed-type boundary conditions. Figure 2 Exact and numerical solutions for Eq. (30), a) profile u(x) and b) maximum error vs. number of nodes. Figure 3 Solution for Eq. (30) with homogeneous Dirichlet-type boundary conditions; a) Exact and approximate solutions and b) maximum error for integral and finite difference schemes. 6 x1.8 + Figure 4 a) Exact and numerical solutions for Eq. (30) with α = 1.2 and f (x) = x3 +x2 − Γ(2.8) 2 0.8 Γ(1.8) x
− u(x), b) maximum error vs. number of nodes .
Figure 5 Dynamic solution for Eq. (3), a) nonhomogeneous and b) homogeneous boundary conditions.
7
Caption to Tables
Table 1 Approximation order for different values of α in steady state.
Case
α
Finite differences
Caputo-GFF
Riemann-Liouville-GFF
1
1.8
0.9876
1.7154
1.7650
2
1.5
0.9950
1.4750
1.4707
3
1.2
0.9817
1.1957
1.1858
Table 1
16
Fractional Green's function, G(z x)
0.4
0.2
0.0
-0.2
-0.4
a
0.0
0.2
0.4
0.6
0.8
1.0
0.6
0.8
1.0
Fractional Green's function, G(z x)
z
0.0
-0.2
-0.4
-0.6
-0.8
b -1.0 0.0
0.2
0.4
z
Figure 1
17
1.0 Exact Solution
0.8
GrFD Caputo-GFF
u(x t)
Riemann-Liouville-GFF
0.6
0.4
0.2
0.0 0.0
0.2
0.4
0.6
Length,
0.8
1.0
x
1E-3
1E-4
GrFD
M
aximum error, E
m
0.01
Caputo-GFF
1E-5
Riemann-Liouville-GFF
20
40
N
umber of nodes,
N
60
80
100
Figure 2
18
0.00
u(x t)
-0.02
-0.04
-0.06
Exact Solution GrFD Caputo-GFF
-0.08
0.0
Riemann-Liouville-GFF
0.2
0.4
0.6
M
aximum error, E
m
Length,
0.8
1.0
x
0.01
1E-3
GrFD Caputo-GFF Riemann-Liouville-GFF
20
40
N
umber of nodes,
N
60
80
100
Figure 3
19
0.15
u(x t)
0.10
0.05 Exact Solution GrFD Caputo-GFF Riemann-Liouville-GFF
0.00 0.0
0.2
0.4
0.6
1.0
x
0.01
GrFD
M
aximum error, E
m
Length,
0.8
Caputo-GFF Riemann-Liouville-GFF
1E-3 20
40
N
umber of nodes,
N
60
80
100
Figure 4
20
2.0 Exact solution GrFD
u(x t)
1.5
Caputo-GFF
t=0.0
Riemann-Liouville-GFF
1.0
0.5
t=1.0
0.0 0.0
0.2
0.4
0.6
Length,
0.8
1.0
x
0.14 Exact solution
0.12
GrFD
t=0.0
Caputo-GFF
u(x t)
0.10
Riemann-Liouville-GFF
0.08
0.06
0.04 t=1.0
0.02
0.00 0.0
0.2
0.4
0.6
Length,
0.8
1.0
x
Figure 5
21
Highlights (for review)
Highlights
• Numerical solutions of fractional RD systems are given using an Green’s function formulations. • The scheme proposed exhibit global approximation orders of O(hα ). • Proposed scheme exhibit better numerical approximation than traditional schemes.
1