# Similarity Solution for Fractional Diffusion Equation

Authors: Jun-Sheng Duan; Ai-Ping Guo; Wen-Zai Yun
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

Fractional diffusion equation in fractal media is an integropartial differential equation parametrized by fractal Hausdorff dimension and anomalous diffusion exponent. In this paper, the similarity solution of the fractional diffusion equation was considered. Through the invariants of the group of scaling transformations we derived the integro-ordinary differential equation for the similarity variable. Then by virtue of Mellin transform, the probability density functionp(r,t), which is just the fundamental solution of the fractional diffusion equation, was expressed in terms of Fox functions.


---

## Body

## 1. Introduction

Standard diffusion ind-dimensional space, where d is a positive integer, is a process described by Gaussian distribution. A main feature of the process is the linear relation between the mean square displacement and time; namely, 〈r2(t)〉∝t. Some anomalous diffusion phenomena that take place in impure media, biological tissues, and porous media can be simulated by the diffusion model in fractals [1–6]. In recent years, the fractal theory has been developed rapidly, and it was found to be closely related to the anomalous diffusion phenomena [3–12].In fractal media, the geometric obstacles existing on all length scales slow down the particle motion in a random walk. The mean square displacement behaves as [2]
(1)R2≡〈r2(t)〉∝t2/dw,
where dw(>2) is the anomalous diffusion exponent. The numerical simulation found that on a large class of fractal structures the general form of the probability density function p(r,t) that the walker is at distance r at time t from its starting point at time t=0 obeys asymptotically a non-Gaussian shape of the form [2, 3]
(2)p(r,t)~t-df/dwexp⁡[-const.×(rR)u],rR≫1,
where u=dw/(dw-1) and df is the fractal Hausdorff dimension.In order to simulate the diffusion phenomena in fractal media, some scholars have introduced fractional diffusion equations [4, 5, 11–13]. In this paper, we consider the fractional diffusion equation [5, 13]:
(3)∂γp(r,t)∂tγ=1rds-1∂∂r(rds-1∂p∂r),r>0,t>0,
where γ=2/dw, ds=2df/dw is the spectral dimension of the fractal, and the fractional time derivative on the left hand side of (3) is defined as the convolution integral [14–20]:
(4)∂γp(r,t)∂tγ=∂∂t∫0t‍(t-τ)-γΓ(1-γ)p(r,τ)dτ,0<γ<1,
where Γ(·) is Euler's gamma function. In the limit case, dw→2 and df→d, (3) reduces to the standard d-dimensional diffusion equation.The fractional calculus has been applied to many fields in science and engineering, such as viscoelasticity, anomalous diffusion, biology, chemistry, and control theory [5, 11–13, 15, 19–22]. Researches on the fractional differential equations attract much attention [15, 23–28]. For linear fractional differential equations, the integral transforms, including the Laplace, Fourier, and Mellin transforms, are usually used to obtain analytic solutions.In this paper using the similarity method [29] we solve (3) with the following initial and boundary conditions and the conservation condition:
(5)p(r,0)=0,r>0,p(∞,t)=0,t>0,ω(df)∫0∞‍p(r,t)rdf-1dr=1,
where ω(df) is a constant, which is defined as
(6)ω(df)=2πdf/2Γ(df/2).We note that the probability density functionp(r,t) is just the fundamental solution of the fractional diffusion equation. The similarity method was used by Gorenflo et al. [30], Wyss [31], and Buckwar and Luchko [32] for solving problems of time fractional partial differential equations in one-dimensional case.


## 2. Derivation of Similarity Solution

First we determine a symmetric group of scaling transformations(7)Tα:r=αr-,t=αht-,p=αlp-,
where α>0 is a parameter and h, l are constants to be determined. Applying the group of scaling transformations (7), the fractional derivative is converted as follows:
(8)∂γp(r,t)∂tγ=αl∂∂t∫0t‍(t-τ)-γΓ(1-γ)p-(r-,α-hτ)dτ=αl-hγ∂∂t-∫0t-‍(t--τ′)-γΓ(1-γ)p-(r-,τ′)dτ′=αl-hγ∂γp-(r-,t-)∂t-γ,
where τ′=α-hτ. Hence the problem (3)–(5) is invariant under the group (7) if and only if
(9)h=2γ,l=-df.
So the symmetric group of scaling transformations is determined:
(10)Tα:r=αr-,t=α2/γt-,p=α-dfp-.Eliminating the parameterα leads to two invariants:
(11)rt-γ/2=r-t--γ/2,rdfp=r-dfp-.
We denote the two invariants of the group of the scaling transformation Tα as
(12)η=rt-γ/2,F=rdfp.Next we use the transformation(13)p(r,t)=r-dfF(η),η=rt-γ/2
to determine the equations for the similarity solution of the problem (3)–(5). Calculating derivative we have
(14)∂p∂r=-dfr-df-1F(η)+r-dft-γ/2F′(η),(15)1rds-1∂∂r(rds-1∂p∂r)=r-dft-γF′′(η)+Δ1r-df-1t-γ/2F′(η)+Δ2r-df-2F(η),
where
(16)Δ1=ds-1-2df,Δ2=df(df-ds+2).
For the left hand side of (3), we introduce the new integral variable
(17)ξ=rτ-γ/2,
we obtain p(r,τ)=r-dfF(ξ), and
(18)∂γp(r,t)∂tγ=∂∂t∫0t‍(t-τ)-γΓ(1-γ)p(r,τ)dτ=-r2/γ-df-1t-γ/2-1ddη∫η+∞‍[(η/ξ)-2/γ-1]-γΓ(1-γ)aaaaaaaaaaaaaaaaaaaaaaa×F(ξ)ξ1-2/γdξ.
Letting
(19)g(w)={(w-2/γ-1)-γΓ(1-γ),0<w<1,0,w>1,
we rewrite (18) as
(20)∂γp(r,t)∂tγ=-r2/γ-df-1t-2/γ-1ddη∫0+∞‍g(ηξ)F(ξ)ξ1-2/γdξ.From (15) and (20), we obtain the integro-ordinary differential equation for the similarity variables:
(21)-ddη∫0+∞‍g(ηξ)F(ξ)ξ-2/γ+1dξ=η-2/γ+1F′′(η)+Δ1η-2/γF′(η)+Δ2η-2/γ-1F(η).
The conditions (5) are converted to
(22)F(+∞)=0,ω(df)∫0+∞‍F(η)η-1dη=1.Considering the integration in (21), we use Mellin transforms for the new problem (21) and (22). The Mellin transform of function f(x) is defined as [33]
(23)f^(s)=ℳ[f(x),s]=∫0+∞‍f(x)xs-1dx.Applying Mellin transform with respect toη to both sides of (21), we get
(24)(s-1)g^(s-1)F^(s-2γ+1)=((s-2γ)(s-2γ-1)-Δ1(s-2γ-1)+Δ2)×F^(s-2γ-1).
Calculating integrations we obtain Mellin transform of the function g(w):
(25)g^(s)=Γ(γ+γs/2)sΓ(γs/2).
Inserting (25) into (24) and then replacing s by s+2/γ+1 we obtain the difference equation for the function F^(s):
(26)Γ(γ(s/2+1/γ+1))Γ(γ(s/2+1/γ))F^(s+2)=(s2+s-Δ1s+Δ2)F^(s).In order to solve the difference equation, we introduces=2q and F^(2q)=T(q), and rewrite (26) into
(27)T(q+1)T(q)=4(q+df2)(q+1+df-ds2)Γ(γq+1)Γ(γq+γ+1).
A particular solution of (27) is
(28)T(q)=C4qΓ(q+df/2)Γ(q+1+(df-ds)/2)Γ(γq+1),
where C is an arbitrary constant. For the solution of (27), we can multiply T(q) by any function Y(q) which satisfies Y(q+1)/Y(q)=1.We notice thatF^(s) is a Mellin transform defined only in some strip 0≤σ1<Re(s)<σ2 from the conditions (22). So (26) is valid only in the overlap of the two strips σ1<Re(s)<σ2 and σ1<Re(s+2)<σ2, and there is no such overlap unless σ1+2<σ2. Thus Y(q) cannot have poles; otherwise, it would have a row of poles separated exactly by one unit. In addition, Y(q) cannot grow faster than |q| as Im⁡(q)→∞ in the inversion strip; otherwise the inversion integral would diverge. Thus Y(q) is a bounded entire function and equals a constant by Liouville's theorem.Therefore,T(q) has only the form of (28) and we have
(29)F^(s)=T(s2)=C2sΓ(df/2+s/2)Γ(1+(df-ds)/2+s/2)Γ(1+γs/2).It follows from (22) that F^(0)=1/ω(df). Thus we have
(30)C=1ω(df)Γ(df/2)Γ(1+(df-ds)/2).The inverse Mellin transform of (29) is
(31)F(η)=C2πi∫c-i∞c+i∞‍2sΓ(df/2+s/2)Γ(1+(df-ds)/2+s/2)Γ(1+γs/2)aaaaaaaaaaaa×η-sds.
Replacing s by -s and using the definition of Fox functions we obtain [34, 35]
(32)F(η)=CH1,22,0(η2|(df/2,1/2),(1+(df-ds)/2,1/2)(1,γ/2)).
Inserting the expressions into (13) and using properties of Fox functions, we obtain the probability density function in terms of the Fox function:
(33)p(r,t)=dw2-dft-df/dwω(df)Γ(df/2)Γ(1+df/2-df/dw)×H1,22,0(rdw2dwt|(0,dw/2),(1-df/dw,dw/2)(1-df/dw,1)).For a large class of fractal structures, the spectral dimension [2] satisfies ds<2; that is, df<dw. So the Fox function in (33) can be expanded into a series by using residue theorem on the simple poles:
(34)Pa={2kdw∣k=0,1,…}∪{2dw(1-dfdw+k)∣k=0,1,…}.
The series representation for the probability density is calculated to be
(35)p(r,t)=21-dft-df/dwω(df)Γ(df/2)Γ(1+df/2-df/dw)×∑k=0∞(-1)kk!×[Γ(1-df/dw-k)Γ(1-df/dw-2k/dw)(rdw2dwt)2k/dwaaaaaaaaaaaaaaaaaa+Γ(df/dw-1-k)Γ(1-df/dw-2/dw(1-df/dw+k))aaaaaaaaaaaaaaaaaa×(rdw2dwt)2/dw(1-df/dw+k)].


## 3. Discussions and Conclusions

In the limit case,dw→2 and df→d, (3) reduces to the d-dimensional standard diffusion equation, and the probability density (35) is simplified to the Gaussian distribution:
(36)p(r,t)=1(4πt)d/2exp⁡(-r24t).In Figures1 and 2, we plot the curves of p(r,0.25) versus r and p(r,0.5) versus r, respectively, for df=1 and different values of dw. In Figures 3 and 4, we plot the curves of p(r,0.25) versus r and p(r,0.5) versus r, respectively, for df=1.5 and different values of dw. The figures display that, as the anomalous diffusion exponent dw increases, the peak value of the probability density function p(r,t) at r=0 decreases. In addition, as the fractal Hausdorff dimension df increases from 1 to 1.5, the peak value of p(r,t) at r=0 decreases.Figure 1
Curves ofp(r,0.25) versus r for df=1 and for dw=2 (solid line), dw=2.5 (dot line), dw=3 (dash line), and dw=3.5 (dot-dash line).Figure 2
Curves ofp(r,0.5) versus r for df=1 and for dw=2 (solid line), dw=2.5 (dot line), dw=3 (dash line), and dw=3.5 (dot-dash line).Figure 3
Curves ofp(r,0.25) versus r for df=1.5 and for dw=2 (solid line), dw=2.5 (dot line), dw=3 (dash line), and dw=3.5 (dot-dash line).Figure 4
Curves ofp(r,0.5) versus r for df=1.5 and for dw=2 (solid line), dw=2.5 (dot line), dw=3 (dash line), and dw=3.5 (dot-dash line).Compared with the similarity method for classic partial differential equations, the similarity method for fractional diffusion equation involves the similarity integral variableξ=rτ-γ/2, and the reduction equation is an integro-ordinary differential equation for the similarity solution. The obtained probability density p(r,t) is just the fundamental solution of the fractional diffusion equation.



---

*Source: 548126-2014-03-12.xml*