Computational solutions of distributed oder reaction-diffusion systems associated with Riemann-Liouville derivatives

This article is in continuation of our earlier article [37] in which computational solution of an unified reaction-diffusion equation of distributed order associated with Caputo derivatives as the time-derivative and Riesz-Feller derivative as space derivative is derived. In this article, we present computational solutions of distributed order fractional reaction-diffusion equations associated with Riemann-Liouville derivatives of fractional orders as the time-derivatives and Riesz-Feller fractional derivatives as the space derivatives. The method followed in deriving the solution is that of joint Laplace and Fourier transforms. The solution is derived in a closed and computational form in terms of the familiar Mittag-Leffler function. It provides an elegant extension of the results given earlier by Chen et al. [1], Debnath [3], Saxena et al. [36], Haubold et al. [15] and Pagnini and Mainardi [30]. The results obtained are presented in the form of two theorems. Some interesting results associated with fractional Riesz derivatives are also derived as special cases of our findings. It will be seen that in case of distributed order fractional reaction-diffusion, the solution comes in a compact and closed form in terms of a generalization of the Kamp\'e de F\'eriet hypergeometric series in two variables, defined by Srivastava and Daoust [46] (also see Appendix B). The convergence of the double series occurring in the solution is also given.

In a recent article, Chen et al. [33] have derived the fundamental and numerical solution of a reaction-diffusion equation associated with the Riesz fractional derivative as the space derivative. Reaction-diffusion models associated with Riemann-Liouville fractional derivative as the time derivative and Riesz-Feller derivative as the space derivative are recently discussed by Haubold et al. [21]. Such equations in case of Caputo fractional derivative are solved by Saxena et al. [27]. The main object of this article is to investigate the computational solutions of fraction reaction-diffusion Equations (1) and (24) below. The results are obtained in a closed and computational forms. Due to the general character of the derived results, many known results given earlier by Chen et al. [33], Haubold et al. [22] and Pagnini and Mainardi [34], Saxena et al. [27], readily follow as special cases of our derived results.

Solution of Fractional Reaction-Diffusion Equations
In this section, we will investigate a computational solution of the one-dimensional fractional reaction-diffusion Equation (1) given below, containing Riemann-Liouville derivative as the time-derivative and a finite number of Riesz-Feller derivative as the space derivatives. The results obtained are in a compact and computational form in terms of the generalized Mittag-Leffler function, defined by Equation (5) below in the form of the following theorem: Theorem 1. Consider the following one-dimensional non-homogeneous unified fractional reaction-diffusion model associated with time-derivative 0 D α t defined by Equation (A1), n ∈ N , and Riesz-Feller space-derivatives x D γ 1 θ 1 , ..., x D γn θn , defined by Equation (A3): where t > 0, x ∈ R; α, θ 1 , ..., θ n , γ 1 , ..., γ n are real parameters with the condition: with initial conditions Here 0 D ν t N (x, t) means the Riemann-Liouville fractional partial derivative of N (x, t) with respect to t of order ν evaluated at t = 0, ν = α − 1, α − 2; x D γ 1 θ 1 , ..., x D γn θn are the Riesz-Feller space fractional derivatives with asymmetries θ 1 , ..., θ n respectively. Further, 0 D α t is the Riemann-Liouville time-fractional derivative of order α; µ 1 , ..., µ n are arbitrary constants; f (x), g(x) and φ(x, t) are given functions. Then for the solution of Equation (1), subject to the above conditions, there holds the formula where E α,β (z) is the generalized Mittag-Leffler function, defined by Wiman (see Erdélyi et al. [35], Dzherbashyan [36]) in the following form: Proof. In order to derive the solution of Equation (1), we introduce the joint Laplace-Fourier transform in the formÑ where (s) > 0, k > 0. If we apply the Laplace transform with respect to the time variable t, Fourier transform with respect to space variable x, use the initial conditions given in Equations (2) where, according to the conventions followed, the symbols ∼ will stand for the Laplace transform with respect to time variable t and * represents the Fourier transform with respect to space variable x. Solving forÑ * (k, s) it yieldsÑ * (k, s) = f * (k) where b = n j=1 µ j ψ θ j γ j (k). On taking the inverse Laplace transform of Equation (7) by means of the formula The required solution of Equation (4) is obtained by taking the inverse Fourier transform of Equation (9). This completes the proof of Theorem 1.

Special Cases
(i) If we take θ 1 = ... = θ n = 0 then by virtue of the relation (A11), we obtain the following corollary.
where t > 0, x ∈ R; α, γ 1 , ..., γ n are real parameters with the constraints with initial conditions where the various terms are as defined in Equation (3). Then there holds the formula where E α,β (z) is the generalized Mittag-Leffler function, defined by Equation (5) and c = n j=1 µ j |k| γ j .
Note that when g(x) = 0 then Theorem 1 yields Equation (9) where the middle term involving g * (k) will be absent.

Corollary 2. Consider the same equation in Equation (10) under the conditions
with the initial conditions where δ(x) is Dirac delta function. Then for the fundamental solution there holds the formula where E α,β (z) is the generalized Mittag-Leffler function defined by Equation (5) and b is given in Equation (6).
For θ = 0, Equation (20) reduces to the result given by Saxena et al. [8] in a slightly different form. On the other hand, if we further set α = 1 we obtain the following result given by Chen et al. [33] in a different form: with the initial conditions where µ is a diffusion constant, µ, t > 0, β are real parameters and δ(x) is Dirac delta function. Then for the fundamental solution of Equation (21) with initial conditions Equation (22) there holds the formula Remark 1. The result obtained by Chen et al. [33] is in terms of the Fourier transform of the Mittag-Leffler function, whereas our Equation (23) is in terms of the H-function in a closed and computable form. It is interesting to observe that for n = 1, Theorem 1 reduces to one given by Haubold et al. [22].

Further Results on Distributed Order Reaction-Diffusion Systems
In this section we will investigate a computational solution of a distributed order reaction-diffusion equation containing two Riemann-Liouville derivatives. The solution is obtained in terms of generalized Mittag-Leffler function due to Prabhakar [37]. The main result can be expressed in terms of Srivastava-Daoust hypergeometric function of two variables [38].
Theorem 2. Consider the following unified one-dimensional non-homogeneous reaction-diffusion equation of fractional order: where t > 0, a, x ∈ R; α, β, θ 1 , ..., θ n are real parameters with the constraints with the initial conditions where the various quantities are as defined in Theorem 1, and f 1 (x), f 2 (x), g 1 (x), g 2 (x) and φ(x, t) are given functions. Then for the solution of Equation (24), subject to the above conditions in Equations (25) and (26), there holds the formula where E γ α,β (z) is the generalized Mittag-Leffler function ( [37]) Proof. If we apply the Laplace transform with respect to the time variable t, Fourier transform with respect to space variable x and use the initial conditions in Equations (25) and (26) and the Equation (28), then the given equation transforms into the form Solving forÑ * (k, s) it yields where b is defined in Equation (6). To invert the Equation (30), we first invert the Laplace transform and then the Fourier transform. Thus to invert the Laplace transform we use the formula given in [2]: where (s) > 0, (α) > 0, (β) > 0, (α − ρ) > −1, (α − β) > 0, | as β b+s α | < 1; and the convolution theorem of the Laplace transform to obtain Now, the application of the inverse Fourier transform gives the following solution: Alternative form of the solution of Equation (27) By virtue of the series representation of the generalized Mittag-Leffler function E α β,γ (z) defined in Equation (28), the expression can be written as where S(·) is the Srivastava-Daoust generalization of the Kampé de Fériet hypergeometric series in two variables [14], for its definition, see (B2) in Appendix B. Hence, Theorem 2 can now be stated in terms of the Srivastava-Daoust hypergeometric function of two variables in the following form: Under the conditions of Theorem 2, the unified one-dimensional fractional reaction-diffusion equation has the solution given by × S 1:0;0 1:0;0 −at α−β , −bt α where (α) > 0, (β) > 0, (α − β) > 0 and b is defined in Equation (6).
Note 4.1. By virtue of the lemma given in the Appendix B, the double infinite power series occurring in Theorem 2 converges for (α) > 0, (α − β) > 0.

Special Cases of Theorem 2
If we set g 1 (x) = g 2 (x) = 0 in Theorem 2 then it reduces to the following: Corollary 5. Consider the following unified one-dimensional non-homogeneous reaction-diffusion equation of fractional order: where t > 0, a, x ∈ R; α, β, θ 1 , ..., θ n , γ 1 , ..., γ n are real parameters with the constraint with the initial conditions where the various quantities are as defined in Theorem 2. Then we have a special case of Equations (33) and (35) with the corresponding changes.
where [α] means the integral part of the number α.
The Laplace transform of the Riemann-Liouville fractional derivative is given by the following: (see Oldham and Spanier ([39], Kilbas et al. [11], Podlubny [40], Samko et al. [41], Mainardi [12], Diethelm [10]) This derivative is useful in deriving the solutions of integral equations of fractional order governing certain physical problems of anomalous reaction and anomalous diffusion. In this connection, one can refer to the monograph by Dzherbashyan [36], Podlubny [40], Samko et al. [41], Oldham and Spanier [39], Kilbas et al. [11], Mainardi [12], Diethelm [10] and a recent paper on the subject [9]. Following Feller [42,43], it is convenient to define the Riesz-Feller space-fractional derivative of order α and skewness θ in terms of its Fourier transform as Further, when θ = 0, we have a symmetric operator with respect to x that can be interpreted as This can be formally deduced by writing −(k) α = −(k 2 ) α/2 . For 0 < α < 2 and |θ| ≤ min{α, 2 − α}, the Riesz-Feller derivative can be shown to possess the following integral representation in x domain: For θ = 0 the Riesz-Feller fractional derivative becomes the Riesz fractional derivative of order α for 1 < α ≤ 2 defined by analytic continuation in the whole range 0 < α ≤ 2, α = 1 (see Gorenflo and Mainardi [44]) as where λ = 1 2 cos(απ/2) ; The Weyl fractional integral operators are defined in the monograph by Samko et al. [41] as Note A1. We note that x D α 0 is a pseudo differential operator. In particular, we have The H-function is defined by means of a Mellin-Barnes type integral in the following manner [13]: and an empty product is interpreted as unity; m, n, p, q and C being the complex number field. A comprehensive account of the H-function is available from Mathai et al. [13] and Kilbas et al. [11]. We also need the following result in the analysis that follows: Haubold et al. [21] have shown that { A j=1 Γ(a j + mθ j + nφ j )}{ B j=1 Γ(b j + mψ j )}{ B j=1 Γ(b j + nψ j )}x m y n { C j=1 Γ(c j + mδ j + n j )}{ D j=1 Γ(d j + mη j )}{ D j=1 Γ(d j + nη j )}m!n!
where the coefficients θ 1 , ..., θ A , ..., η 1 , ..., η D > 0. For the sake of brevity (a) is taken to denote the sequence of A parameters a 1 , ..., a A with a similar interpretation for (b), ..., (d ). Srivastava and Daoust [38] have shown that the series (B2) converges for all x, y ∈ C when and For a detailed account of the convergence conditions of the double series, see Srivastava and Daoust [38].

Conflicts of Interest
The authors declare no conflict of interest.