223_234 Caputo.pdf 223 ANNALS OF GEOPHYSICS, VOL. 46, N. 2, April 2003 Diffusion with space memory modelled with distributed order space fractional differential equations Michele Caputo Dipartimento di Fisica, Università di Roma «La Sapienza», Roma, Italy Abstract 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 Key words distributed order – fractional order – differential equations – constitutive equations – diffusion – space fractional derivative 1. Introduction Fractional order derivatives to model phys- ical phenomena have been widely used in the last few decades in basic mathematics (Graffi , 1936; Oldham and Spanier, 1974; Ross, 1975; Caputo, 1994, 1995; Kiryakova, 1994; Podlubny, 1994a,b, 1999; Flores and Osler, 1999; Luchko and Gorenfl o, 1999; Bagley and Torvik, 2000a,b) and in applied research (Le Mehaute and Cré- py, 1983; Jacquelin, 1991; Wyss, 1986; Caputo, 1989, 2000; Körnig and Müller, 1989; Mainardi, 1996; Caputo and Plastino, 2002; Belfi ore and Caputo, 2000). Fractional differential equations consisting of sums of fractional order derivatives have been extensively studied by Podlubny (1994a,b, 1999). Fractional differential equations where the order of differentiation has been integrat- ed over a given range, and therefore there is no single order of differentiation, have been used and studied (Caputo, 1967, 1995) to represent phenomena where the order of differentia- tion 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 224 Michele Caputo 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 inte- grated 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 frequen- cy 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. In general, we may state that the formalism of FODE is now spread to many linear studies of dissipative and dispersive properties of many anelastic and dielectric media (e.g., Caputo, 1967, 1989, 1996; Le Mehaute and Crépy, 1983; Körnig and Müller, 1989; Jacquelin, 1991; Mainardi, 1996). 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 dia- grammatic analysis of groundwater flow in heterogeneous porous media, Hu and Cushman (1991) used a non equilibrium statistical me- chanical 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) pre- sented 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 experi- mental 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 225 Diffusion with space memory modelled with distributed order space fractional differential equations 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 phe- nomena. 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. 2. 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: Q u t uq t P u t k( , ) ( , ) ( , ) / .− + =0 0 ∂ ∂ ∂ρ ∂q x t x x t t( , ) / ( , ) /+ = 0 p x t k x t( , ) ( , ).= ρ (2.3) q b a A n p x dn p xn n a b = − − [ ] −∫( / ( )) ( ) / /α ∂ ∂ β∂ ∂ Q u t b a us P u t us p t us us us uP u t p t n b a ( , ) ( / ( ))( ) ( , ) ( ) ( , ) (( ) ( ) ) / log ( , ) ( , ) = − − −[ ] −[ ] − −[ ] −α β 1 0 0 uq t P u t k t b a us P u t us p t us us us uP u t p t n b a ( , ) ( , ) / ( / ( ))( ) ( , ) ( ) ( , ) (( ) ( ) ) / log ( , ) ( , ) 0 0 0 1 − = = − − −[ ] −[ ] − −[ ] − ∂ ∂ α β ukq t kp t k b a us p t us us us P u t t P u t k b a us us us us u k b a n b a ( , ) ( , ) ( / ( ))( ) ( , ) (( ) ( ) ) / log ( , ) / ( , ) ( / ( ))( ) (( ) ( ) ) / log 0 0 01− − − ] −[ ] = + − − − −[ ] −β α ∂ ∂ α β LT ukq t kp t k b a us p t us us us k b a us us us us u k V u v P u t v b a n b a , , , / , / log / / / log , , / 0 0 0 0 1 ( ) − ( ) − −( )( )[{ ( ) ( )] ( ) − ( )( )[ ]} − −( )( )( ) ( ) − ( )( )[{ + ]} = ( ) − ( ) − − β α ν α β ν ααk /([{ βb a us us us us u kn b a / log−( )( ) ( ) − ( )( ) + )]} Q u t b a us P u tn a b ( , ) ( / ( )) ( ) ( , )= − − +[[∫α us p t dn p x( ) ( , ) /− ] −− β∂ ∂1 0 ukq t kp t k b a us p t us us us k b a us us us us u k P u t P u k b a us b a n b a n ( , ) ( , ) ( / ( ))( ) ( , ) (( ) ( ) ) / log exp / ( ))( ) (( ) ( ) ) / log ( , ) ( , ) exp ( / ( ))( ) 0 0 0 0 1− − −[ ]{ −[ ]} [{ − − + ]} = = − − −β α α β α (((( ) / log us us us u k b a +[{ −( ) ) + ]}β P u t P u k b a us us us us u k ukq t kp t k b a us p t us us us k b a us n b a b a n ( , ) ( , ) exp ( / ( ))( ) (( ) ( ) ) / log ( , ) ( , ) ( / ( ))( ) ( , ) (( ) ( ) ) / log exp ( / ( ))( ) = − +[{ − + }] + − +[{ − − ] − ][ } − − 0 0 0 01 α β β α α (((( ) ( ) ) / log us us us u k b a−[{ + ]}β Q u t uq t P u t k t( , ) ( , ) ( , ) /= − =0 0∂ ∂ 226 Michele Caputo with m = 1, 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 (2.6') where n is the new LT variable of the LT with respect to t. From (2.6') follows: (2.7) where V u v P u tt v( , ) ( , ),= LT and P u p xx u( , ) ( , ),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 LTt 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 Æ a the integral in eq. (2.3) gives ∂ ∂a np x t x( , ) / , with 0 < a < 0.5, which is the case studied by Caputo and Plastino (2002). We note that in eq. (2.3) the fractional de- rivative (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) ∂ ∂ α ∂ ∂p t k b a A n p x dnn n a b / ( / ( )) ( ) /= − [ ] ++ +∫ 1 1 β∂ ∂k p x/+ 2 2 P u t LT p x tx u( , ) ( , ),= ∂ ∂n np x t x( , ) / = ∂ ∂n x n x w p w t w dw/ ( ) ( ) ( , ) /−[ ] − [ ]−∫1 1 0 Γ LTx u m n n m m m m m f x u u F u u f u f f u f , ( ) ( ) ( ) ( ) ( ) ( ) ( ) .... ( ) ( ) + − − − − [ ] = − +[ − − − − ] 1 2 1 1 0 0 0 0' ∂ ∂ α β α ∂ ∂ β ∂ P u t t k us us b a us u u P u t k p t u p t x us us b a us k up t p t b a b a ( , ) / (( ) ( ) ) / ( ) log( ) ( , ) ( , ) ( , ) / (( ) ( ) ) / ( ) log ( , ) ( , ) − −[ − + ] = = − +[ ] − − − + + + − + + 1 1 2 2 1 1 10 0 0 0 // ∂x[ ] vV u v P u k us us b a us u V u v LT u t u t k p t u p t x us us b a b a t v b a ( , ) ( , ) (( ) ( ) ) / log ( , ) ( , ) ( , ) ( , ) ( , ) / (( ) ( ) ) / ( ) log( , − − −[ −( ) ( ) + ] = = − +[ ] − − + + − + + 0 0 0 1 1 2 1 1 1 α β α ∂ ∂ Ξ Ξ usus k up t p t x ) ( , ) ( , ) / + − +[ ]β ∂ ∂0 0 V u v P u v ku k us us b a us k us us b a us p t u p t x v k us us b a b a t v b ( , ) ( , ) / ( (( ) ( ) ) / ( ) log ) ((( ) ( ) ) / ( ) log ( , ) ( , ) / / (( ) ( ) , = − + −[ − ] − −{ − ) + ][ } − − − − + + − + 0 0 0 2 1 1 1 1 1 β α α ∂ ∂ α LT ++ + + − −[ ] + − + +[ − − − −[ ] a t v b a b a us ku k up t p t x v k us us b a us ku ) / ( ) log ( , ) ( , ) / / (( ) ( ) ) / ( ) log , β β ∂ ∂ α β 2 1 1 2 0 0LT 227 Diffusion with space memory modelled with distributed order space fractional differential equations (2.8) where the asterisk indicates convolution and its subscript indicates that the convolution is made with respect to the variable t; performing the LTt v, −1 we fi nd (2.8') As shown in the Appendix the LTt 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 re- spectively. 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 P u t V u v P u v k us us b a us ku k us us b a us p t u p t t v t v b a u b a ( , ) ( , ) ( , ) / ( (( ) ( ) ) / ( ) log (( ) ( ) ) / ( ) log ) ( , ) ( , ) / , , = = = − −[{ − − ]} − −{ − + − − + + + + − LT LT 1 1 1 1 1 1 1 0 1 0 0 α β α ∂ ∂xx v ku us us b a us ku t t v b a u [ ]} − −[{ − − ]} + − + + * / ( (( ) ( )) / ( ) log ) ,LT 1 1 11 α β − +[ ] − −[{ − − ]} − + + β ∂ ∂ α β k up t p t x v k us us b a us ku t v t v b a LT LT , , ( , ) ( , ) / * / (( ) ( ) ) / ( ) log 0 0 11 1 1 2 P u t V u v P u tk us us b a us tku t v b a ( , ) ( , ) ( , ) exp ( (( ) ( ) ) / ( ) log( )) ,= = − − +{ } + − + + LT 1 1 1 2 0 α β − − − +[ + ] ( ) − ( )( )[ + + − + + α ∂ ∂ α k us us b a us p t u p t x tk us us b a b a ((( ) ( ) ) / ( ) log ) ( , ) ( , ) / exp / 1 1 1 1 1 0 0 ( ) log , , / exp / log . b a us ku k up t p t x tk us us b a us ku b a − + ] − ( ) + ∂ ( ) ∂[ ] ( ) − ( )( ) −( ) −[ ]+ + β β α β 2 1 1 2 0 0 b a us kp t us us b a us tk u tku us us x u b a b a −( ) } − −{[ − } + −{ − + + + + log ( , ) (( ) ( ) ) / ( ) log exp (( ) ( ) ) / ,α β α 0 1 1 1 2 1 1 LT b a us k p t x u us us b a us tk u tku us us t x u b a b a −( ) }] − − −{ }[ + ( ) − ( )( ){ − − + + − − log ( , ) / * (( ) ( ) ) / ( ) log exp / , α ∂ ∂ β α 0 1 1 1 1 2 LT b a us kp t LT u tk us us b a us t ku x u b a −( ) }] − { − − +[ ]} + − + + log ( , ) exp (( ) ( ) ) / ( ) log ,β α β 0 1 1 1 2 − [{ + − − + ]} − + + β ∂ ∂ α β k p t x tk us us b a us ku t t x u b a ( , ) / * exp (( ) ( ) ) / ( ) log ,0 1 1 1 2 LT p x t LT P u t p x LT tk u u u b a u x u x x u b a , , , exp / log , ,( ) = ( ) = ( ) + −( ) −( ) )({ } + − ∗ − + + 1 1 2 1 1 0 β α − − −[ ] +[{ + ] + −( )([ − + + − + + LT k u u b a u p t u p t x tk u u u x u b a b a , (( ) / ( ) log ) ( , ) ( , ) / exp / 1 1 1 1 2 1 1 0 0 α ∂ ∂ β α ( ) log ) , , / exp / log b a u k up t p t x tk us us b a us ku x u b a − ]} − ( ) + ∂ ( ) ∂[ ]{ ( ) − ( )( ) −( ) +[ ]}+ + LT , -1 β α β 0 0 1 1 2 p x t P u t p x tk u tk us us x u x x u b a ( , ) ( , ) ( , ) exp (( ) ( )) / , , = = + −{ − ∗ − + + LT LT 1 1 2 1 1 0 β α 228 Michele Caputo 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 Æ 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 Æ 0 and fi nd that ua in (2.9') substitutes the term u u b a ub a1 1− −−( ) −( )/ log whose limit, when b Æ a, is ua . With the extreme value theorem applied to eq. (2.8') it is verifi ed that p(0,0) = 0. 3. 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 (3.1) Using the expression 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) Once more one may verify the computations as- suming b = a + e , computing the limit of (3.2) when e Æ 0 and fi nding that the exponent of the exponential becomes tk rs aaα π( )( )+1 cos and the argument of the sine becomes tk rs aaα π( ) +1 sin 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 (3.2') With the extreme value theorem applied to eq. (2.7) it is verifi ed that p(x,0) = C. p x t LT P u t p x tk us us b a us x u x x u b a ( , ) ( , ) ( , ) exp (( ) ( ) ) / ( ) log . , , = = − −{ }[ ] − ∗ − + + 1 1 1 1 0 LT α p x t p x rx tk rs nb rs a rs rs b rs a b a rs tk rs b rs x b a b a b ( , ) ( ( , ) / ) exp( )sin (( ) sin ( ) sin log ( ) cos ( ) cos ) / ( )(log ) exp cos = = − [{( − ) − +( − ] − + } −( ) + ( ∗ + ∞ + + + + ∫ 0 1 0 1 1 1 2 2 1 π α π π π π π π α π )) )([[ + − ( ) − ( ) ]( −( ) +( )]) + + + 1 1 1 2 2 a b a a rs rs b rs a b a rs dr cos log sin sin / log . π π π π π p x t C rx r tk rs nb rs a rs rs b rs a b a rs b a b a ( , ) ( / ) (( exp( )) / )sin (( ) sin ( ) sin ) log ( ) cos ( ) cos ) / ( )(log ) = = − − [{( − − +( − ] − + } + ∞ + + + ∫ π α π π π π π π 1 1 0 1 1 1 2 2 exp ( ( ) cos ( ) cos ) log (( ) sin ( ) sin ) / ( )(log ) . tk rs b rs a rs rs b rs a b a rs dr b a b a α π π π π π π − + +[[ − − ] − + ]) + + + + 1 1 1 1 2 2 229 Diffusion with space memory modelled with distributed order space fractional differential equations 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 LTx u b atk us us us. exp ( ( ) ( ) ) / (log ) − + +− + ]{[1 1 1 2 2α π and fi nally that LTx u tku, exp( ) − [ ]1 2β is known since it is the solution of the classic initial value case without memory. The effect of the memory is obviously re- presented 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. 4. 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 LTx 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 (4.2) and introduce from eqs. (A.7) and (A.8) of the Appendix, the values of l, m, c, d (4.3) we obtain p x t P u t p x LT tk u tk us us b a us u x x u b a ( , ) ( , ) ( , ) exp (( ) ( ) ) / ( ) log , , = = + −{[ − }] − ∗ − + + LTx 1 1 2 1 1 0 β α p x t p x tk us us us LT tku x x u b a x x u ( , ) ( , ) exp (( ) ( ) ) / (log ) exp( ) , , = = −{[ + [ ][ ∗ − + + ∗ − 0 1 1 1 2 2 1 2 LT α π β p x t kp t us us b a us tk us us b a us t x u b a b a ( , ) ( , ) (( ) ( ) ) / ( ) log exp (( ) ( ) ) / ( ) log . ,= −{[ − } −{ − }] ∗ − + + + + α α 0 1 1 1 1 1 LT LTx u u t rx c l d m d dr , ( , ) ( / ) exp( ) exp ( sin cos ) − ∞ = = − +∫ 1 0 1 Ψ π c tk rs b rs a rs rs b rs a b a rs b a b a = − + +[ − − ] − + + + + + α π π π π π π ( ( ) cos ( ) cos ) log (( ) sin ( ) sin ) ) / ( )(log ) 1 1 1 1 2 2 d rs b rs a rs rs b rs a b a rs b a b a = − +[ − − ] − + + + + + (( ) sin ( ) sin ) log (( ) cos ( ) cos ) ) / ( )(log ) 1 1 1 1 2 2 π π π π π π l rs b rs a rs rs b rs a b a rs b a b a = − + +[ − − ] − + + + + + ( ( ) cos ( ) cos ) log (( ) sin ( ) sin ) ) / ( )(log ) 1 1 1 1 2 2 π π π π π π m rs b rs a rs rs b rs a b a rs b a b a = − +[ − − ] − + + + + + (( ) sin ( ) sin ) log (( ) cos ( ) cos ) ) / ( )(log ) 1 1 1 1 2 2 π π π π π π exp ( ( ) cos ( ) cos ) log (( ) sin ( ) sin ) ) / tk rs b rs a rs rs b rs a b a b a α π π π π π − + +[{ − − ] + + + + 1 1 1 1 LTx u u t rx, ( , ) ( / ) exp( ) − ∞ = −   ∫ 1 0 1Ψ π 230 Michele Caputo (4.4) Assuming p(0,t) = d (t), we obtain from eqs. (4.2) the Green function of the problem. 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 p3(x,t) the solution (3.2'), the solution of the boundary value problem is C – p3(x,t). 5. 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. b a rs rs b rs a rs rs b rs a b a b a −( ) +( )} − +[{[ ( ) − − ] + + + + log ( ( ) cos cos ) log (( ) sin ( ) sin ) ) / 2 2 1 1 1 1 π π π π π π b a rs kt rs b rs a rs rs b rs a b a rs b a b a −( ) +( )} +[{ − − + − − + } + + + + + log sin (( ) sin ( ) sin ) log (( ) cos ( ) cos ) ) / ( )(log ) 2 2 1 1 1 1 2 2 π α π π π π π π + − +[{ − − ] − + } ( ) +([{ + + + + + (( ) sin ( ) sin ) log (( ) cos ( ) cos ) ) / ( )(log ) cos sin rs b rs a rs rs b rs a b a rs kt rs b b a b a b 1 1 1 1 2 2 1 π π π π π π α π − − + − ] − + }]) + + + ( ) sin ) log (( ) cos ( ) cos ) ) / ( )(log ) . rs a rs rs b rs a b a rs dr a b a 1 1 1 2 2 π π π π π Concerning the initial value problem the solution (3.2') may be compared with that of Caputo and Plastino (2002) who used the con- stitutive 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+asinap 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 fi - nancial phenomena (Caputo and Kolari, 2001; Caputo, 2001b). p x t C xr r ktr r a ktr r a dr a a ( , ) ( / ) (( exp( )) / )sin( ( cos )) exp( ( cos ) = = − − + ∞ − + ∫π α π α π β 1 2 0 2 1 231 Diffusion with space memory modelled with distributed order space fractional differential equations 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 r1 when r1 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 LTu 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 LTx 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). Glossary a, b limits of the distribution of the fractional order derivatives. k (s-2 m2) ratio of the fl uid pressure to the fl uid density (see eq. (2.2)). n fractional order of differentiation (see eqs. (2.3) and (2.4)) (dimensionless). p(x,t) (kg s-2 m-1) fl uid pressure. q(x,t) (kg s-1 m-2) fl uid mass fl ow rate in the porous medium. r1, R radius of the inner and outer circles, respectively, of the integration path of eq. (A.1) shown in fi g. A.1. s (m) variable normalizing the dimensions. t (s) time. x (m) distance from the boundary plane. a (s mn ) coeffi cient of the Darcy’s law modifi ed (see eq. (2.3)). ak (s-1m2+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. The difference of the two, in all terms appearing in eq. (2.9'), may be reduced to the form (A.2) ∞ ∞ = − − − = − − = − ] =    ∫∫ 0 0 1 2( / ) exp( ) ( , ) exp( ) ( , )Ψ Ψi rx r dr rx r drπ θ π θ π LT− − ∞ + ∞ = − =∫1 1 2( ( )) ( / ) exp( ) ( )Ψ Ψu i xu u du i i π ξ ξ = − = − = −[1 2( / ) exp( ) ( , ) ( , )Ψ Ψi rx r r drπ θ π θ π ]] =    ∞ ∫ 0 232 Michele Caputo where (l + im)exp(c + id) and (l – im)exp(c – id) represent the values of the integrand of (2.9') when the argument of m in the complex plane is q = p and q = - p respectively. The integral along the path of fi g. 1 gives then (A.3)p x t P u t rx l d m d c drx u( , ) ( , ) ( / ) exp( )( sin cos ) exp( ) .,= = − + − ∞ ∫LT 1 0 1 π = − +{ } = ∞ ∫( / ) exp( ) (exp )sin (exp ) cos1 2 0 πi rx il c d im c d dr = − +[ ] ∞ ∫( / ) exp( ) exp sin cos1 0 π rx c l d m d dr = − − −[ ] +( / ) exp( exp (exp exp( ) exp exp1 2πi rx l c id id im c idid id dr+ −[ ]{ } = ∞ ∫ exp( ) 0 Fig. A.1. Path of integration in the complex plane LTu x, −1 in eq. (2.9'). b and iu are the real and imaginary parts respectively. = − + + − − −[ ] = = − + − −[ ] + + + −[ ]{ } = ∞ ∞ ∫ ∫ ( / ) exp( ( ) exp( ) ( ) exp( ) ( / ) exp( (exp( ) exp( ) exp( ) exp( ) 1 2 1 2 0 π π i rx l im c id l im c id dr i rx l c id c id im c id c id dr o (A.2) 233 Diffusion with space memory modelled with distributed order space fractional differential equations 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.7) and d tk rs b r b rs a a rs b a rs b a = − +[ − − + ] − + + + α π π π π π π π ( ) (sin log cos ) ( ) ( cos sin log ) / ( )(log ). 1 1 2 2 c tk rs b rs a rs rs b rs a b a rs b a b a = − + +[ − − ] − + + + + + α π π π π π π ( ( ) cos ( ) cos ) log (( ) sin ( ) sin ) ) / ( )(log ) 1 1 1 1 2 2 d tk rs b rs a rs rs b rs a b a rs b a b a = − +[ − − ] − + + + + + α π π π π π π (( ) sin ( ) sin ) log (( ) cos ( ) cos ) ) / ( )(log ) 1 1 1 1 2 2 l rs b rs a rs rs b rs a b a rs b a b a = − + +[ − − ] − + + + + + ( ( ) cos ( ) cos ) log (( ) sin ( ) sin ) ) / ( )(log ) 1 1 1 1 2 2 π π π π π π exp (( ) ( ) ) / ( ) logtk us us b a usb aα 1 1+ +− −{ } c tk rs b rs a rs rs b rs a b a rs b a b a = − + +[ − − ] − + + + + + α π π π π π π ( ( ) cos ( ) cos ) log (( ) sin ( ) sin )) / ( )(log ) 1 1 1 1 2 2 m rs b rs a rs rs b rs a b a rs b a b a = − +[ − − ] − + + + + + (( ) sin ( ) sin ) log (( ) cos ( ) cos ) ) / ( )(log ). 1 1 1 1 2 2 π π π π π π l r rs b rs a rs rs b rs a b a rs b a b a = − − + +[ − − ] − + − + + + + 1 1 1 1 1 2 2 ( ( ) cos ( ) cos ) log (( ) sin ( ) sin ) ) / ( )(log ) π π π π π π m r rs b rs a rs rs b rs a b a rs b a b a = − − +[ − − ] − + − + + + + 1 1 1 1 1 2 2 (( ) sin ( ) sin ) log (( ) cos ( ) cos ) ) / ( )(log ). π π π π π π (A.6) 234 Michele Caputo REFERENCES BAGLEY, R.L. (1991): The thermorheologically complex material, Int. J. Eng. Sci., 29 (7), 797-806. BAGLEY, R.L. and P.J TORVIK (2000a): On the existence of the order domain and the solution of distributed order equations - Part I, Int. J. Appl. Math., 2 (7), 865-882. BAGLEY, R.L. and P.J. TORVIK (2000b): On the existence of the order domain and the solution of distributed order equations - Part II, Int. J. Appl. Math., 2 (8), 965-987. BELFIORE, L. and M. CAPUTO (2000): The experimental set-valued index of refraction of dielectric and anelastic media, Ann. Geofi s., 43 (2), 207-216. CAPUTO, M. (1967): Linear models of dissipation whose Q is almost frequency independent-II, Geophys. J. R. Astron. Soc., 13, 529-539. CAPUTO, M. (1969): Elasticità e Dissipazione (Zanichelli Publisher, Bologna). CAPUTO, M. (1989): The rheology of an anelastic medium studied by means of the observation of the splitting of its eigenfrequencies, J. Acoust. Soc. Am., 85 (5), 1984-1987. CAPUTO, M. (1994): The Riemann sheet solutions of linear anelasticity, Ann. Mat. Pura Appl., IV (CLXVI), 335- 342. CAPUTO, M. (1995): Mean fractional-order-derivatives differential equations and fi lters, Ann. Univ. Ferrara Sci. Mat., 41, 73-83. CAPUTO, M. (1996): Modern rheology and dielectric induction: multivalued index of refraction, splitting of eigenvalues and fatigue, Ann. Geofi s., 39 (5), 941-966. CAPUTO, M. (2000): Models of flux in porous media with memory, Water Resour. Res., 36 (3), 693-705. CAPUTO, M. (2001a): Distributed order differential equations modelling dielectric induction and diffusion, in Fractional Calculus and Applied Analysis, 4 (4), 421-442.; also in Atti Accad. Sci. Ferrara, 77, 161-184 (1999-2000). CAPUTO, M. (2001b): Economy equilibrium equations with memory, in Nuovi Progressi nella Fisica Matematica dall’Eredità di Dario Graffi , Accademia Nazionale dei Lincei, 177, 191-220. CAPUTO, M. and J. KOLARI (2001): An analytical model of the Fisher equation with memory function, Alternative Perspectives of Finance and Accounting, 1, 1-14. CAPUTO, M. and W. PLASTINO (2002): Diffusion with space memory, in Geodesy, the Challenge of the 3rd Millenium, edited by E. GRAFAREND, F. KRUMM and V. SCHWARZE (Springer Verlag, Berlin-Heidelberg), 429-435. CHRISTAKOS, G., D.T. HRISTOPULOS and D.T. MILLER (1995): Stochastic diagrammatic analysis of ground water flow in heterogeneous porous media, Water Resour. Res., 31 (7), 1687-1703. CISOTTI, V. (1911): L’ereditarietà lineare ed i fenomeni dispersivi, Nuovo Cimento, 6 (2), 234-244. CUSHMAN, J.H and M. MORONI (2001): Statistical mechanics with 3D-PTV experiments in the study of anomalous dispersion: part I, Theory, Phys. Fluids, 13 (1),75-80. FLORES, E. and T.J. OSLER (1999): The tautochrone under arbitrary potentials using fractional derivatives, Am. J. Phys., 67 (8), 718-721. GRAFFI, D. (1936): Sopra alcuni fenomeni ereditari del- l’elettrologia: I, II, Rend. Ist. Sci. Lett. Arti, 2 (69), 128-181. HU, X. and J. CUSHMAN (1991): Non equilibrium statistical mechanical derivation of a non-local Darcy’s law for unsaturated/saturated fl ow, Stochastic Hydrol. Hydraul., 8, 109-116. JACQUELIN, J. (1991): Synthese de circuits électriques equivalents a partir de mesures d’impédences complexes, in 5ème Forum sur les Impedences Électrochimiques, 287-295. KABALA, Z.J. and G. SPOSITO (1991): A stochastic model of reactive solute transport with time-varying velocity in a heterogeneous aquifer, Water Resour. Res., 27 (3), 341-350. KIRYAKOVA, V. (1994): Generalised fractional calculus and applications, Pitman Research Notes in Mathematics, n. 301 (Longman Harlow). KÖRNIG, H. and G. MÜLLER (1989): Rheological model and interpretation of postglacial uplift, Geophys. J. R. Astron. Soc., 98, 245-253. LE MEHAUTE, A. and G. CRÉPY (1983): Introduction to transfer motion in fractal media: the geometry of kinetics, Solid State Ionic, 9-10, 17-30. LUCHKO, Y. and R. GORENFLO (1998): The initial value problem for some fractional differential equations with the Caputo derivatives, Preprint n. A8-98, Fachber. Math. Informatik, Serie A, Mathematik (Freie Univesität Berlin). MAINARDI, F. (1996): Fractional relaxation-oscillation and fractional diffusion-wave phenomena, Chaos, Solitons and Fractals, 7, 1461-1477. MAINARDI, F., YU. LUCHKO and G. PAGNIN (2001): The fundamental solution of space-time fractional diffusion equation, Fractional Calculus and Applied Analysis, 4 (2), 153-192. OLDHAM, K.B. and J. SPANIER (1974): The Fractional Calculus (Academic Press). PODLUBNY, I. (1994a): Analytical solution of linear differential equations of the fractional order, in Proceedings of the 12th IMACS (Intern. Assoc. for Mathematics and Comp. in Simulation) World Congress, Atlanta, U.S.A. PODLUBNY, I. (1994b): Numerical solution of initial value problems for ordinary fractional-order differential equations, in Proceedings IMACS (Intern. Assoc. for Mathematics and Comp. in Simulation) World Congress, Atlanta, U.S.A. PODLUBNY, I. (1999): Fractional Differential Equations (Academic Press, New York). ROSS, B. (1975): A brief history and exposition of the fundamental theory of fractional calculus, Lect. Notes Math. (Springer Verlag), vol. 457. WYSS, W. (1986): Fractional diffusion equation, J. Math. Phys., 27, 2782-2785.