Diffusion with space memory modelled with distributed order space fractional differential equations

Distributed order fractional differential equations (Caputo, 1995, 2001; Bagley and Torvik, 2000a,b) were fi rst used in the time domain; they are here considered in the space domain and introduced in the constitutive equation of diffusion. The solution of the classic problems are obtained, with closed form formulae. In general, the Green functions act as low pass fi lters in the frequency domain. The major difference with the case when a single space fractional derivative is present in the constitutive equations of diffusion (Caputo and Plastino, 2002) is that the solutions found here are potentially more fl exible to represent more complex media (Caputo, 2001a). The difference between the space memory medium and that with the time memory is that the former is more fl exible to represent local phenomena while the latter is more fl exible to represent variations in space. Concerning the boundary value problem, the difference with the solution of the classic diffusion medium, in the case when a constant boundary pressure is assigned and in the medium the pressure is initially nil, is that one also needs to assign the fi rst order space derivative at the boundary. Mailing address: Prof. Michele Caputo, Dipartimento di Fisica, Università di Roma «La Sapienza», P.le A. Moro 2, 00185 Roma, Italy; e-mail: mcaputo@g24ux.phys. uniroma1.it

Fractional differential equations consisting of sums of fractional order derivatives have been extensively studied by Podlubny (1994aPodlubny ( ,b, 1999)).Fractional differential equations where the order of differentiation has been integrated over a given range, and therefore there is no single order of differentiation, have been used and studied (Caputo, 1967(Caputo, , 1995) ) to represent phenomena where the order of differentiation varies in a given range.Bagley and Torvik (2000a,b) studied extensively these Distributed Order Differential Equations (DODE), developed a very interesting mathematical theory and gave series expansion solutions.
In ordinary differential equations, where the derivatives appear with integer order, the order of the equation is given by the maximum order of differentiation appearing in the equation.The same could be true also when the orders of differentiations are real; however when this order is smaller than unity it makes no difference on the properties of the solutions.The same applies when the order of differentiation is distributed in a range 0 < a < b < 1 and therefore the name DODE emphasises the property that the order of differentiation is distributed in a continuum.Caputo (2001a) solved the classic problems of anelastic, dielectric media and of diffusion when the classic constitutive equations of these media are described in the time domain with DODE.
The scope of this note is to fi nd more fl exible models of the diffusion of fl uid introducing space domain DODE in the constitutive equa-tion of diffusion which should allow more fl ex-ible representations of non local phenomena.Specifi cally we will use DODE in modelling solutions of problems in diffusion which have been already considered in the literature with constitutive equations containing time domain fractional derivatives of a given order (Wyss, 1986;Mainardi, 1996;Caputo, 2000) and also with space domain Fractional Order Differential Equations (FODE) (Caputo and Plastino, 2002;Mainardi et al., 2001), assuming now that the fractional order of space differentiaton is integrated over a given range generating thus a DODE.
The solutions to the problems of diffusion are here found in a half space where the boundary values are applied to the plane limiting it; among others is the computation of the Green function of the pressure in the case when the pressure in the half space is initially constant with the fi rst and second order space derivative.
The problem is solved with the use of the Laplace Transform (LT) and closed form formulae are obtained in the space-time domain.
It has been shown by Cisotti (1911) that the memory introduced in the dielectric constant implies that the phenomenology of the medium is irreversible.It has also been shown also that in anelastic media the stress strain relations with FODE are causal (Bagley, 1991).
In an extensive study of energy storage in electric networks, Jacquelin (1991) used the complex frequency dependent impedance represented by transforming into a FODE the relation between the induced electric current and the electric fi eld.Jacquelin's (1991) discussion is practically extended to almost all possible circuits; he illustrates his results in the frequency domain at steady state, when the Fresnell terms of the fractional order derivatives are nil, and his technique, based on the observations at many frequencies, relies on convolution for the retrieval of the response in the time domain.
The studies of diffusion in porous media are very important in various fi elds of science and also for social needs especially in the case of water and pollutants.The diffusion is currently modelled to represent the diversifi ed forms of deviations from the classic constitutive law and several complex mathematical methods for modelling the phenomena have been introduced, ranging from non linearity to statistical mechanics and memory formalisms.Christakos et al. (1995) used stochastic diagrammatic analysis of groundwater flow in heterogeneous porous media, Hu and Cushman (1991) used a non equilibrium statistical mechanical derivation of a non-local Darcy's law for unsaturated/saturated flow, Kabala and Sposito (1991) studied the diffusion using a stochastic model of reactive solute transport with time-varying velocity in heterogeneous aquifers.
Recently, Cushman and Moroni (2001) presented a theory with statistical mechanics origin to simulate Fickian, quasi-Fickian, convolution Fickian and fractional Fickian fl uxes fi nding a dispersive fl ux which is a convolution in space.
We may also add that appropriate experimental data are not yet available to verify most theories.
In this note, we will consider an extension of the constitutive relation of diffusion to the case when a space memory mechanism operating in the medium is represented by FODE whose order covers a continuum in a given range; in which case the constitutive equation becomes a DODE.
The time fractional order derivative of the pressure represents the local variations and is particularly valid when considering local phenomena.In an infi nite medium it more appropriate to introduce the space fractional order derivative to represent the effect of the medium and its space interaction with the fl uid.Therefore, the fl ow is not directly related to the instantaneous pressure gradient in the measurement site but to the spatial fractional derivative i.e. to the pressure gradient investigated by the fl uid in the path from the starting point to the measurement site.
We will consider the porous medium fi lling a half space where the initial conditions are assigned and the boundary conditions are assigned in the plane bounding it and fi nd the closed form solution of the initial and boundary value problems.

The constitutive equation with space derivative model and the solution of the governing equation
In order to fi nd the general solution of the problem, that is the pressure distribution in the porous media, we begin setting the constitutive equations.The fi rst equation is the continuity equation between the time variation of the density and the divergence of the fl ux (2.1) Another constitutive equation is that relating the pressure to the variation of the density from its undisturbed condition (2.2) Successively, to take into account the observed deviations of the fl ow from those implied by the classic diffusion equation we introduce, as follows, a space memory formalism in the classic Darcy's law: , where u is the LT variable, we obtain from eq. ( 2.4) the following equation: (2.6) where .Proceeding to the solution we take the LT of eq.(2.5) with respect to the t variable and obtain where n is the new LT variable of the LT with respect to t.From (2.6') follows: (2.7) , 0 0 = LT .Substituting n = iw in (2.7) we see that in the diffusion the pressure is affected by a low pass fi lter.
The solution p(x,t) is then obtained by inverting both LT.The LT t v , −1 of eq.(2.7) gives 0 £ n < 1, where clearly the memory mechanism is applied to the pressure and not to the pressure gradient.The defi nition of derivative of fractional order n is (Caputo, 1969) where w is the integration variable.The factor A(n), as in the work of Caputo (1995), is a weight of the fractional derivative of order n in the range [a,b].The factor 1/(b-a) normalizes the integration with respect to n and implies that when b AE a the integral in eq. ( 2.3) gives ∂ ∂ , with 0 < a < 0.5, which is the case studied by Caputo and Plastino (2002).
We note that in eq. ( 2.3) the fractional derivative (the memory mechanism) operates on the pressure rather than on the pressure gradient (the fl ux), which was done in previous work (Caputo and Plastino, 2002); in fact it seems easier to conceive the space memory limiting its mechanism to a rate of variation of lower order (order n (0 < n < 1)) of the pressure rather than to a rate of higher order (order 1 + n) as when applying it to the gradient of the pressure.
Replacing p(x, t)/k from eq. (2.2) in eq.(2.1) and taking into account the derivative with respect to x variable in eq. ( 2.3) we obtain a single equation in p(x,t) . (2.4) In order to normalize the dimensions, we will here assume A(n) = s n − , where s has the dimension of a length.Since 0 < a < n < b < 1/2, choosing s near to unity implies that all fractional derivatives of order in the range [a,b] have practically the same weight factor.
In order to solve eq. ( 2.4) we fi rst take its Laplace Transform (LT) respect to x variable; using the LT theorem (Caputo, 1969) (2.5) where the asterisk indicates convolution and its subscript indicates that the convolution is made with respect to the variable t; performing the As shown in the Appendix the LT t v , −1 of eq.(2.8') gives fi nally (2.9) where the asterisks indicates convolution and the subscripts t or x indicate that the convolution is made with respect to the variables t or x respectively.
Separating initial conditions from boundary conditions we fi nd (2.9') (2.9') is the formal LT domain solution of the problem and includes the boundary condition p(0,t), in the terms of the 3nd, and 8th lines and the boundary condition p(0,t)/ x in the 5th, 6th and 10th lines the initial condition p(x,0) appears in the terms of 1st and 2nd line.The inverse LTs appearing in (2.9') are computed in the Appendix and give the fi nal solution in the time and space domain expressed in terms of Bestimmte integrals.These inverse LT have vague similarities with the Mittag Leffl er function and stronger similarities with the functions found by Caputo (2001a) and appearing in the solution of the diffusion problem when the memory of the corresponding constitutive equation is expressed in the form of a DODE in the time domain.
It is to be noted that the boundary conditions depend on the pressure on the boundary p(0,t) but also on its fi rst space derivatives p(0,t)/ x which may be arbitrarily assumed.The presence of a p(0,t)/ x in (2.9') is due to the use of the LT in the x domain which is applied to a derivative of  order 1 + n and therefore, according to theorem (2.5), implies in (2.9') the presence of the 1st order derivatives of p(x,0) with respect to x.The presence of b p(0,t)/ x in (2.9') is also due to the use of the LT in the x domain which is applied to a derivative of order 2 and therefore, according to theorem (2.5), implies in (2.9') the presence of the 1st derivative of p(x,0) with respect to x.
We also note that in eq.(2.9') we have two types of convolution, one relative to the time variable and one relative to the space variable.
Moreover, we note that the fi rst term in eq.(2.9') takes into account the initial values of the pressure in the medium while the other terms take into account the boundary values.The computation of the term with the initial value implies the convolution relative to the space variable only while the computation of the terms relative to the boundary values imply convolutions relative to the time variable only.
As a check it is seen that eq.(2.9'), when b AE a, gives the particular case when the constitutive eq. ( 2.3) is a space FODE which was studied by Caputo and Plastino (2002); in order to see this in detail it is suffi cient to substitute b = a + e, compute the limit of (2.9') for e AE 0 and fi nd that u a in (2.9') substitutes the With the extreme value theorem applied to eq. (2.8') it is verifi ed that p(0,0) = 0.

The initial value problem
A case of interest is when the pressure and its fi rst respect to the x variable on the boundary (x = 0) are nil for any t while the pressure in the medium p(x,t) has initial (t = 0) constant value p(x,0).We will fi rst assume a 0 and b = 0 and solve later the general case a 0 and b 0.
The solution is then readily obtained from eq. (2.9') considering only the fi rst term as in the work of Caputo and Plastino (2002).
Assuming that the pressure in the medium has initial (t = 0) constant value C 0, performing the convolution and grouping the exponential terms in (3.2) we obtain With the extreme value theorem applied to eq. (2.7) it is verifi ed that p(x,0) = C. ) The solution of the problem with b 0, a 0 is found considering that the portion of eq.(2.9') of interest in this section is (3.2'') which may be written and seeing also that in the Appendix is found and fi nally that LT x u tku , exp( ) β is known since it is the solution of the classic initial value case without memory.
The effect of the memory is obviously represented by the exponential and sine terms which contain a in eq.(3.2').
An analysis of the integrand of eq.(3.2') shows again that in order to ensure the convergence, in the exponential term containing the time t, it must be 0 < b < 1/2 which implies a negative decreasing value of the exponent when r is relatively large.
In (3.2') we see that the increase of p(x,t), with increasing, x is affected by a only indirectly, that is, at a given time p(x,t) increases with x at a rate dependent of a only through the sine and the exponential terms, in the integrand of (3.2'), which are not depending on x.

The boundary value problem
In order to give a better description of the solution (2.9'), we will specialize it to the case when p(0,t) is assigned and p(x,0) = p(0,t)/ x= = 0. Note that p(0,t)/ x = 0 does not imply that the fl ux is nil at the boundary since the fl ux is given by eq. ( 2.3).Equation (2.9') is then (4.1) In order to compute the LT x u , −1 in eq. ( 4.1) we assume that Y(u,t) is the analyitic function of the complex variable u resulting from the product of the two functions in braces of eq. ( 4.1) and consider eq.(A.2) and (A.3) of the Appendix  The discussion of the case when b 0, a 0 made in the Section 4 may be applied also to the case studied in this section.The solution of the boundary value problem with p(0,t) = C is readily obtained also from eq. (3.2').In fact indicating with p 3 (x,t) the solution (3.2'), the solution of the boundary value problem is C -p 3 (x,t).

Conclusions
The major difference between the diffusion in medium with a constitutive equation of the type (2.3) and that whose constitutive equation contains a fractional derivative of given order (Caputo and Plastino, 2002), is that the solutions found here are potentially more fl exible to represent more complex systems (Caputo, 2001a).
The difference between the space memory medium and that with the time memory is that this is fl exible to represent local phenomena while the latter is more fl exible to represent variations in space.Concerning the initial value problem the solution (3.2') may be compared with that of Caputo and Plastino (2002) who used the constitutive eq. ( 2.4) with a single fractional order derivative of order a instead of a DODE.The solution of Caputo and Plastino (2002) is (5.1) which is valid for 0 < a < 1/2.
The solution (5.1) is similar to (3.2'), it is the presence of the additional free parameter b which possibly makes (3.2') more fl exible than (5.1) to represent more complex types of diffusion.
With more detail the two solutions differ for the presence in (3.2') of the terms, (rs) 1+b cosbp -(rs) 1+a cosap, (rs) 1+b sinbp -(rs) 1+a sinap and of their weighted differences, which represents the average effect of the memory when the order of differentiation is integrated over the range [a,b].
The discussion made for the initial value problem solution, and its comparison with that of the case when in the constitutive eq.(2.4) a FODE is used instead of a DODE, is valid also for the boundary value problem solution represented by eqs.(4.1) through (4.5).
We conclude that besides the difference in the number of parameters which may be used to represent the more complex types of diffusion, there is no great difference between the solution of the DODE constitutive eq.(2.4) used here and the solution found by Caputo and Plastino (2002) in the case when a FODE constitutive eq.(2.4) is used.
The introduction of fractional order derivatives in the equations governing economy and fi nance have been proved successful also modeling financial phenomena (Caputo and Kolari, 2001;Caputo, 2001b).Appendix.
In order to compute the inverse LT of eq.(2.9') we set u = r exp (iq) (A.1) in the complex plane of fi g. 1 and integrate along the path shown there.We then note that the contribution to the integral from the inner circle of radius r 1 when r 1 0 → , is nil since the path of integration of the fi nite function is nil; also when the radius of the outer circle R, when R → ∞ the contribution of the integral is nil since the function itself is nil.The integral is then reduced to the integration along the two lines parallel to the negative real axis and to the line in the real plane and parallel to the imaginary axis which give the LT u x , −1 of eq.(2.9').In the points of the line from -• to 0 (q = p) the integrand is the complex conjugate of the value in the same points on the line from 0 to • (q = -p).
In order to compute the LT x u , −1 in eq. ( 4.1) we assume that Y(u,t) is the analytic function of the complex variable u resulting from the product of the two functions in braces of eq.(4.1).2.3)).ak (s -1 m 2+n ) pseudodiffusivity.r (x,t) (kg m -3 ) fl uid density.w, e imaginary and real parts in the plane of the integral in eq.(A.1) and shown in fi g.A.1.

Glossary
The difference of the two, in all terms appearing in eq.(2.9'), may be reduced to the form We now proceed to the computation of the inverse LT appearing in (2.9') beginning with the fi rst term (the factor convolving p(x,0) in lines 1 and 2) of eq.(2.9') (A.4) we fi nd l = 1, m = 0, (A.5) For the second term (the factor convolving ap(0,t) in lines 3 and 8) of eq.(2.9') we fi nd c, d, l, m (A.7) For the third term (the factor convolving a p(0,t)/ x in lines 5 and 10) of eq.(2.9') we fi nd again c and d as in eq.(A.
b AE a, is u a .
of (A.2) of Appendix for the integrand and substituting the values of c, d, l, m of (A.5) and (A.6) in the Appendix in eq.(A.3) of the same Appendix; we fi nd(3.2)Oncemore one may verify the computations assuming b = a + e, computing the limit of (3.2) when e AE 0 and fi nding that the exponent of the exponential becomes

(
eqs. (A.7) and (A.8) of the Appendix, the values of l, m, c, d t) = d(t), we obtain from eqs. (4.2) the Green function of the problem.